text:
# Plotting examples
# =================
# ## Plotters
#  All plot functions in GridVisualize.jl have a `Plotter` keyword argument
#  which defaults to `nothing`.  This allows to pass a module as plotting backend
#  without creating a dependency. Fully supported are `PyPlot`, `GLMakie` and `PlutoVista`.
#  `WGLMakie` and `CairoMakie` work in principle but in the moment don't deliver
#  all necessary functionality. For `Plots` we miss the possibility to work with triangle meshes (this is under development, though)
#
# Also supported is [`VTKView`](https://github.com/j-fu/VTKView.jl)  which is exeprimental and works only on linux.
#
#
# ## Grid plots
# Here, we define some sample grids for plotting purposes.

using ExtendableGrids
using GridVisualize

# ### 1D grids
function grid1d(;n=50)
    X=collect(0:1/n:1)
    g=simplexgrid(X)
end

function plotting_grid1d(;Plotter=default_plotter(), kwargs...)
    gridplot(grid1d(); Plotter=Plotter, resolution=(500,200),kwargs...)
end
# ![](plotting_grid1d.svg)

# ### 2D grids

function grid2d(;n=20)
    X=collect(0:1/n:1)
    g=simplexgrid(X,X)
end

function plotting_grid2d(;Plotter=default_plotter(), kwargs...)
    gridplot(grid2d(); Plotter=Plotter, kwargs...)
end
# ![](plotting_grid2d.svg)

# ### 3D grids
# The kwargs `xplane`, `yplane` and `zplane` which allow to control
# cutplanes which peel off some elements from the grid in 3d and allow to
# explore the inner triangulation.
#
# For Makie and VTKView, the cutplane values can be controlled interactively.

function grid3d(;n=15)
    X=collect(0:1/n:1)
    g=simplexgrid(X,X,X)
end

function plotting_grid3d(;Plotter=default_plotter(), kwargs...)
    gridplot(grid3d(); Plotter=Plotter, kwargs...)
end
# ![](plotting_grid3d.svg)

# ## Function plots

# ### Function on 1D grid
function func1d(;n=50)
    g=grid1d(n=n)
    g,map(x->sinpi(2*x[1]),g)
end

function plotting_func1d(;Plotter=default_plotter(), kwargs...)
    g,f=func1d()
    scalarplot(g,f; Plotter=Plotter,resolution=(500,300),kwargs...)
end
# ![](plotting_func1d.svg)

# ### Function on 2D grid
function func2d(;n=20)
    g=grid2d(n=n)
    g,map((x,y)->sinpi(2*x)*sinpi(3.5*y),g)
end

function plotting_func2d(;Plotter=default_plotter(), kwargs...)
    g,f=func2d()
    scalarplot(g,f; Plotter=Plotter,levels=10, kwargs...)
end
# ![](plotting_func2d.svg)

# ### Function on 3D grid
#
# Plotting a function then goes as follows:
# `xplane`, `yplane` and `zplane` now define cut planes where
# the function projection is plotted as a heatmap.
# The additional `flevel` keyword argument allows
# to control an isolevel.
#
# For Makie and VTKView, the cutplane values and the flevel can be controlled interactively.

function func3d(;n=15)
    g=grid3d(n=n)
    g, map((x,y,z)->sinpi(2*x)*sinpi(3.5*y)*sinpi(1.5*z),g)
end

function plotting_func3d(;Plotter=default_plotter(), kwargs...)
    g,f=func3d()
    scalarplot(g,f; Plotter=Plotter, zplane=0.49,xplane=0.49,flevel=0.25, kwargs...)
end
# ![](plotting_func3d.svg)

# ## Vector and stream plots
# ### 2D vector
function vec2d(;n=20)
    g=grid2d(n=n)
    g,vcat(map((x,y)->sinpi(2*x)*sinpi(3.5*y),g)',map((x,y)->cospi(2*x)*cospi(3.5*y),g)')
end

function plotting_vec2d(;Plotter=default_plotter(), kwargs...)
    g,f=vec2d()
    vectorplot(g,f; Plotter=Plotter,kwargs...)
end
# ![](plotting_vec2d.svg)

# ### 2D stream
# Stream plots are currently only available with PyPlot.
function plotting_stream2d(;Plotter=default_plotter(), kwargs...)
    g,f=vec2d()
    streamplot(g,f; Plotter=Plotter,spacing=0.05,kwargs...)
end
# ![](plotting_stream2d.svg)

# ## Multiscene plots
# We can combine multiple plots into one scene according to
# some layout grid given by the layout parameter.
#
# This is not currently supported by the PlutVista backend.
#
# The ',' key for GLMakie and the '*' key for VTKView allow to
# switch between gallery view (default) and focused view of only
# one subscene.
function plotting_multiscene(;Plotter=default_plotter(),resolution=(1000,500))
    p=GridVisualizer(;Plotter=Plotter,layout=(2,4),clear=true,resolution=resolution)
    gridplot!(p[1,1],grid1d(), title="1D grid")
    scalarplot!(p[2,1],grid1d(), sin, title="1D grid function", label="sin",markershape=:diamond,color=:red,legend=:rb)
    scalarplot!(p[2,1],grid1d(), cos, title="1D grid function", label="cos",linestyle=:dash,markershape=:none,color=:green,clear=false)
    gridplot!(p[1,2],grid2d(),title="2D grid")
    scalarplot!(p[2,2],func2d()...,colormap=:bamako,title="2D grid function")
    gridplot!(p[1,3],grid3d(),zplane=0.49,title="3D grid")
    scalarplot!(p[2,3],func3d()...,zplane=0.49,flevel=0.5,colormap=:bamako, title="3D grid function")
    vectorplot!(p[1,4],vec2d()...,title="2D quiver")
    GridVisualize.streamplot!(p[2,4],vec2d()...,title="2D stream")
    reveal(p)
end
# ![](plotting_multiscene.svg)


text:
#=
lm-block.jl
block operations for linear maps
2019-06-15, Jeff Fessler, University of Michigan
=#

export block_lm

using LinearMaps
using LinearAlgebra: UniformScaling, I
using Test: @test

#const BlockVector{T} = Vector{Union{LinearMap,AbstractMatrix{T}}}
if !@isdefined(BlockVector)
	BlockVector{T} = Vector{<:Any} # include I etc.
end

"""
`ob = block_lm(blocks, ...)`

Construct `LinearMap` object composed of blocks.
such as `block_diag(A_1, A_2, ..., A_M)`.
Each block can be anything like a `LinearMap` such as an `AbstractMatrix`.

See `block_lm(:test)` for example usage.

in
* `blocks::BlockVector`	array of the blocks

* option
* `how::Symbol` options:
  + `:col` for `[A_1; A_2; ...; A_M]`
  + `:diag` for block diagonal (default)
  + `:kron` for `kron(eye(Mkron), blocks[1])`
  + `:row` for `[A1, A2, ..., A_M]`
  + `:sum` for `A1 + A2 + ... + A_M`
* `T::DataType` default `promote_type(eltype.(blocks)...)`
* `Mkron::Int' required for `:kron` type
* `tomo::Bool` special support for tomography-type objects; default `false`
	(not implemented)

out
* `ob` `LinearMap`
"""
function block_lm(
		blocks::BlockVector;
		T::DataType = promote_type(eltype.(blocks)...),
		how::Symbol = :diag,
		tomo::Bool = false,
		Mkron::Int = 0,
	)

	if how === :col
		 return block_lm_col(blocks, T, tomo)
	elseif how === :diag
		return block_lm_diag(blocks, T)
	elseif how === :kron
		return block_lm_kron(blocks, T, Mkron)
	elseif how === :row
		return block_lm_row(blocks, T)
	elseif how === :sum
		return block_lm_sum(blocks, T)
	else
		throw(ArgumentError("unknown block type $how"))
	end
end

"""
`block_lm_col()`
"""
function block_lm_col(blocks::BlockVector, T::DataType, tomo::Bool)

	MM = length(blocks)
	dims = zeros(Int, MM, 2)
	for mm=1:MM
		B = blocks[mm]
		if isa(B, UniformScaling)
			throw("I not yet supported")
		end
		dims[mm,:] .= size(B)
		dims[mm,2] != dims[1,2] && throw("all blocks must have same #cols for :col")
	end

	# start/end indices for selecting parts of x and y
	istart = cumsum([1; dims[1:end-1,1]])
	iend = istart + dims[:,1] .- 1

	dim = [sum(dims[:,1]), dims[1,2]]

#	if tomo # prep for mat2cell later
#		for mm=1:MM
#			odim = blocks[mm].odim
#			tomo_ndim = length(odim)
#			tomo_dim = tomo_ndim # insist on last dimension
#			if mm==1
#				for id = 1:tomo_ndim
#					if id != tomo_dim
#						mat2cell_arg[id] = odim[id]
#					end
#				end
#			end
#			mat2cell_arg[tomo_dim][mm] = odim(tomo_dim)
#		end
#	end

	return LinearMap{T}(
			x -> vcat([blocks[mm] * x for mm=1:MM]...),
			y -> sum([blocks[mm]' * (@view y[istart[mm]:iend[mm]]) for mm=1:MM]),
			dim[1], dim[2])
end

"""
#
# block_lm_col_gram()
#
function [T, reuse] = block_lm_col_gram(ob, W, reuse, varargin)

blocks = ob.blocks
T = cell(size(blocks))
for (mm,A) in enumerate(blocks)
	if isnumeric(A)
		if isempty(W)
			T[mm] = A' * A
		else
			warn 'todo: this may not work, need piece of W!'
			T[mm] = A' * W * A
		end
	else
		if isempty(W)
			T[mm] = build_gram(A, [], reuse, varargin{:})
		else
			if isvar('W.blocks[mm]')
				T[mm] = build_gram(A, W.blocks[mm], ...
					reuse, varargin{:})
			else
				fail('block_lm_col_gram needs block diag W')
			end
		end
	end
end
T = block_lm(T, how=:sum)
"""

"""
`block_lm_row()`
trick: just reuse col via transpose!
cannot use hcat(blocks...) because there might be a mix of Matrix/LinearMaps
"""
function block_lm_row(blocks::BlockVector, T::DataType)
	tblocks = [b' for b in blocks] # trick: transpose
	return block_lm(tblocks, how=:col)' # trick: transpose
end

"""
`block_lm_diag()`
"""
function block_lm_diag(blocks::BlockVector, T::DataType)

	MM = length(blocks)
#	@show dims = vcat(map(size, blocks))
	dims = zeros(Int, MM, 2)
	for mm=1:MM
		dims[mm,:] .= size(blocks[mm])
	end

	# start/end indices for selecting parts of x and y
	istart = cumsum([1; dims[1:end-1,1]])
	iend = istart + dims[:,1] .- 1
	jstart = cumsum([1; dims[1:end-1,2]])
	jend = jstart + dims[:,2] .- 1

	dim = sum(dims, dims=1)

	return LinearMap{T}(
		x -> vcat([blocks[mm] * x[jstart[mm]:jend[mm]] for mm=1:MM]...),
		y -> vcat([blocks[mm]' * y[istart[mm]:iend[mm]] for mm=1:MM]...),
		dim[1], dim[2])
end

"""
#
# block_lm_diag_mtimes_block()
# caution: it is quite unclear whether these are useful!
#
function y = block_lm_diag_mtimes_block(arg, is_transpose, x, istart, nblock)
if is_transpose
	y = block_lm_diag_block_back(arg, x, istart, nblock)
else
	y = block_lm_diag_block_forw(arg, x, istart, nblock)
end

#
# block_lm_diag_block_forw()
#
function y = block_lm_diag_block_forw(arg, x, istart, nblock)

if nrow(x) ~= dim(2)
	throw('x size=%d vs dim(2)=%d', nrow(x), dim(2))
end
y = []
for mm=istart:nblock:length(blocks)
	t = blocks[mm] * x([jstart(mm):jend(mm)], :)
	y = [y; t]
end

#
# block_lm_diag_block_back()
#
function x = block_lm_diag_block_back(arg, y, istart, nblock)

if nrow(y) ~= dim(1), error 'bad y size', end
x = []
for mm=istart:nblock:length(blocks)
	t = blocks[mm]' * y([istart(mm):iend(mm)], :)
	x = [x; t]
end

#
# block_lm_diag_gram()
#
function [T, reuse] = block_lm_diag_gram(ob, W, reuse, varargin)

blocks = ob.blocks
T = cell(size(blocks))
for (mm,A) in enumerate(blocks)
	if isnumeric(A)
		if isempty(W)
			T[mm] = A' * A
		else
			T[mm] = A' * W * A
		end
	else
		T[mm] = build_gram(A, W, reuse, varargin{:})
	end
end
T = block_lm(T, 'type', 'diag')
"""

#
# block_lm_kron()
#
function block_lm_kron(blocks::BlockVector, T::DataType, Mkron::Int)

	length(blocks) != 1 && throw(ArgumentError("kron expects one block"))

	dims = repeat(collect(size(blocks[1]))', Mkron, 1)
	# start/end indices for selecting parts of x and y
	istart = cumsum([1; dims[1:end-1,1]])
	iend = istart + dims[:,1] .- 1
	jstart = cumsum([1; dims[1:end-1,2]])
	jend = jstart + dims[:,2] .- 1
	dim = sum(dims,dims=1)

	return LinearMap{T}(
			x -> vcat([blocks[1] * x[jstart[mm]:jend[mm]] for mm=1:Mkron]...),
			y -> vcat([blocks[1]' * y[istart[mm]:iend[mm]] for mm=1:Mkron]...),
			dim[1], dim[2])
end

"""
#
# block_lm_kron_mtimes_block(): y = A{ib} * x
#
function y = block_lm_kron_mtimes_block(arg, is_transpose, x, iblock, nblock)

bl1 = blocks{1}; % base block, already put through Gblock
warn 'todo: size check not done'

if ~is_transpose % forw
#	if nrow(x) ~= dim(2)
#		fail('x size=%d vs dim(2)=%d', nrow(x), dim(2))
#	end
	y = []
	for mm=1:Mkron
		t = bl1{iblock} * x([jstart(mm):jend(mm)], :)
		y = [y; t]
	end

else % back
	y = x
#	if nrow(y) ~= dim(1), error 'bad y size', end
	x = []
	for mm=1:Mkron
		tmp = [istart(mm):iend(mm)]; % i list (if all data)
		% todo: we need a certain subset of that list
		fail('todo: kron subset backprojector not done')
		t = bl1{iblock}' * y(tmp, :)
		x = [x; t]
	end
	y = x
end
"""

#
# block_lm_sum()
#
function block_lm_sum(blocks::BlockVector, T::DataType)
	dim = size(blocks[1])
	MM = length(blocks)
	for mm=1:MM
		!isequal(size(blocks[mm]), dim) &&
			throw(DimensionMismatch("blocks must have same size for :sum"))
	end

	return LinearMap{T}(
			x -> sum([blocks[mm] * x for mm=1:MM]),
			y -> sum([blocks[mm]' * y for mm=1:MM]),
			dim[1], dim[2])
end

"""
#
# block_lm_free()
#
function block_lm_free(arg)
if chat
	printm 'freeing block_lm object static memory'
end
for mm=1:length(blocks)
	try
		free(blocks[mm])
	catch
	end
end
"""

"""
#
# block_lm_update()
#
function out = block_lm_update(ob, varargin)
# todo: figure out how to update, e.g., new_zmap(s), ...
"""

"""
`block_lm(:test)`
self test
"""
function block_lm(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))

	# test :col
	A = ones(4,3)
	B = LinearMap(x -> A*x, y -> A'*y, 4, 3)
#	C = I # todo later
	C = [2A; 3A]
	blocks = [A, B, C]
	Tc = block_lm(blocks, how=:col)
	Tc * ones(3)
	Tc' * ones(16)
	@test Matrix(Tc)' == Matrix(Tc')

	# test :diag
	A = ones(3,2)
	M = ones(4,3)
	B = LinearMap(x -> M*x, y -> M'*y, 4, 3)
#	C = I # todo later
	C = [2M; 3M]
	blocks = [A, B, C]
	Td = block_lm(blocks, how=:diag)
	Td * ones(8)
	Td' * ones(15)
	@test Matrix(Td)' == Matrix(Td')

	# test :kron
	Tk = block_lm([ones(3,2)], how=:kron, Mkron=4)
	Tk * ones(4*2)
	Tk' * ones(4*3)
	@test Matrix(Tk)' == Matrix(Tk')

	# test :row
	A = ones(4,3)
	B = LinearMap(x -> A*x, y -> A'*y, 4, 3)
#	C = I # todo later
	C = [2A 3A]
	blocks = [A, B, C]
	Tr = block_lm(blocks, how=:row)
	Tr * ones(12)
	Tr' * ones(4)
	@test Matrix(Tr)' == Matrix(Tr')

	# test :sum
	A = ones(4,3)
	B = LinearMap(x -> A*x, y -> A'*y, 4, 3)
#	C = I # todo later
	C = 2A
	blocks = [A, B, C]
	Ts = block_lm(blocks, how=:sum)
	Ts * ones(3)
	Ts' * ones(4)
	@test Matrix(Ts)' == Matrix(Ts')

	@test_throws ArgumentError block_lm(blocks, how=:bad)
	@test_throws String block_lm([I], how=:col)

	true
end

# block_lm(:test)

ext: jl


text:
function Tile(arg0::GeodeticPoint, arg1::GeodeticPoint, arg2::GeodeticPoint, arg3::GeodeticPoint)
    return Tile((GeodeticPoint, GeodeticPoint, GeodeticPoint, GeodeticPoint), arg0, arg1, arg2, arg3)
end

function get_center(obj::Tile)
    return jcall(obj, "getCenter", GeodeticPoint, ())
end

function get_interpolated_point(obj::Tile, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "getInterpolatedPoint", GeodeticPoint, (jdouble, jdouble), arg0, arg1)
end

function get_vertices(obj::Tile)
    return jcall(obj, "getVertices", Vector{GeodeticPoint}, ())
end


text:
module prof

using RustModels
using ProfileView
import Markov.rsm

include("../test/toymodels.jl")

srand(425)
models=classicaltoys()
dirm=models[1]
theta0=[0.6,0.8]
coef!(dirm,theta0)
data=rand(dirm,100,100)

Profile.clear()
@profile mle(models[2],data);
ProfileView.view()

end #module


text:
import IntervalArithmetic: mid, ±, isinterior, entireinterval, where_bisect,
    diam, bisect

export bisect

diam(a::AbstractVector) = maximum(diam.(a))

function bisect(X::AbstractVector, α=where_bisect)
    i = argmax(diam.(X))  # find longest side
    return bisect(X, i, α)
end

function bisect(X::AbstractVector, i::Integer, α=where_bisect)
    x1, x2 = bisect(X[i], α)
    X1 = setindex(X, x1, i)
    X2 = setindex(X, x2, i)
    return (X1, X2)
end

contains_root(region) = all(0 .∈ region)
isinterior(a::AbstractVector, b::AbstractVector) = all(isinterior.(a, b))

strict_isinterior(a, b) = isinterior(a, b) && !any(a .=== entireinterval.(a))

# these definitions break the concept of Vectors as sets of their elements
# (rather than Cartesian boxes)
_isempty(a) = any(isempty.(a))


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/11/2014

module AbstractCASCoordImpl

export AbstractCASCoord

abstract AbstractCASCoord

end

module AbstractCASCoordInterfaces

export
    setRecord,
    getRecord,
    getAll

function setRecord()
end

function getRecord()
end

function getAll()
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
using Documenter, Schedulers

makedocs(sitename="Schedulers", modules=[Schedulers])

deploydocs(
    repo = "git@github.com:ChevronETC/Schedulers.jl.git"
)


text:
function parse_line(line)
  low = 0
  high = 0
  char = ""
  passw = ""
  start = 1
  idx = 1
  len = length(line)
  while idx <= len
    if line[idx + 1] == '-'
      break
    end

    idx += 1
  end

  low = parse(Int32, line[start:idx])
  idx += 2
  start = idx
  while idx <= len
    if line[idx + 1] == ' '
      break
    end

    idx += 1
  end

  high = parse(Int32, line[start:idx])
  char = line[idx + 2]
  passw = line[idx + 5:end]

  low, high, char, passw
end

function load(input_path)
  map(parse_line, readlines(input_path))
end

function valid_pw_1(spec)
  low, high, char, passw = spec
  count = 0
  idx = 1
  while (count <= high + 1) && (idx <= length(passw))
    if passw[idx] == char
      count += 1
    end

    idx += 1
  end

  low <= count <= high
end

function solve_1(data)
  length(filter(valid_pw_1, data))
end

function valid_pw_2(spec)
  low, high, char, passw = spec
  xor(passw[low] == char, passw[high] == char)
end

function solve_2(data)
  length(filter(valid_pw_2, data))
end

function solve(data)
  println(@time solve_1(data))
  println(@time solve_2(data))
end

data = load("input.txt")
solve(data)

max_stars_repo_name: wbthomason/advent

max_issues_repo_name: wbthomason/advent

max_forks_repo_name: wbthomason/advent


text:
# Mon 31 Dec 2012 01:41:39 PM EST
#
# NeHe Tut 19 - Make some colored stars and play w/ alpha blending a bit more
#
# Q - quit
# L - turn lights on/off
# F - change texture filter (linear, nearest, mipmap)
# PageUp/Down - move camera closer/further away
# Up/Down - increase/decrease x-rotation speed
# Left/Right - increase/decrease y-rotation speed
# Home/End - increase/decrease the amount of "friction" on particles
# Enter - turn rainbow effect on/off
# Space - step through colors for rainbow effect

# load necessary GLUT/OpenGL routines

global OpenGLver="1.0"
using OpenGL
using GLUT

# initialize variables

global window

global tilt    = 90.0
global zoom    = -15.0
global spin    = 0.0

global twinkle = false

global tex     = Array(Uint32,1) # generating 1 texture

width          = 640
height         = 480

colors = [1.0  0.5  0.5;
          1.0  0.75 0.5;
          1.0  1.0  0.5;
          0.75 1.0  0.5;
          0.5  1.0  0.5;
          0.5  1.0  0.75;
          0.5  1.0  1.0;
          0.5  0.75 1.0;
          0.5  0.5  1.0;
          0.75 0.5  1.0;
          1.0  0.5  1.0;
          1.0  0.5  0.75]

global rainbow       = true
global slowdown      = 2.0
global xspeed        = 0
global yspeed        = 0
global zoom          = -40.0
global color         = 1
global delay         = 0

global MAX_PARTICLES = 1000

global keystate = [false false false false false]

type particle
    active::Bool
    life::Float64
    fade::Float64
    red::Float64
    green::Float64
    blue::Float64
    xPos::Float64
    yPos::Float64
    zPos::Float64
    xSpeed::Float64
    ySpeed::Float64
    zSpeed::Float64
    xGrav::Float64
    yGrav::Float64
    zGrav::Float64
end

global particles = [particle(true,            # Julia doesn't like it when you try to initialize an empty array of
                      1.0,                    # a composite type and try to fill it afterwards, so we
                      rand(1:100)/1000+0.003,  # start with a 1-element vector and tack on values
                      colors[1,1],
                      colors[1,2],
                      colors[1,3],
                      0.0,
                      0.0,
                      0.0,
                      (rand(1:50)-26.0)*10.0,
                      (rand(1:50)-26.0)*10.0,
                      (rand(1:50)-26.0)*10.0,
                      0.0,
                      -0.8,
                      0.0)]

for loop = 2:MAX_PARTICLES
    active = true
    life   = 1.0
    fade   = rand(1:100)/1000+0.003
    red    = colors[loop%size(colors,1)+1,1]
    green  = colors[loop%size(colors,1)+1,2]
    blue   = colors[loop%size(colors,1)+1,3]
    xPos   = 0.0
    yPos   = 0.0
    zPos   = 0.0
    xSpeed = (rand(1:50)-26.0)*10.0
    ySpeed = (rand(1:50)-26.0)*10.0
    zSpeed = (rand(1:50)-26.0)*10.0
    xGrav  = 0.0
    yGrav  = -0.8
    zGrav  = 0.0
    particles = push!(particles, particle(active,life,fade,red,green,blue,xPos,yPos,zPos,xSpeed,ySpeed,zSpeed,xGrav,yGrav,zGrav))
end

# load textures from images

function LoadGLTextures()
    global tex

    img, w, h = glimread(expanduser("~/.julia/GLUT/Examples/NeHe/tut19/Particle.bmp"))

    glGenTextures(1,tex)
    glBindTexture(GL_TEXTURE_2D,tex[1])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)
end

# function to init OpenGL context

function initGL(w::Integer,h::Integer)
    global tex

    glViewport(0,0,w,h)
    LoadGLTextures()
    glClearColor(0.0, 0.0, 0.0, 0.0)
    glClearDepth(1.0)
    glDisable(GL_DEPTH_TEST)
    glShadeModel(GL_SMOOTH)
    glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST)
    glHint(GL_POINT_SMOOTH_HINT, GL_NICEST)

    # enable texture mapping and alpha blending
    glEnable(GL_TEXTURE_2D)
    glEnable(GL_BLEND)
    glBlendFunc(GL_SRC_ALPHA, GL_ONE)
    glColor(1.0, 1.0, 1.0, 0.5)

    glBindTexture(GL_TEXTURE_2D,tex[1])

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)
end

# prepare Julia equivalents of C callbacks that are typically used in GLUT code

function ReSizeGLScene(w::Int32,h::Int32)
    if h == 0
        h = 1
    end

    glViewport(0,0,w,h)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    return nothing
end

_ReSizeGLScene = cfunction(ReSizeGLScene, Void, (Int32, Int32))

function DrawGLScene()
    global zoom
    global particles
    global MAX_PARTICLES
    global tex
    global rainbow
    global slowdown
    global xspeed
    global yspeed
    global zoom
    global color
    global delay
    global keystate

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    for loop = 1:MAX_PARTICLES
        if particles[loop].active
            x = particles[loop].xPos
            y = particles[loop].yPos
            z = particles[loop].zPos+zoom
            glColor(particles[loop].red, particles[loop].green, particles[loop].blue, particles[loop].life)
            glBegin(GL_TRIANGLE_STRIP)
                glTexCoord(1, 1)
                glVertex(x+0.5, y+0.5, z)
                glTexCoord(0, 1)
                glVertex(x-0.5, y+0.5, z)
                glTexCoord(1, 0)
                glVertex(x+0.5, y-0.5, z)
                glTexCoord(0, 0)
                glVertex(x-0.5, y-0.5, z)
            glEnd()
            particles[loop].xPos   += particles[loop].xSpeed/(slowdown*1000)
            particles[loop].yPos   += particles[loop].ySpeed/(slowdown*1000)
            particles[loop].zPos   += particles[loop].zSpeed/(slowdown*1000)
            particles[loop].xSpeed += particles[loop].xGrav
            particles[loop].ySpeed += particles[loop].yGrav
            particles[loop].zSpeed += particles[loop].zGrav
            particles[loop].life   -= particles[loop].fade
            if particles[loop].life < 0.0
                particles[loop].life   = 1.0
                particles[loop].fade   = (rand(1:100))/1000+0.003
                particles[loop].xPos   = 0.0
                particles[loop].yPos   = 0.0
                particles[loop].zPos   = 0.0
                particles[loop].xSpeed = xspeed+rand(1:60)-32.0
                particles[loop].ySpeed = yspeed+rand(1:60)-30.0
                particles[loop].zSpeed = rand(1:60)-30.0
                particles[loop].red    = colors[color,1]
                particles[loop].green  = colors[color,2]
                particles[loop].blue   = colors[color,3]
            end
            if keystate[1] == true && particles[loop].yGrav < 1.5
                particles[loop].yGrav +=0.01
                keystate[1] = false
            elseif keystate[2] == true && particles[loop].yGrav > -1.5
                particles[loop].yGrav -=0.01
                keystate[2] = false
            elseif keystate[3] == true && particles[loop].xGrav < 1.5
                particles[loop].xGrav +=0.01
                keystate[3] = false
            elseif keystate[4] == true && particles[loop].xGrav > -1.5
                particles[loop].xGrav -=0.01
                keystate[4] = false
            elseif keystate[5] == true
                particles[loop].xPos   = 0.0
                particles[loop].yPos   = 0.0
                particles[loop].zPos   = 0.0
                particles[loop].xSpeed = (rand(1:52)-26.0)*10
                particles[loop].ySpeed = (rand(1:50)-25.0)*10
                particles[loop].zSpeed = (rand(1:50)-25.0)*10
                keystate[5] = false
            end
        end
    end

    delay +=1

    glutSwapBuffers()

    return nothing
end

_DrawGLScene = cfunction(DrawGLScene, Void, ())

function keyPressed(the_key::Char,x::Int32,y::Int32)
    global delay
    global color
    global keystate
    global rainbow

    if the_key == int('q')
        glutDestroyWindow(window)
    elseif the_key == int(' ')
        delay = 0
        color +=1
        if color > 12
            color = 1
        end
    elseif the_key == int('2')
        keystate[1] = true
    elseif the_key == int('8')
        keystate[2] = true
    elseif the_key == int('4')
        keystate[3] = true
    elseif the_key == int('6')
        keystate[4] = true
    elseif the_key == int('\t')
        keystate[5] = true
    elseif the_key == int('\n')
        rainbow = (rainbow ? false : true)
    end

    return nothing # keyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_keyPressed = cfunction(keyPressed, Void, (Char, Int32, Int32))

function specialKeyPressed(the_key::Int32,x::Int32,y::Int32)
    global zoom
    global xspeed
    global yspeed
    global slowdown
    global keystate

    if the_key == GLUT_KEY_PAGE_UP
        zoom +=0.1
    elseif the_key == GLUT_KEY_PAGE_DOWN
        zoom -=0.1
    elseif the_key == GLUT_KEY_UP
        yspeed +=1.0
    elseif the_key == GLUT_KEY_DOWN
        yspeed -=1.0
    elseif the_key == GLUT_KEY_LEFT
        xspeed +=1.0
    elseif the_key == GLUT_KEY_RIGHT
        xspeed -=1.0
    elseif the_key == GLUT_KEY_HOME
        slowdown -=0.01
    elseif the_key == GLUT_KEY_END
        slowdown +=0.01
    end

    return nothing # specialKeyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_specialKeyPressed = cfunction(specialKeyPressed, Void, (Int32, Int32, Int32))

# run GLUT routines

glutInit()
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH)
glutInitWindowSize(width, height)
glutInitWindowPosition(0, 0)

window = glutCreateWindow("NeHe Tut 19")

glutDisplayFunc(_DrawGLScene)
glutFullScreen()

glutIdleFunc(_DrawGLScene)
glutReshapeFunc(_ReSizeGLScene)
glutKeyboardFunc(_keyPressed)
glutSpecialFunc(_specialKeyPressed)

initGL(width, height)

glutMainLoop()


text:
"""
    unsafe_createfeaturedefn(name::AbstractString)

Create a new feature definition object to hold field definitions.

The `FeatureDefn` maintains a reference count, but this starts at zero, and
should normally be incremented by the owner.
"""
unsafe_createfeaturedefn(name::AbstractString)::FeatureDefn =
    FeatureDefn(GDAL.ogr_fd_create(name))

"""
    reference(featuredefn::FeatureDefn)

Increments the reference count in the FeatureDefn by one.

The count is used to track the number of `Feature`s referencing this definition.

### Returns
The updated reference count.
"""
reference(featuredefn::FeatureDefn)::Integer =
    GDAL.ogr_fd_reference(featuredefn.ptr)

"""
    dereference(featuredefn::FeatureDefn)

Decrements the reference count by one, and returns the updated count.
"""
dereference(featuredefn::FeatureDefn)::Integer =
    GDAL.ogr_fd_dereference(featuredefn.ptr)

"""
    nreference(featuredefn::AbstractFeatureDefn)

Fetch the current reference count.
"""
nreference(featuredefn::AbstractFeatureDefn)::Integer =
    GDAL.ogr_fd_getreferencecount(featuredefn.ptr)

"Destroy a feature definition object and release all memory associated with it"
function destroy(featuredefn::FeatureDefn)::Nothing
    GDAL.ogr_fd_destroy(featuredefn.ptr)
    featuredefn.ptr = C_NULL
    return nothing
end

"Destroy a feature definition view"
function destroy(featuredefn::IFeatureDefnView)::Nothing
    featuredefn.ptr = C_NULL
    return nothing
end

"""
    release(featuredefn::FeatureDefn)

Drop a reference, and destroy if unreferenced.
"""
function release(featuredefn::FeatureDefn)::Nothing
    GDAL.ogr_fd_release(featuredefn.ptr)
    return nothing
end

"""
    getname(featuredefn::AbstractFeatureDefn)

Get name of the OGRFeatureDefn passed as an argument.
"""
getname(featuredefn::AbstractFeatureDefn)::String =
    GDAL.ogr_fd_getname(featuredefn.ptr)

"""
    nfield(featuredefn::AbstractFeatureDefn)

Fetch number of fields on the passed feature definition.
"""
nfield(featuredefn::AbstractFeatureDefn)::Integer =
    GDAL.ogr_fd_getfieldcount(featuredefn.ptr)

"""
    getfielddefn(featuredefn::FeatureDefn, i::Integer)

Fetch field definition of the passed feature definition.

### Parameters
* `featuredefn`: the feature definition to get the field definition from.
* `i`:  index of the field to fetch, between `0` and `nfield(featuredefn)-1`.

### Returns
an handle to an internal field definition object or NULL if invalid index. This
object should not be modified or freed by the application.
"""
getfielddefn(featuredefn::FeatureDefn, i::Integer)::FieldDefn =
    FieldDefn(GDAL.ogr_fd_getfielddefn(featuredefn.ptr, i))

getfielddefn(featuredefn::IFeatureDefnView, i::Integer)::IFieldDefnView =
    IFieldDefnView(GDAL.ogr_fd_getfielddefn(featuredefn.ptr, i))

"""
    findfieldindex(featuredefn::AbstractFeatureDefn,
        name::Union{AbstractString, Symbol})

Find field by name.

### Returns
the field index, or -1 if no match found.

### Remarks
This uses the OGRFeatureDefn::GetFieldIndex() method.
"""
function findfieldindex(
    featuredefn::AbstractFeatureDefn,
    name::Union{AbstractString,Symbol},
)::Integer
    return GDAL.ogr_fd_getfieldindex(featuredefn.ptr, name)
end

"""
    addfielddefn!(featuredefn::FeatureDefn, fielddefn::FieldDefn)

Add a new field definition to the passed feature definition.

To add a new field definition to a layer definition, do not use this function
directly, but use OGR_L_CreateField() instead.

This function should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn. The OGRFieldDefn passed in is copied,
and remains the responsibility of the caller.
"""
function addfielddefn!(
    featuredefn::FeatureDefn,
    fielddefn::FieldDefn,
)::FeatureDefn
    GDAL.ogr_fd_addfielddefn(featuredefn.ptr, fielddefn.ptr)
    return featuredefn
end

"""
    deletefielddefn!(featuredefn::FeatureDefn, i::Integer)

Delete an existing field definition.

To delete an existing field definition from a layer definition, do not use this
function directly, but use `OGR_L_DeleteField()` instead.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn.
"""
function deletefielddefn!(featuredefn::FeatureDefn, i::Integer)::FeatureDefn
    result = GDAL.ogr_fd_deletefielddefn(featuredefn.ptr, i)
    @ogrerr result "Failed to delete field $i in the feature definition"
    return featuredefn
end

"""
    reorderfielddefns!(featuredefn::FeatureDefn, indices::Vector{Cint})

Reorder the field definitions in the array of the feature definition.

To reorder the field definitions in a layer definition, do not use this function
directly, but use `OGR_L_ReorderFields()` instead.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn.

### Parameters
* **fd**: handle to the feature definition.
* **indices**: an array of `GetFieldCount()` elements which is a permutation of
    `[0, GetFieldCount()-1]`. `indices` is such that, for each field definition
    at position `i` after reordering, its position before reordering was
    `indices[i]`.
"""
function reorderfielddefns!(
    featuredefn::FeatureDefn,
    indices::Vector{Cint},
)::FeatureDefn
    result = GDAL.ogr_fd_reorderfielddefns(featuredefn.ptr, indices)
    @ogrerr result "Failed to reorder $indices in the feature definition"
    return featuredefn
end

"""
    getgeomtype(featuredefn::AbstractFeatureDefn)

Fetch the geometry base type of the passed feature definition.

For layers without any geometry field, this method returns `wkbNone`.

This returns the same result as `OGR_FD_GetGeomType(OGR_L_GetLayerDefn(hLayer))`
but for a few drivers, calling `OGR_L_GetGeomType()` directly can avoid lengthy
layer definition initialization.

For layers with multiple geometry fields, this method only returns the geometry
type of the first geometry column. For other columns, use
    `OGR_GFld_GetType(OGR_FD_GetGeomFieldDefn(OGR_L_GetLayerDefn(hLayer), i))`.
"""
getgeomtype(featuredefn::AbstractFeatureDefn)::OGRwkbGeometryType =
    GDAL.ogr_fd_getgeomtype(featuredefn.ptr)

"""
    setgeomtype!(featuredefn::FeatureDefn, etype::OGRwkbGeometryType)

Assign the base geometry type for the passed layer (same as the fd).

All geometry objects using this type must be of the defined type or a derived
type. The default upon creation is `wkbUnknown` which allows for any geometry
type. The geometry type should generally not be changed after any OGRFeatures
have been created against this definition.
"""
function setgeomtype!(
    featuredefn::FeatureDefn,
    etype::OGRwkbGeometryType,
)::FeatureDefn
    GDAL.ogr_fd_setgeomtype(featuredefn.ptr, etype)
    return featuredefn
end

"""
    isgeomignored(featuredefn::AbstractFeatureDefn)

Determine whether the geometry can be omitted when fetching features.
"""
isgeomignored(featuredefn::AbstractFeatureDefn)::Bool =
    Bool(GDAL.ogr_fd_isgeometryignored(featuredefn.ptr))

"""
    setgeomignored!(featuredefn::FeatureDefn, ignore::Bool)

Set whether the geometry can be omitted when fetching features.
"""
function setgeomignored!(featuredefn::FeatureDefn, ignore::Bool)::FeatureDefn
    GDAL.ogr_fd_setgeometryignored(featuredefn.ptr, ignore)
    return featuredefn
end

"""
    isstyleignored(featuredefn::AbstractFeatureDefn)

Determine whether the style can be omitted when fetching features.
"""
isstyleignored(featuredefn::AbstractFeatureDefn)::Bool =
    Bool(GDAL.ogr_fd_isstyleignored(featuredefn.ptr))

"""
    setstyleignored!(featuredefn::FeatureDefn, ignore::Bool)

Set whether the style can be omitted when fetching features.
"""
function setstyleignored!(featuredefn::FeatureDefn, ignore::Bool)::FeatureDefn
    GDAL.ogr_fd_setstyleignored(featuredefn.ptr, ignore)
    return featuredefn
end

"""
    ngeom(featuredefn::AbstractFeatureDefn)

Fetch number of geometry fields on the passed feature definition.
"""
ngeom(featuredefn::AbstractFeatureDefn)::Integer =
    GDAL.ogr_fd_getgeomfieldcount(featuredefn.ptr)

"""
    getgeomdefn(featuredefn::FeatureDefn, i::Integer = 0)

Fetch geometry field definition of the passed feature definition.

### Parameters
* `i`  geometry field to fetch, between `0` (default) and `ngeomfield(fd)-1`.

### Returns
an internal field definition object or `NULL` if invalid index. This object
should not be modified or freed by the application.
"""
getgeomdefn(featuredefn::FeatureDefn, i::Integer = 0)::GeomFieldDefn =
    GeomFieldDefn(GDAL.ogr_fd_getgeomfielddefn(featuredefn.ptr, i))

getgeomdefn(featuredefn::IFeatureDefnView, i::Integer = 0)::IGeomFieldDefnView =
    IGeomFieldDefnView(GDAL.ogr_fd_getgeomfielddefn(featuredefn.ptr, i))

"""
    findgeomindex(featuredefn::AbstractFeatureDefn, name::AbstractString = "")

Find geometry field by name.

The geometry field index of the first geometry field matching the passed field
name (case insensitively) is returned.

### Returns
the geometry field index, or -1 if no match found.
"""
function findgeomindex(
    featuredefn::AbstractFeatureDefn,
    name::AbstractString = "",
)::Integer
    return GDAL.ogr_fd_getgeomfieldindex(featuredefn.ptr, name)
end

"""
    addgeomdefn!(featuredefn::FeatureDefn, geomfielddefn::AbstractGeomFieldDefn)

Add a new field definition to the passed feature definition.

To add a new geometry field definition to a layer definition, do not use this
function directly, but use OGRLayer::CreateGeomField() instead.

This method does an internal copy of the passed geometry field definition,
unless bCopy is set to `false` (in which case it takes ownership of the field
definition.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn.
"""
function addgeomdefn!(
    featuredefn::FeatureDefn,
    geomfielddefn::AbstractGeomFieldDefn,
)::FeatureDefn
    # `geomfielddefn` is copied, and remains the responsibility of the caller.
    GDAL.ogr_fd_addgeomfielddefn(featuredefn.ptr, geomfielddefn.ptr)
    return featuredefn
end

"""
    deletegeomdefn!(featuredefn::FeatureDefn, i::Integer)

Delete an existing geometry field definition.

To delete an existing field definition from a layer definition, do not use this
function directly, but use OGRLayer::DeleteGeomField() instead.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn.
"""
function deletegeomdefn!(featuredefn::FeatureDefn, i::Integer)::FeatureDefn
    result = GDAL.ogr_fd_deletegeomfielddefn(featuredefn.ptr, i)
    @ogrerr result "Failed to delete geom field $i in the feature definition"
    return featuredefn
end

"""
    issame(featuredefn1::AbstractFeatureDefn, featuredefn2::AbstractFeatureDefn)

Test if the feature definition is identical to the other one.
"""
function issame(
    featuredefn1::AbstractFeatureDefn,
    featuredefn2::AbstractFeatureDefn,
)::Bool
    return Bool(GDAL.ogr_fd_issame(featuredefn1.ptr, featuredefn2.ptr))
end

"""
    unsafe_createfeature(featuredefn::AbstractFeatureDefn)

Returns the new feature object with null fields and no geometry

Note that the OGRFeature will increment the reference count of it's defining
OGRFeatureDefn. Destruction of the OGRFeatureDefn before destruction of all
OGRFeatures that depend on it is likely to result in a crash.

Starting with GDAL 2.1, returns NULL in case out of memory situation.
"""
function unsafe_createfeature(featuredefn::AbstractFeatureDefn)::Feature
    return Feature(GDAL.ogr_f_create(featuredefn.ptr))
end

"""
    getfeaturedefn(feature::Feature)

Fetch feature definition.
"""
getfeaturedefn(feature::Feature)::IFeatureDefnView =
    IFeatureDefnView(GDAL.ogr_f_getdefnref(feature.ptr))


text:
using Iris
using Test
using LinearAlgebra
using SparseArrays

@testset "LU Factorizations" begin
    A = I + sprand(1000,1000,.2)
    y = rand(1000)
    alu = @test_nowarn lu(A,USolver)
    @test norm(A*(alu\y)-y) ≤ 1e-12
    # alu = @test_skip lu(A,PSolver)
    # @test_skip norm(A*(alu\y)-y) ≤ 1e-12
    # alu = @test_skip lu(A,MSolver)
    # @test_skip norm(A*(alu\y)-y) ≤ 1e-12
end


text:
"""
    ContinuousMarkovChain{T1,T2}

Continuous-time Markov chain with finite state space.

# Fields
- `π0::T1`: initial state distribution
- `Q::T2`: rate matrix.
"""
struct ContinuousMarkovChain{T1<:AbstractVector{<:Real},T2<:AbstractMatrix{<:Real}} <:
       AbstractMarkovChain
    π0::T1
    Q::T2
end

MeasureTheory.sampletype(::ContinuousMarkovChain) = TemporalHistory{Int}

## Access

nb_states(mc::ContinuousMarkovChain) = length(mc.π0)

initial_distribution(mc::ContinuousMarkovChain) = mc.π0

rate_matrix(mc::ContinuousMarkovChain) = mc.Q

rate_negdiag(mc::ContinuousMarkovChain) = -diag(rate_matrix(mc))

function transition_matrix(mc::ContinuousMarkovChain)
    P = rate_matrix(mc) ./ rate_negdiag(mc)
    P[diagind(P)] .= zero(eltype(P))
    return P
end

## Simulation

function Base.rand(rng::AbstractRNG, mc::ContinuousMarkovChain, tmin, tmax)
    D = rate_negdiag(mc)
    P = transition_matrix(mc)
    transitions = [Dists.Categorical(P[s, :]) for s = 1:nb_states(mc)]  # TODO: @view
    waiting_times = [Dists.Exponential(1 / D[s]) for s = 1:nb_states(mc)]

    h = TemporalHistory(Float64[], Int[], tmin, tmax)
    s = rand(rng, Dists.Categorical(mc.π0))
    push!(h, tmin, s)
    t = tmin
    while t < tmax
        Δt = rand(rng, waiting_times[s])
        if t + Δt < tmax
            t += Δt
            s = rand(rng, transitions[s])
            push!(h, t, s)
        else
            break
        end
    end
    return h
end

Base.rand(mc::ContinuousMarkovChain, tmin, tmax) = rand(GLOBAL_RNG, mc, tmin, tmax)

## Prior

"""
    ContinuousMarkovChainPrior

Define a Dirichlet prior on the initial distribution and a Gamma prior on the transition rates from each state.

# Fields
- `π0α`: Dirichlet parameter
- `Pα`: Gamma rate parameters
- `Pβ`: Gamma shape parameters
"""
struct ContinuousMarkovChainPrior <: AbstractMarkovChainPrior
    π0α::Vector{Float64}
    Pα::Matrix{Float64}
    Pβ::Vector{Float64}
end

## Prior logpdf

function MeasureTheory.logdensity(
    prior::ContinuousMarkovChainPrior,
    mc::ContinuousMarkovChain,
)
    l = logdensity(Dists.Dirichlet(prior.π0α), Dists.Categorical(mc.π0))
    Q = rate_matrix(mc)
    for i = 1:nb_states(mc), j = 1:nb_states(mc)
        if i != j
            l += logdensity(Dists.Gamma(prior.Pα[i, i], 1 / prior.Pβ[i, j]), Q[i, j])
        end
    end
    return l
end

## Fitting

"""
    ContinuousMarkovChainStats

Sufficient statistics for the likelihood of a ContinuousMarkovChain.
"""
struct ContinuousMarkovChainStats <: Dists.SufficientStats
    initialization::Vector{Float64}
    transition_count::Matrix{Float64}
    duration::Vector{Float64}
end

function Dists.suffstats(::Type{ContinuousMarkovChain}, h::TemporalHistory{<:Integer})
    states = h.marks
    S = maximum(states)
    initialization = collect(1:S) .== states[1]
    duration = zeros(Float64, S)
    transition_count = zeros(Int, S, S)
    for t = 1:length(states)-1
        duration[states[t]] += h.times[t+1] - h.times[t]
        transition_count[states[t], states[t+1]] += 1
    end
    duration[states[length(states)]] += h.tmax - h.times[length(states)]
    return ContinuousMarkovChainStats(initialization, transition_count, duration)
end

function Dists.suffstats(::Type{ContinuousMarkovChain}, m̂, D̂)
    error("not implemented")
end

function Dists.suffstats(
    ::Type{<:ContinuousMarkovChain},
    prior::ContinuousMarkovChainPrior,
    args...,
)
    ss = suffstats(ContinuousMarkovChain, args...)
    ss.initialization .+= (prior.π0α .- 1)
    ss.transition_count .+= (prior.Pα .- 1)
    ss.duration .+= prior.Pβ
    return ss
end

function Dists.fit_mle(::Type{<:ContinuousMarkovChain}, ss::ContinuousMarkovChainStats)
    π0 = ss.initialization
    π0 ./= sum(π0)
    Q = ss.transition_count ./ ss.duration
    Q[diagind(Q)] .= -vec(sum(Q, dims = 2)) + Q[diagind(Q)]
    return ContinuousMarkovChain(π0, Q)
end

## Asymptotics

function stationary_distribution(mc::ContinuousMarkovChain)
    π = real.(eigvecs(rate_matrix(mc)')[:, end])
    return π / sum(π)
end

## Misc

function discretize_chain(mc::ContinuousMarkovChain, Δt)
    return DiscreteMarkovChain(initial_distribution(mc), exp(rate_matrix(mc) * Δt))
end


text:
using DSGE
using HDF5
path = dirname(@__FILE__)

# Test in model optimization
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

file = "$path/../reference/optimize.h5"
# For regenerating test file
#=params_test = h5read(file, "params")
data_test = h5read(file, "data") =#

x0 = h5read(file, "params")
data = h5read(file, "data")'
minimizer = h5read(file, "minimizer")
minimum = h5read(file, "minimum")
H_expected = h5read(file, "H")

# See src/estimate/estimate.jl
update!(m, x0)
n_iterations = 3

x0 = Float64[p.value for p in m.parameters]

out, H = optimize!(m, data; iterations=n_iterations)

# Re-generate test file
#=h5open("$path/../reference/optimize.h5", "w") do file
    file["params"] = params_test
    file["data"] = data_test
    file["minimizer"] = out.minimizer
    file["minimum"] = out.minimum
    file["H"] = H
end=#

@testset "Check optimize minimizers are the same" begin
    @test @test_matrix_approx_eq minimizer out.minimizer
    @test minimum ≈ out.minimum atol=5e-7
    @test @test_matrix_approx_eq H_expected H
end

nothing


text:
module AtomBase

include("Atom.jl")
include("Atoms.jl")

export AbstractAtom, Atom, AbstractProperties, Properties, Energy, Atoms
export addAtom!, addProperty!

end


text:
struct OrbitalF90
    npy  :: Cint
    naky :: Cint
    ey   :: Cdouble
    my   :: Cint
    pa   :: Ptr{Cdouble}
    qa   :: Ptr{Cdouble}
    ra   :: Ptr{Cdouble}
end

struct RWFNOrbital
    orbital_n :: Int
    orbital_kappa :: Int
    r :: Vector{Float64}
    p :: Vector{Float64}
    q :: Vector{Float64}
    # ....
    original :: OrbitalF90
    #energy :: Float64
end

function read_rwfn(filename)
    isfile(filename) || error("Unable to open $(filename)")

    orbitals = Ref{Ptr{OrbitalF90}}()
    norbitals = Ref{Cint}()
    status = ccall( (:rwfnread, libgrasp_so),
        Cint, (Cstring, Ref{Cint}, Ref{Ptr{OrbitalF90}}),
        filename, norbitals, orbitals
    )
    if status != zero(status)
        error("rwfnread returned with status $(status)")
    end
    orbitals = unsafe_wrap(Array{OrbitalF90}, orbitals.x, norbitals.x, own=true)
    map(orbitals) do orb
        RWFNOrbital(
            orb.npy, orb.naky,
            unsafe_wrap(Array{Cdouble}, orb.ra, (orb.my,), own=true),
            unsafe_wrap(Array{Cdouble}, orb.pa, (orb.my,), own=true),
            unsafe_wrap(Array{Cdouble}, orb.qa, (orb.my,), own=true),
            orb
        )
    end
end


text:
using ITensors

let
  N = 100
  cutoff = 1E-8
  tau = 0.1
  ttotal = 5.0

  # Make an array of 'site' indices
  s = siteinds("S=1/2", N; conserve_qns=true)

  # Make gates (1,2),(2,3),(3,4),...
  gates = ITensor[]
  for j in 1:(N - 1)
    s1 = s[j]
    s2 = s[j + 1]
    hj =
      op("Sz", s1) * op("Sz", s2) +
      1 / 2 * op("S+", s1) * op("S-", s2) +
      1 / 2 * op("S-", s1) * op("S+", s2)
    Gj = exp(-1.0im * tau / 2 * hj)
    push!(gates, Gj)
  end
  # Include gates in reverse order too
  # (N,N-1),(N-1,N-2),...
  append!(gates, reverse(gates))

  # Initialize psi to be a product state (alternating up and down)
  psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")

  c = div(N, 2) # center site

  # Compute and print <Sz> at each time step
  # then apply the gates to go to the next time
  for t in 0.0:tau:ttotal
    Sz = expect(psi, "Sz"; sites=c)
    println("$t $Sz")

    t ≈ ttotal && break

    psi = apply(gates, psi; cutoff)
    normalize!(psi)
  end

  return nothing
end


text:
module Dubins

using Memento

const LOGGER = getlogger(@__MODULE__)
Memento.setlevel!(LOGGER, "info")

__init__() = Memento.register(LOGGER)

"Suppresses information and warning messages output by Dubins, for fine grained control use the Memento package"
function silence()
    Memento.info(LOGGER, "Suppressing information and warning messages for the rest of this session.  Use the Memento package for more fine-grained control of logging.")
    Memento.setlevel!(Memento.getlogger(Dubins), "error")
end

"alows the user to set the logging level without the need to add Memento"
function logger_config!(level)
    Memento.config!(Memento.getlogger("Dubins"), level)
end

include("typedefs.jl")
include("paths.jl")
include("path_fcns.jl")

end


text:
# Autogenerated wrapper script for ITSOL_2_jll for i686-w64-mingw32-libgfortran3
export libITSOL_2

using LAPACK_jll
using OpenBLAS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ITSOL_2")
JLLWrappers.@declare_library_product(libITSOL_2, "libITSOL_2.dll")
function __init__()
    JLLWrappers.@generate_init_header(LAPACK_jll, OpenBLAS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libITSOL_2,
        "bin\\libITSOL_2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Forward

export Dual, ind_var

struct Dual{T<:AbstractFloat} <: AbstractFloat
    value::T
    grad::T
end

#===== Operations =====#
# https://en.wikipedia.org/wiki/Automatic_differentiation#Automatic_differentiation_using_dual_numbers

import Base: +, -, *, /, ^
u::Dual + v::Dual = Dual(u.value + v.value, u.grad + v.grad)
u::Dual - v::Dual = Dual(u.value - v.value, u.grad - v.grad)
u::Dual * v::Dual = Dual(u.value * v.value, v.value * u.grad + u.value * v.grad)
u::Dual / v::Dual = Dual(u.value * v.value, (v.value * u.grad - u.value * v.grad) / (v.value ^ 2))

u::Dual ^ k::Int = Dual(u.value ^ k, k * (u.value ^ (k - 1)) * u.grad)

Base.sin(u::Dual) = Dual(sin(u.value), u.grad * cos(u.value))
Base.cos(u::Dual) = Dual(cos(u.value), -u.grad * sin(u.value))

Base.exp(u::Dual) = begin
    exp_u = exp(u.value)
    Dual(exp_u, u.grad * exp_u)
end
Base.log(u::Dual) = Dual(log(u.value), u.grad / u.value)

#===== Conversions =====#

function Base.convert(::Type{Dual{T}}, x::Dual) where T
    Dual(convert(T, x.value), convert(T, x.grad))
end

function Base.convert(::Type{Dual{T}}, x::AbstractFloat) where T
    Dual(convert(T, x), zero(T))
end

function Base.promote_rule(::Type{Dual{T}}, ::Type{U}) where {T, U}
    Dual{promote_type(T, U)}
end

"Independent variable"
ind_var(x::AbstractFloat) = Dual(x, one(x))

end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_capacity(path::AbstractString, inputs::Dict, setup::Dict, EP::Model))

Function for writing the diferent capacities for the different generation technologies (starting capacities or, existing capacities, retired capacities, and new-built capacities).
"""
function write_capacity(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	# Capacity decisions
	dfGen = inputs["dfGen"]
	capdischarge = zeros(size(inputs["RESOURCES"]))
	for i in inputs["NEW_CAP"]
		if i in inputs["COMMIT"]
			capdischarge[i] = value(EP[:vCAP][i])*dfGen[!,:Cap_Size][i]
		else
			capdischarge[i] = value(EP[:vCAP][i])
		end
	end

	retcapdischarge = zeros(size(inputs["RESOURCES"]))
	for i in inputs["RET_CAP"]
		if i in inputs["COMMIT"]
			retcapdischarge[i] = first(value.(EP[:vRETCAP][i]))*dfGen[!,:Cap_Size][i]
		else
			retcapdischarge[i] = first(value.(EP[:vRETCAP][i]))
		end
	end

	capcharge = zeros(size(inputs["RESOURCES"]))
	retcapcharge = zeros(size(inputs["RESOURCES"]))
	for i in inputs["STOR_ASYMMETRIC"]
		if i in inputs["NEW_CAP_CHARGE"]
			capcharge[i] = value(EP[:vCAPCHARGE][i])
		end
		if i in inputs["RET_CAP_CHARGE"]
			retcapcharge[i] = value(EP[:vRETCAPCHARGE][i])
		end
	end

	capenergy = zeros(size(inputs["RESOURCES"]))
	retcapenergy = zeros(size(inputs["RESOURCES"]))
	for i in inputs["STOR_ALL"]
		if i in inputs["NEW_CAP_ENERGY"]
			capenergy[i] = value(EP[:vCAPENERGY][i])
		end
		if i in inputs["RET_CAP_ENERGY"]
			retcapenergy[i] = value(EP[:vRETCAPENERGY][i])
		end
	end
	dfCap = DataFrame(
		Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone],
		StartCap = dfGen[!,:Existing_Cap_MW],
		RetCap = retcapdischarge[:],
		NewCap = capdischarge[:],
		EndCap = value.(EP[:eTotalCap]),
		StartEnergyCap = dfGen[!,:Existing_Cap_MWh],
		RetEnergyCap = retcapenergy[:],
		NewEnergyCap = capenergy[:],
		EndEnergyCap = dfGen[!,:Existing_Cap_MWh]+capenergy[:]-retcapenergy[:],
		StartChargeCap = dfGen[!,:Existing_Charge_Cap_MW],
		RetChargeCap = retcapcharge[:],
		NewChargeCap = capcharge[:],
		EndChargeCap = dfGen[!,:Existing_Charge_Cap_MW]+capcharge[:]-retcapcharge[:]
	)
	if setup["ParameterScale"] ==1
		dfCap.StartCap = dfCap.StartCap * ModelScalingFactor
		dfCap.RetCap = dfCap.RetCap * ModelScalingFactor
		dfCap.NewCap = dfCap.NewCap * ModelScalingFactor
		dfCap.EndCap = dfCap.EndCap * ModelScalingFactor
		dfCap.StartEnergyCap = dfCap.StartEnergyCap * ModelScalingFactor
		dfCap.RetEnergyCap = dfCap.RetEnergyCap * ModelScalingFactor
		dfCap.NewEnergyCap = dfCap.NewEnergyCap * ModelScalingFactor
		dfCap.EndEnergyCap = dfCap.EndEnergyCap * ModelScalingFactor
		dfCap.StartChargeCap = dfCap.StartChargeCap * ModelScalingFactor
		dfCap.RetChargeCap = dfCap.RetChargeCap * ModelScalingFactor
		dfCap.NewChargeCap = dfCap.NewChargeCap * ModelScalingFactor
		dfCap.EndChargeCap = dfCap.EndChargeCap * ModelScalingFactor
	end
	total = DataFrame(
			Resource = "Total", Zone = "n/a",
			StartCap = sum(dfCap[!,:StartCap]), RetCap = sum(dfCap[!,:RetCap]),
			NewCap = sum(dfCap[!,:NewCap]), EndCap = sum(dfCap[!,:EndCap]),
			StartEnergyCap = sum(dfCap[!,:StartEnergyCap]), RetEnergyCap = sum(dfCap[!,:RetEnergyCap]),
			NewEnergyCap = sum(dfCap[!,:NewEnergyCap]), EndEnergyCap = sum(dfCap[!,:EndEnergyCap]),
			StartChargeCap = sum(dfCap[!,:StartChargeCap]), RetChargeCap = sum(dfCap[!,:RetChargeCap]),
			NewChargeCap = sum(dfCap[!,:NewChargeCap]), EndChargeCap = sum(dfCap[!,:EndChargeCap])
		)

	dfCap = vcat(dfCap, total)
	CSV.write(joinpath(path, "capacity.csv"), dfCap)
	return dfCap
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using Turing
using Turing: Sampler

alg = NUTS(1000, 200, 0.65)
println(alg)
sampler = Sampler(alg)

alg = NUTS(1000, 0.65)
println(alg)
sampler = Sampler(alg)

alg = NUTS(1000, 200, 0.65, :m)
println(alg)
sampler = Sampler(alg)


text:
include("getdigit_n.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n
using .MyPrintf

const DIGIT = 100
const MAXN = 100

function factorial(m)
    return m < 2 ? BigFloat(1) : reduce(*, [BigFloat(i) for i in 2:m])
end

function ramanujan(n)
    res = BigFloat(0)
    for k = 0:n
        numerator = factorial(4 * k) * (BigFloat(1103) + BigFloat(26390 * k))
        denominator = (factorial(k) * (BigFloat(396))^(BigFloat(k)))^(BigFloat(4))

        res += numerator / denominator
    end

    return BigFloat(9801) / (res * BigFloat(2) * sqrt(BigFloat(2)))
end

function main()
    n = GetDigit_n.getn(ramanujan, MAXN, DIGIT, true)
    precision = GetDigit_n.getdigit(ramanujan, n, DIGIT, true)
    precision = GetDigit_n.getdigit_2(ramanujan, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, ramanujan(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


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

using Markdown
using InteractiveUtils

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

# ╔═╡ ac87ced2-986f-4e6c-80b2-104b25c171c2
begin
	using CSV, JLD2, FITSIO, FileIO# , HDF5
	using DataFrames #, Query
	using PyCall
	using PlutoUI, PlutoTeachingTools
	eval(Meta.parse(code_for_check_type_funcs))
end

# ╔═╡ e0e60e09-3808-4d6b-a773-6ba59c02f517
md"""
# Astro 528, Lab 3, Exercise 2
# Benchmarking File I/O
# (no Python dependance)
"""

# ╔═╡ 8c1d81ab-d215-4972-afeb-7e00bf6063c2
md"""
For many applications, its important that we be able to read input data from a file and/or to write our outputes to files so they can be reused later.  Disk access is typically much slower than accessing system memory.  Therefore, disk access can easily become the limiting factor for a project.  In this set of exercises, you'll see examples of how to perform basic file I/O.

You'll be provided with most of the code you need, so that you can focus on comparing how much disk space and time is required by different file formats.  Near the end of the lab, you'll be asked to think about when each type of file format would be a good choice for you to use in your research projects.
"""

# ╔═╡ 1607eac9-e76f-4d1f-a9ce-981ce3be9bea
md"""
### Download some data
First, we're going to download some data from the web.  Julia has a built in `download` function that can be handy for this.  It relies on your system having some utilities already installed (e.g., `curl`, `wget` or `fetch`).  If you run this on a local system and run into trouble, then you can leave the cell below, and manually download the file to the data subdirectory.
"""

# ╔═╡ f27e1e8f-15eb-4754-a94c-7f37c54b871e
begin
	path = basename(pwd())=="test" ? "../data/" : "data/"
	if !isdir(path)  mkdir(path)  end  # make sure there's a data directory
	url = "https://personal.psu.edu/~ebf11/data/kplr_dr25_inj1_plti.csv"
	filename_csv = joinpath(path,basename(url)) # extract the filename and prepend "data/"
	time_to_download = NaN
	if !isfile(filename_csv)  # skip downloading if file already exists
	    time_to_download = @elapsed download(url,filename_csv)
	end
end

# ╔═╡ 80f02c3a-6751-48df-92ec-13f5c7d8c71e
if !isnan(time_to_download)
	md"Downloading the file took $time_to_download seconds."
end

# ╔═╡ 624c9038-3008-4e78-a149-60796dacf9c0
md"""
Previously, everything you needed for an assignment was included in a GitHub repository.  So why did I make you download the file?

Notice the size of the file.  Git is great for tracking source code, but it wasn't really designed for working with large files (especially large _binary_ files).  Since  we're not going to be editing it, we'll simply download it once.
"""

# ╔═╡ c808d0e1-4829-432f-b91f-968d4be20ecf
md"""
(We're skipping 2a-2d, since we're avoiding using Python in this notebook.)
"""

# ╔═╡ 191ba96e-2573-4bc1-a352-46a66e0a5c4f
md"""
## CSV files
Let's say that we'd like to read/write data as [CSV files](https://en.wikipedia.org/wiki/Comma-separated_values), so that it's easy for other programs to read in.  Using the CSV package, we can read a CSV file into a DataFrame (or several other tabular data structures) and write CSV files from a DataFrame with code like the following.
"""

# ╔═╡ 82a757ad-566d-4c1d-8b3d-366ffd980fb4
begin
	# Read & write a small test file so compilation time not included below
	n_small_df = 1024
	small_csv_filename = "random_numbers.csv"
	small_df = DataFrame(a=rand(n_small_df),b=rand(n_small_df) )
	CSV.write(joinpath(path,small_csv_filename),small_df)
	small_df_from_csv = CSV.read(joinpath(path,small_csv_filename),DataFrame,types=Dict(:TCE_ID=>String))
end;

# ╔═╡ 32124f8a-f5cf-41d1-97a9-ce1d05145fde
md"""
2d. How long do you think it will take to read 146,294 rows and 25 columns of data data stored in CSV format?
"""

# ╔═╡ 507e1516-5433-49eb-831d-32426f30895e
response_2d = missing

# ╔═╡ eac67cc9-754b-4f7d-add8-93900a1b5b49
display_msg_if_fail(check_type_isa(:response_2d,response_2d,[AbstractString,Markdown.MD]))

# ╔═╡ 73c06bb1-be49-46bd-b7f1-c45cc56af7b4
md"""
I'm ready to benchmark reading a CSV file. $(@bind ready_read_csv CheckBox())
$(@bind go_read_csv Button("Rerun the benchmarks."))
"""

# ╔═╡ ebfaa677-829b-4bf9-bdbb-19f3c87dd3a4
if  ready_read_csv
	go_read_csv
	time_read_csv = @elapsed df_csv = CSV.read(filename_csv, DataFrame)
	already_benchmarked_csv_read = true
	md"It took $time_read_csv seconds to read the CSV file."
end

# ╔═╡ b0c3875b-ef07-4e92-a0a8-55f42b266c6b
if already_benchmarked_csv_read
	ready_read_csv
	df = CSV.read(filename_csv,DataFrame,types=Dict(:TCE_ID=>String),missingstring="NA")
end

# ╔═╡ 945f5a55-3026-4497-9ece-8af878c87788
md"""
2e.  How did the time required to read the data in CSV format compare to your expectation?  If you were suprsised, try to explain what happened.
"""

# ╔═╡ 57397ee4-9efc-48b3-b640-d2b7a10da633
response_2e = missing

# ╔═╡ 8059a6a3-384a-4344-8a23-650ee0be10c2
display_msg_if_fail(check_type_isa(:response_2e,response_2e,[AbstractString,Markdown.MD]))

# ╔═╡ 64224c6b-c5a0-44f2-b2a0-7f77759cb848
md"""
Now, we'll try writing the same data out to a new CSV file.
"""

# ╔═╡ f5b93929-2c59-4360-8c41-97a1324ba455
md"2f. How long do you think it will take to write the same data?  Once you've made your prediction, mouse over the hint box.  If you were suprsised, try to explain what happened."

# ╔═╡ 122196fa-45ca-4031-85eb-4afd4782de9e
response_2f = missing

# ╔═╡ e9dc1456-616b-4e4b-b209-9f6ba4c48607
display_msg_if_fail(check_type_isa(:response_2f,response_2f,[AbstractString,Markdown.MD]))

# ╔═╡ 3e550b71-4750-460b-be18-911a848a8f49
if  ready_read_csv
	go_read_csv
	filename_csv_out = replace(filename_csv, ".csv" => "_2.csv")
	time_write_csv = @elapsed CSV.write(filename_csv_out, df)
end;

# ╔═╡ f15d37a7-d962-4da0-977f-76729a3313be
hint(md"It took $time_write_csv seconds to write the CSV file.")

# ╔═╡ 28f195c4-4f61-4873-85d6-b4e3aaa3660f
md"""
## Binary formats: HDF5/JLD2

There are numerous binary file formats that one could use.  Here, we'll try using JLD2 which is a subset of the [HDF5](https://www.hdfgroup.org/solutions/hdf5/) file format.  This means that when [Julia's JLD2 package](https://github.com/JuliaIO/JLD2.jl) writes jld2 files, they can be read by other programs that can read HDF5 files.  However, a generic HDF5 file is not a valid JLD2 file.  If you want to read a HDF5 file, then you can use Julia's [HDF5.jl package](https://github.com/JuliaIO/HDF5.jl).  The [FileIO.jl](https://github.com/JuliaIO/FileIO.jl) package provides a common interface for reading and writing from multiple file formats, including these and several others.

As before, we'll call each function once using a small DataFrame, just so they get compiled before we benchmark them.
"""

# ╔═╡ 3837e439-250b-4577-b0d7-93352ec19f6e
begin
	filename_jld2_small = replace(small_csv_filename, ".csv" => ".jld2")
	@save joinpath(path,filename_jld2_small) small_df
	small_df_from_jld2 = load(joinpath(path,filename_jld2_small), "small_df")
end;

# ╔═╡ df9b701f-d314-4512-b2ea-1f6ae015166c
response_2g = missing

# ╔═╡ 377f7527-a338-4526-bb24-9766c635e719
display_msg_if_fail(check_type_isa(:response_2g,response_2g,[AbstractString,Markdown.MD]))

# ╔═╡ 25f24754-16d4-4343-b54a-cf8ea1358ce9
begin
	small_jld2_filesize = filesize(joinpath(path,filename_jld2_small))/1024
	small_csv_filesize = filesize(joinpath(path,small_csv_filename))/1024
	hint(md"For the small random data set, the JLD2 file's size is $small_jld2_filesize KB versus $small_csv_filesize KB for the CSV file.")
end

# ╔═╡ 8255b5a9-cbe6-4604-bedd-0e366f311096
md"""
Let's check the filesizes for the JLD2 file to the CSV file.
The small CSV file we created is $small_csv_filesize KB.

2g.  How large would you guess the JLD2 file will be?  Once you've made your prediction, mouse over the hint box.  If you were suprised, try to explain what happened.
"""

# ╔═╡ 570fd826-23fd-46ee-bdb4-58fb0c45719a
md"""
Now let's time how long it takes to save the data from CSV into a JLD2 file.
"""

# ╔═╡ 691410bb-0472-4800-a90d-29ddf947de3e
begin
	filename_jld2 = replace(filename_csv, ".csv" => ".jld2")
	with_terminal() do
		@time @save filename_jld2 df
	end
end

# ╔═╡ 8cbb1c90-bd94-44b5-80b6-81d38f3e6252
md"2h.  How long do you think it will take to load the data from the JLD2 file? "

# ╔═╡ c3065acf-6205-455f-ba74-ca51f3f6761b
response_2h = missing

# ╔═╡ 9c392be9-1505-40bc-a290-68085a1c2700
display_msg_if_fail(check_type_isa(:response_2h,response_2h,[AbstractString,Markdown.MD]))

# ╔═╡ d738bdcc-5f83-4dfa-a17f-e9ea23db2986
md"""
Now time how long it takes to load the data from the JLD2 file.
"""

# ╔═╡ d55ce157-099c-4c1b-94db-62918f04e5fe
md"""
I'm ready to benchmark reading a JLD2 file. $(@bind ready_read_jld2 CheckBox())
$(@bind go_read_jld2 Button("Rerun the benchmarks."))
"""

# ╔═╡ 206f464b-55fe-46aa-85b7-8f0246a0aaad
if ready_read_jld2
	go_read_jld2
	time_read_jld2 = @elapsed df_jld2_read_back_in = load(filename_jld2, "df")
	md"It took $time_read_jld2 seconds to read the JLD2 file."
end

# ╔═╡ 6f72d1b2-63f6-4301-8272-bb2d6d2d049e
md"""
2i.  How does the time required to read and write the JLD2 file compare to the time required to read the CSV formatted files?  If you were suprised, try to explain what happend.
"""

# ╔═╡ 01201c37-0b79-46b1-a001-e716f5b3ba67
response_2i = missing

# ╔═╡ e3d37bc1-a119-4add-a111-899ee0caea05
display_msg_if_fail(check_type_isa(:response_2i,response_2i,[AbstractString,Markdown.MD]))

# ╔═╡ 39eff7cf-6f6c-4aac-91a7-b776e0a62c45
csv_filesize = filesize(filename_csv) / 1024^2;

# ╔═╡ 6dff3f21-fac0-42e1-910a-f969a231374f
md"""
Next, we'll compare the filesizes for the JLD2 file to the CSV file. The CSV file is $csv_filesize MB.

2j.  How large would you guess the JLD2 file will be?  Once you've made your prediction, mouse over the hint box.  If you were suprised, try to explain what happened.
"""

# ╔═╡ b26b8253-e6cd-49f4-81c5-2a3c2963a37c
response_2j = missing

# ╔═╡ e5dc123f-1596-4311-9398-f0cfe80a5342
display_msg_if_fail(check_type_isa(:response_2j,response_2j,[AbstractString,Markdown.MD]))

# ╔═╡ 7b9d26c2-899e-45e9-b664-39d5f1adfe3f
begin
	jld2_filesize = filesize(filename_jld2) / 1024^2;
	hint(md"The JLD2 file size is $jld2_filesize MB.")
end

# ╔═╡ fc01d57f-c90b-4231-96be-ddd48656d55e
md"""
## Flexible Format: FITS

Astronomers often use the [FITS file format](https://en.wikipedia.org/wiki/FITS).  Like [HDF
5](https://www.hdfgroup.org/solutions/hdf5/), it's a very flexible (e.g., it can store both text and binary data) and thus complicated file
 format.
Therefore, most languages call a common [FITSIO library written in C](https://heasarc.gsfc.nasa.gov/fitsio/), rather than implementing code themselves.  Indeed, that's what [Julia's FITSIO.jl package](https://github.com/JuliaAstro/FITSIO.jl) does.

Unfortunately, the FITSIO package isn't as polished as the others.  It expects a `Dict` rather than a `DataFrame`, and it can't handle missing values.  So I've provided some helper functions at the bottom of the notebook.  Also, FITS files have complicated headers, so I'll provide a function to read all the tabular data from a simple FITS file.  As usual, we'll use each function once, so that Julia compiles them before we start timing.
"""

# ╔═╡ 3b232365-f2fe-4edb-a39f-3e37c8cbb666
md"Now we can time how long it takes to write and read the data as FITS files."

# ╔═╡ 568862b3-6fce-426a-a9e2-e558adf3932a
md"""
I'm ready to benchmark reading a FITS file. $(@bind ready_read_fits CheckBox())
$(@bind go_read_fits Button("Rerun the benchmarks."))
"""

# ╔═╡ ae79ee89-d788-4612-8b1d-fc22d85c7744
md"2k.  How do the read/write times and file sizes for FITS compare to CSV and JLD2?"

# ╔═╡ 3745c237-ba48-4f2c-959e-441484244764
response_2k = missing

# ╔═╡ 6d061411-6f19-4119-aefb-cc380198b0ce
display_msg_if_fail(check_type_isa(:response_2k,response_2k,[AbstractString,Markdown.MD]))

# ╔═╡ 8def87d2-f10b-4a82-b353-a6477eeead9b
md"## Implications for your project"

# ╔═╡ c74b3105-f480-4688-b85f-3e7dff70da3b
md"""
2l.  How does the time required to read any of the above formats compare to the time required to download the files ($time_to_download seconds)?
Will your project need to transfer large files over the internet?  If so, very roughly how large do you expect they will be?
"""

# ╔═╡ 55438c09-1d94-4ff7-90c3-0cc6064a091e
response_2l = missing

# ╔═╡ bcc796c9-db11-4a09-a5f9-215127ac0938
display_msg_if_fail(check_type_isa(:response_2l,response_2l,[AbstractString,Markdown.MD]))

# ╔═╡ 92b2ecdd-9491-4f16-8cb9-bacbcf180280
md"""
2m.  Will your project need to read in any input files?  If so, what format are they in?

I assume everyone's project code will write at least some output files?  Very roughly, how large do you expect that they'll be?   What format(s) would be a good choice for your project?"""

# ╔═╡ 83216915-cdcf-4f0f-829f-a5ff4c4b8da0
response_2m = missing

# ╔═╡ 452872f9-2009-4733-b2d3-28f262ae19b7
display_msg_if_fail(check_type_isa(:response_2m,response_2m,[AbstractString,Markdown.MD]))

# ╔═╡ 29415ddc-e002-4f56-a169-95f7b1c36be9
md"# Helper Functions"

# ╔═╡ fb23d6c6-b812-4fe1-b224-0014bedbd43f
ChooseDisplayMode()

# ╔═╡ 1e53aa10-dff6-40d5-89e2-da194ffc2052
TableOfContents()

# ╔═╡ 14cca8ce-cc61-4fae-b871-21c3fd23d0ea
"Convert a DataFrame to a Dict"
function convert_dataframe_to_dict_remove_missing(df::DataFrame)
#=
	"Convert a DataFrame to a Dict, replacing missing values with 0 or an empty string."
	d = Dict(map(n->"$n"=>               # create a dictionary
              ( any(ismissing.(df[!,n])) ? # if column contains a missing
                    map(x-> !ismissing(x) ? # search for missings
                        x :                 # leave non-missing values alone
                        ( (eltype(df[!,n]) <: Number) ? zero(eltype(n)) : "")
                        , df[!,n])            # but replace missing with 0 or ""
                : df[!,n] ), # if nothing is missing, just use column as is
            names(df) ))
=#
	d = Dict(zip(names(df),collect.(eachcol(df))))
end

# ╔═╡ 28fc8de4-749b-4093-b32f-c398f8d27d3d
"Write a DataFrame to a FITS file, replacing missing values with 0 or an empty string."
function write_dataframe_as_fits(filename::String, df::DataFrame)
    try
       dict = convert_dataframe_to_dict_remove_missing(df)
       fits_file = FITS(filename,"w")
       write(fits_file, dict )
       close(fits_file)
    catch
        @warn("There was a problem writing a dataframe to " * filename * ".")
    end
end

# ╔═╡ 57b422e5-0ad0-4674-bdd3-a8358bc7aaeb

"Read the columns of the first table from a FITS file into a Dict"
function read_fits_tables(filename::String)
    dict = Dict{String,Any}()
    fits_file = FITS(filename,"r")
    # fits_file[1] is image data, we're interested in the table
    @assert length(fits_file) >= 2
    header = read_header(fits_file[2])
    for i in 1:length(header)
        c = get_comment(header,i)
        if !occursin("label for field",c)
            continue
        end
        h = header[i]
        @assert typeof(h) == String
        try
            dict[h] = read(fits_file[2],h)
        catch
            @warn "# Problem reading table column " * h * "."
        end
    end
    close(fits_file)
    return dict
end

# ╔═╡ e05f16d6-eb50-49a9-bf14-95d63c9da7ff
begin
	filename_fits_small = replace(small_csv_filename, ".csv" => ".fits")
	write_dataframe_as_fits(joinpath(path,filename_fits_small),small_df)
	read_fits_tables(joinpath(path,filename_fits_small))
end;

# ╔═╡ 1992595f-9976-4b18-bf9c-df8a73d30dc8
begin
	filename_fits_small # make sure have already compiled functions
	filename_fits = replace(filename_csv, ".csv" => ".fits")
	with_terminal() do
		@time write_dataframe_as_fits(filename_fits,df)
	end
end

# ╔═╡ 90d5244e-17be-4601-b922-8c254f1248bf
begin
	fits_filesize = filesize(filename_fits) /1024^2
	hint(md"The FITS file size is $fits_filesize MB versus $jld2_filesize MB for the JLD2 and $csv_filesize MB for the CSV.")
end

# ╔═╡ 52f9edd0-79b0-4a9a-9930-3a05d3aa2447
if ready_read_fits
	go_read_fits
	time_read_fits = @elapsed read_fits_tables(filename_fits)
	md"It took $time_read_fits seconds to read the FITS file."
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
PlutoTeachingTools = "661c6b06-c737-4d37-b85c-46df65de6f69"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.2.2"
FITSIO = "~0.16.7"
FileIO = "~1.10.1"
JLD2 = "~0.4.13"
PlutoTeachingTools = "~0.1.3"
PlutoUI = "~0.7.9"
PyCall = "~1.92.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CFITSIO]]
deps = ["CFITSIO_jll"]
git-tree-sha1 = "c860f5545064216f86aa3365ec186ce7ced6a935"
uuid = "3b1b4be9-1499-4b22-8d78-7db3344d1961"
version = "1.3.0"

[[CFITSIO_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "Pkg"]
git-tree-sha1 = "2fabb5fc48d185d104ca7ed7444b475705993447"
uuid = "b3e40c51-02ae-5482-8a39-3ace5868dcf4"
version = "3.49.1+0"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FITSIO]]
deps = ["CFITSIO", "Printf", "Reexport", "Tables"]
git-tree-sha1 = "85b66005c9d16d3d27c9d06fd3be5e25074128da"
uuid = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
version = "0.16.7"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.10.1"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLD2]]
deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "59ee430ac5dc87bc3eec833cc2a37853425750b4"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.13"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "2ca267b08821e86c5ef4376cffed98a46c2cb205"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoTeachingTools]]
deps = ["LaTeXStrings", "Markdown", "PlutoUI", "Random"]
git-tree-sha1 = "e2b63ee022e0b20f43fcd15cda3a9047f449e3b4"
uuid = "661c6b06-c737-4d37-b85c-46df65de6f69"
version = "0.1.4"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.3"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "54f37736d8934a12a200edea2f9206b03bdf3159"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.7"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.0"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─e0e60e09-3808-4d6b-a773-6ba59c02f517
# ╟─8c1d81ab-d215-4972-afeb-7e00bf6063c2
# ╟─1607eac9-e76f-4d1f-a9ce-981ce3be9bea
# ╠═f27e1e8f-15eb-4754-a94c-7f37c54b871e
# ╠═80f02c3a-6751-48df-92ec-13f5c7d8c71e
# ╟─624c9038-3008-4e78-a149-60796dacf9c0
# ╟─c808d0e1-4829-432f-b91f-968d4be20ecf
# ╟─191ba96e-2573-4bc1-a352-46a66e0a5c4f
# ╠═82a757ad-566d-4c1d-8b3d-366ffd980fb4
# ╟─32124f8a-f5cf-41d1-97a9-ce1d05145fde
# ╠═507e1516-5433-49eb-831d-32426f30895e
# ╟─eac67cc9-754b-4f7d-add8-93900a1b5b49
# ╟─73c06bb1-be49-46bd-b7f1-c45cc56af7b4
# ╟─ebfaa677-829b-4bf9-bdbb-19f3c87dd3a4
# ╠═b0c3875b-ef07-4e92-a0a8-55f42b266c6b
# ╟─945f5a55-3026-4497-9ece-8af878c87788
# ╠═57397ee4-9efc-48b3-b640-d2b7a10da633
# ╟─8059a6a3-384a-4344-8a23-650ee0be10c2
# ╟─64224c6b-c5a0-44f2-b2a0-7f77759cb848
# ╟─f5b93929-2c59-4360-8c41-97a1324ba455
# ╠═122196fa-45ca-4031-85eb-4afd4782de9e
# ╟─e9dc1456-616b-4e4b-b209-9f6ba4c48607
# ╟─3e550b71-4750-460b-be18-911a848a8f49
# ╟─f15d37a7-d962-4da0-977f-76729a3313be
# ╟─28f195c4-4f61-4873-85d6-b4e3aaa3660f
# ╠═3837e439-250b-4577-b0d7-93352ec19f6e
# ╠═8255b5a9-cbe6-4604-bedd-0e366f311096
# ╠═df9b701f-d314-4512-b2ea-1f6ae015166c
# ╟─377f7527-a338-4526-bb24-9766c635e719
# ╠═25f24754-16d4-4343-b54a-cf8ea1358ce9
# ╟─570fd826-23fd-46ee-bdb4-58fb0c45719a
# ╠═691410bb-0472-4800-a90d-29ddf947de3e
# ╟─8cbb1c90-bd94-44b5-80b6-81d38f3e6252
# ╠═c3065acf-6205-455f-ba74-ca51f3f6761b
# ╟─9c392be9-1505-40bc-a290-68085a1c2700
# ╟─d738bdcc-5f83-4dfa-a17f-e9ea23db2986
# ╟─d55ce157-099c-4c1b-94db-62918f04e5fe
# ╟─206f464b-55fe-46aa-85b7-8f0246a0aaad
# ╟─6f72d1b2-63f6-4301-8272-bb2d6d2d049e
# ╠═01201c37-0b79-46b1-a001-e716f5b3ba67
# ╟─e3d37bc1-a119-4add-a111-899ee0caea05
# ╟─39eff7cf-6f6c-4aac-91a7-b776e0a62c45
# ╟─6dff3f21-fac0-42e1-910a-f969a231374f
# ╠═b26b8253-e6cd-49f4-81c5-2a3c2963a37c
# ╟─e5dc123f-1596-4311-9398-f0cfe80a5342
# ╟─7b9d26c2-899e-45e9-b664-39d5f1adfe3f
# ╟─fc01d57f-c90b-4231-96be-ddd48656d55e
# ╠═e05f16d6-eb50-49a9-bf14-95d63c9da7ff
# ╟─3b232365-f2fe-4edb-a39f-3e37c8cbb666
# ╠═1992595f-9976-4b18-bf9c-df8a73d30dc8
# ╟─568862b3-6fce-426a-a9e2-e558adf3932a
# ╟─52f9edd0-79b0-4a9a-9930-3a05d3aa2447
# ╟─90d5244e-17be-4601-b922-8c254f1248bf
# ╟─ae79ee89-d788-4612-8b1d-fc22d85c7744
# ╠═3745c237-ba48-4f2c-959e-441484244764
# ╟─6d061411-6f19-4119-aefb-cc380198b0ce
# ╟─8def87d2-f10b-4a82-b353-a6477eeead9b
# ╟─c74b3105-f480-4688-b85f-3e7dff70da3b
# ╠═55438c09-1d94-4ff7-90c3-0cc6064a091e
# ╟─bcc796c9-db11-4a09-a5f9-215127ac0938
# ╟─92b2ecdd-9491-4f16-8cb9-bacbcf180280
# ╠═83216915-cdcf-4f0f-829f-a5ff4c4b8da0
# ╟─452872f9-2009-4733-b2d3-28f262ae19b7
# ╟─29415ddc-e002-4f56-a169-95f7b1c36be9
# ╟─fb23d6c6-b812-4fe1-b224-0014bedbd43f
# ╠═1e53aa10-dff6-40d5-89e2-da194ffc2052
# ╠═ac87ced2-986f-4e6c-80b2-104b25c171c2
# ╠═14cca8ce-cc61-4fae-b871-21c3fd23d0ea
# ╠═28fc8de4-749b-4093-b32f-c398f8d27d3d
# ╠═57b422e5-0ad0-4674-bdd3-a8358bc7aaeb
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
baremodule MemoryAccesses

function record end
function init end
function clear end
function isfull end

function reference end

function copy end
function view end

function table end

function location_time end

module Internal

using ..MemoryAccesses: MemoryAccesses

import Tables
using AutoHashEquals: @auto_hash_equals
using Dates: Dates, Period, Microsecond, Nanosecond
using Requires: @require

if !@isdefined(Returns)
    Returns(x) = (_args...; _kwargs...) -> x
end

include("record.jl")
include("reference.jl")
include("table.jl")
include("analysis.jl")
include("plotter.jl")

function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plots_show.jl")

    MemoryAccesses.init()
end

end  # module Internal

end  # baremodule MemoryAccesses


text:
"""
Epoch of a Spike Train `x`, where `min <= x[i] < max`.

kwargs `isminzero`, `ismaxzero` and `shift` set epoch zero value to `min+shift` or `max+shift`.

return:
- y: epoch of the Spike Train
- n: number of elements in the epoch, or divided by duration(max-min) of the epoch, based on kwargs `israte`
- w: epoch window(min,max)
- i: indices of epoch elements in the original Spike Train, such that y = x[i]

See also: [`epochspiketrains`](@ref), [`epochspiketrainresponse`](@ref)
"""
function epochspiketrain(x,min::Real,max::Real;isminzero::Bool=false,ismaxzero::Bool=false,shift::Real=0,israte::Bool=false)
    if ismaxzero && isminzero
        error("Zero setting conflicts, only one of them is allowed to be true.")
    end
    i = findall(min .<= x .< max)
    w = (min, max)
    n = length(i)
    y = x[i]
    if isminzero
        y .-= min+shift
    end
    if ismaxzero
        y .-= max+shift
    end
    if israte
        n /= ((max-min)*SecondPerUnit)
    end
    return (;y,n,w,i)
end
function epochspiketrain(x,mins,maxs;isminzero::Bool=false,ismaxzero::Bool=false,shift::Real=0,israte::Bool=false)
    yn = length(mins)
    if yn != length(maxs)
        error("Length of mins and maxs do not match.")
    end
    ys = Vector{Vector{Float64}}(undef,yn)
    ns = Vector{Float64}(undef,yn)
    ws = Vector{Tuple{Float64,Float64}}(undef,yn)
    is = Vector{Vector{Int}}(undef,yn)
    for i in 1:yn
        ys[i],ns[i],ws[i],is[i] = epochspiketrain(x,mins[i],maxs[i],isminzero=isminzero,ismaxzero=ismaxzero,shift=shift,israte=israte)
    end
    return (y=ys,n=ns,w=ws,i=is)
end
epochspiketrain(x,minmaxs::AbstractMatrix;isminzero::Bool=false,ismaxzero::Bool=false,shift::Real=0,israte::Bool=false) = @views epochspiketrain(x,minmaxs[:,1],minmaxs[:,2],isminzero=isminzero,ismaxzero=ismaxzero,shift=shift,israte=israte)
"Epochs of a Spike Train in between binedges"
function epochspiketrain(x,binedges;isminzero::Bool=false,ismaxzero::Bool=false,shift::Real=0,israte::Bool=false)
    nbinedges = length(binedges); nbinedges<2 && error("Have $nbinedges binedges, need at least two binedges.")
    epochspiketrain(x,binedges[1:end-1],binedges[2:end],isminzero=isminzero,ismaxzero=ismaxzero,shift=shift,israte=israte)
end

"Epochs for each Spike Trains"
function epochspiketrains(xs,binedges;isminzero::Bool=false,ismaxzero::Bool=false,shift::Real=0,israte::Bool=false)
    tn = length(xs)
    yss = Vector{Vector{Vector{Float64}}}(undef,tn)
    nss = Vector{Vector{Float64}}(undef,tn)
    ws = nothing
    iss = Vector{Vector{Vector{Int}}}(undef,tn)
    for i in 1:tn
        yss[i],nss[i],ws,iss[i] = epochspiketrain(xs[i],binedges,isminzero=isminzero,ismaxzero=ismaxzero,shift=shift,israte=israte)
    end
    return (y=yss,n=nss,w=ws,i=iss)
end

"""
Response of each epoch of a Spike Train, could be mean firing rate or number of spikes based on kwarg `israte`.

See also: [`epochspiketrainresponse_ono`](@ref)
"""
epochspiketrainresponse(x,minmaxs::AbstractMatrix;israte::Bool=true) = epochspiketrain(x,minmaxs,israte=israte)[2]
epochspiketrainresponse(x,mins,maxs;israte::Bool=true) = epochspiketrain(x,mins,maxs,israte=israte)[2]

"""
Response of each epoch of a Spike Train, could be mean firing rate or number of spikes based on kwarg `israte`.

!!! note
    This is a faster(~500x) version compared to [`epochspiketrainresponse`](@ref), but only works when `x`, `mins` and `maxs` are ascendingly ordered, and each `maxs-mins` range are none-overlapping.
"""
epochspiketrainresponse_ono(x,minmaxs::AbstractMatrix;israte::Bool=true,isnan2zero::Bool=true) = @views epochspiketrainresponse_ono(x,minmaxs[:,1],minmaxs[:,2];israte,isnan2zero)
function epochspiketrainresponse_ono(x,mins,maxs;israte::Bool=true,isnan2zero::Bool=true)
    n = length(mins)
    if n != length(maxs)
        error("Length of mins and maxs do not match.")
    end
    ns = zeros(n)
    ni=1
    for v in x
        @label start
        if mins[ni]<=v
            if v<maxs[ni]
                ns[ni] += 1
            else
                ni+=1
                ni>n && break
                @goto start
            end
        end
    end
    if israte
        ns = ns ./ ((maxs.-mins).*SecondPerUnit)
        isnan2zero && replace!(ns,NaN=>0)
    end
    return ns
end

function subrvr_onotest(rv::RealVector,mins::RealVector,maxs::RealVector;israte::Bool=true,isnan2zero::Bool=true,isbin::Bool=true)
    n = length(mins)
    if n != length(maxs)
        error("Length of mins and maxs do not match.")
    end
    ns = zeros(n)  # store FR in condition on/off bins
    ni=1  # bin counter
    sp=0
    for v in rv
        @label start
        if (mins[ni]<=v) && (v<maxs[ni])

        else
            ni+=1
            ni>n && break
            @goto start
        end
        ns[ni] += 1
    end
    if israte
        ns = ns ./ ((maxs.-mins).*SecondPerUnit)
        isnan2zero && replace!(ns,NaN=>0)
    end
    return ns
end

function subrmr_ono(rv::RealVector,mins::RealVector,maxs::RealVector;israte::Bool=true,isnan2zero::Bool=true)
end

"Generate a Poisson Spike Train"
function poissonspiketrain(r,dur;t0=0.0,rp=2.0)
    isid = Exponential(1000/r)
    st = Float64[-rp]
    while true
        i = rand(isid)
        if i > rp
            s = i + st[end]
            if s <= dur
                push!(st,s)
            else
                break
            end
        end
    end
    deleteat!(st,1)
    return st.+t0
end
function poissonspiketrain(ir::Function,dur;t0=0.0,rp=2.0)
    st=Float64[-rp]
    for t in 0.05:0.1:dur
        rand()<=ir(t)/10000 && (t-st[end])>rp && push!(st,t)
    end
    deleteat!(st,1)
    return st.+t0
end

"Flat Spike Trains to SpikeTimes and Trials, optionally sort trial based on `sv`"
function flatspiketrains(xs,sv=[])
    tn = length(xs)
    if isempty(sv)
        issort=false
    elseif length(sv)==tn
        issort=true
    else
        @warn """Length of "xs" and "sv" do not match, sorting ignored."""
        issort=false
    end
    if issort
        sxs=xs[sortperm(sv)]
        ssv=sort(sv)
    else
        sxs=xs
        ssv=sv
    end
    x=Float64[];y=Float64[];s=[]
    for i in 1:tn
        v = sxs[i];n=length(v)
        n==0 && continue
        append!(x,v);append!(y,fill(i,n))
        issort && append!(s,fill(ssv[i],n))
    end
    return x,y,s
end

"Vertical stack same length vectors to matrix"
function vstack(xs)
    tn = length(xs)
    n = length(xs[1])
    mat = Matrix{Float64}(undef,tn,n)
    for i in 1:tn
        mat[i,:] = xs[i]
    end
    return mat
end

"Vertical Mean and SEM of a matrix"
function vmeanse(mat::AbstractMatrix;normfun=nothing)
    n = size(mat,1)
    if !isnothing(normfun)
        for i=1:n
            mat[i,:]=normfun(mat[i,:])
        end
    end
    m = dropdims(mean(mat,dims=1),dims=1)
    se = dropdims(std(mat,dims=1),dims=1)/sqrt(n)
    return (;m,se)
end

"PSTH of Spike Trains"
function psthspiketrains(xs,binedges;israte::Bool=true,ismean::Bool=true,normfun=nothing)
    nss = epochspiketrains(xs,binedges,israte=israte)[2]
    halfbinwidth = (binedges[2]-binedges[1])/2
    x = binedges[1:end-1].+halfbinwidth
    mat = vstack(nss)
    if ismean
        m,se = vmeanse(mat,normfun=normfun)
        return (;m,se,x)
    else
        return (;mat,x)
    end
end