text:
function test_fltype()
  Tf = Float16
  Tc = Complex{Tf}
  Tf == fltype(Tf) && Tf == fltype(Tc) && Tf == fltype((Tc, Tf)) && Tf == fltype((Tf, Tc))
end

function test_cxtype()
  Tf = Float16
  Tc = Complex{Tf}
  Tc == cxtype(Tf) && Tc == cxtype(Tc)
end

function test_innereltype(T=Complex{Float32})
  a = [zeros(T, 3), zeros(T, 2, 5)]
  T == innereltype(a)
end

function test_superzeros()
  T = Float64
  dims1 = (2, 3)
  dims2a = (7,)
  dims2b = dims1
  dims2 = (dims2a, dims2b)

  a1 = rand(T, dims1)
  a2 = [(1+im)*rand(T, dims2a), (2+im)*rand(T, dims2b)]

  a1a = zeros(T, dims1)
  a2a = [zeros(Complex{T}, dims2a), zeros(Complex{T}, dims2b)]

  a1z = superzeros(a1)
  a2z = superzeros(a2)
  a1d = superzeros(dims1)
  a2d = superzeros(Complex{T}, dims2)

  @superzeros innereltype(a1) a1 a1ma dum
  @superzeros innereltype(a2) a2 a2ma dum
  @superzeros innereltype(a1) dims1 a1md dum
  @superzeros innereltype(a2) dims2 a2md dum

  ( a1a == a1z && a1a == a1d && a1a == a1ma && a1a == a1md &&
    a2a == a2z && a2a == a2d && a2a == a2ma && a2a == a2md )
end

function test_supertuplezeros(; T1=Float64, T2=Complex{Float64}, dims1=(1,), dims2=(3, 3))
  T = (T1, T2)
  dims = (dims1, dims2)
  a = superzeros(T, dims)
  ( eltype(a[1]) == T1 && eltype(a[2]) == T2 &&
    size(a[1]) == dims1 && size(a[2]) == dims2 )
end

# This test could use some further work.
function test_radialspectrum(dev::Device, n, ahkl, ahρ; debug=false, atol=0.1, rfft=false)
  g = TwoDGrid(dev, n, 2π)
  if rfft==true
    ah = @. ahkl(g.kr, g.l)
  else
    ah = @. ahkl(g.k, g.l)
  end
  ah[1, 1] = 0.0

  ρ, ahρ_estimate = FourierFlows.radialspectrum(ah, g; refinement=16)

  if debug
    println(sum(ahρ.(ρ) - ahρ_estimate))
    return ah, ahρ_estimate, ahρ.(ρ)
  else
    normalizeddiff = sum(abs.(ahρ.(ρ) - ahρ_estimate)) / length(ρ)
    return isapprox(normalizeddiff, 0.0, atol=atol)
  end
end

function integralsquare(func, grid::OneDGrid)
    sum(abs2.(func))*grid.dx
end

function integralsquare(func, grid::TwoDGrid)
    sum(abs2.(func))*grid.dx*grid.dy
end

function test_parsevalsum(func, grid; realvalued=true)
  # Compute integral in physical space
  integral = integralsquare(func, grid)
  # Compute integral in wavenumber space using Parseval's theorem
  if realvalued==true;      funch = rfft(func)
  elseif realvalued==false; funch = fft(func)
  end
  parsevalsum = FourierFlows.parsevalsum(abs2.(funch), grid)
  isapprox(integral, parsevalsum; rtol=rtol_utils)
end

function test_parsevalsum2(func, grid; realvalued=true)
  # Compute integral in physical space
  integral = integralsquare(func, grid)
  # Compute integral in wavenumber space using Parseval's theorem
  if realvalued==true;      funch = rfft(func)
  elseif realvalued==false; funch = fft(func)
  end
  parsevalsum2 = FourierFlows.parsevalsum2(funch, grid)
  isapprox(integral, parsevalsum2; rtol=rtol_utils)
end

"""
Compute the jacobian J(a, b) and compare the result with `analytic`. Use `atol` for the comparison
to ensure validity when `analytic=0`.
"""
test_jacobian(a, b, analytic, grid) = isapprox(FourierFlows.jacobian(a, b, grid), analytic;
                                               atol=grid.nx*grid.ny*10*eps())

"""
Test the zeros macro.
"""
function test_zeros(T=Float64, dims=(13, 45))
  a1, b1 = zeros(T, dims), zeros(T, dims)
  @zeros T dims a2 b2
  a1 == a2 && b1 == b2
end

abstract type TestVars <: AbstractVars end

physicalvars = [:a, :b]
 fouriervars = [:ah, :bh]

varspecs = cat(FourierFlows.getfieldspecs(physicalvars, :(Array{T,2})),
  FourierFlows.getfieldspecs(fouriervars, :(Array{Complex{T},2})), dims=1)

eval(FourierFlows.varsexpression(:VarsFields, physicalvars, fouriervars))
eval(FourierFlows.varsexpression(:VarsFieldsParent, physicalvars, fouriervars; parent=:TestVars))

eval(FourierFlows.varsexpression(:VarsSpecs, varspecs; typeparams=:T))
eval(FourierFlows.varsexpression(:VarsSpecsParent, varspecs; parent=:TestVars, typeparams=:T))

function test_varsexpression_fields(g::AbstractGrid{T}) where T
  @zeros T (g.nx, g.ny) a b
  @zeros Complex{T} (g.nkr, g.nl) ah bh
  v1 = VarsFields(a, b, ah, bh)

  (
     typeof(v1.a) == Array{T,2} &&
    typeof(v1.ah) == Array{Complex{T},2} &&
       size(v1.a) == size(v1.b) &&
      size(v1.ah) == size(v1.bh) &&
       size(v1.a) == (g.nx, g.ny) &&
      size(v1.ah) == (g.nkr, g.nl)
  )
end

function test_varsexpression_fields_parent(g)
  @zeros T (g.nx, g.ny) a b
  @zeros Complex{T} (g.nkr, g.nl) ah bh
  v2 = VarsFieldsParent(a, b, ah, bh)
  (
     typeof(v2.a) == Array{T,2} &&
    typeof(v2.ah) == Array{Complex{T},2} &&
       size(v2.a) == size(v2.b) &&
      size(v2.ah) == size(v2.bh) &&
       size(v2.a) == (g.nx, g.ny) &&
      size(v2.ah) == (g.nkr, g.nl)
  )
end

function test_varsexpression_specs(g::AbstractGrid{T}) where T
  @zeros T (g.nx, g.ny) a b
  @zeros Complex{T} (g.nkr, g.nl) ah bh
  v1 = VarsSpecs(a, b, ah, bh)

  (
     typeof(v1.a) == Array{T,2} &&
    typeof(v1.ah) == Array{Complex{T},2} &&
       size(v1.a) == size(v1.b) &&
      size(v1.ah) == size(v1.bh) &&
       size(v1.a) == (g.nx, g.ny) &&
      size(v1.ah) == (g.nkr, g.nl)
  )
end

function test_varsexpression_specs_parent(g)
  @zeros T (g.nx, g.ny) a b
  @zeros Complex{T} (g.nkr, g.nl) ah bh
  e2 = VarsSpecsParent(a, b, ah, bh)
  (
     typeof(v2.a) == Array{T,2} &&
    typeof(v2.ah) == Array{Complex{T},2} &&
       size(v2.a) == size(v2.b) &&
      size(v2.ah) == size(v2.bh) &&
       size(v2.a) == (g.nx, g.ny) &&
      size(v2.ah) == (g.nkr, g.nl)
  )
end

function test_supersize()
  a = rand(16, 16)
  dimsa = size(a)

  b = [rand(1), rand(3, 34)]
  dimsb = ((1,), (3, 34))

  dimsa == supersize(a) && dimsb == supersize(b)
end

function test_arraytype(dev::Device)
  dev==CPU() ? ArrayType(dev)==Array : ArrayType(dev)==CuArray
end

function test_arraytypeTdim(dev::Device, T=Float64, dim=1)
  dev==CPU() ? ArrayType(dev, T, dim)<:Array{T, dim} : ArrayType(dev, T, dim)<:CuArray{T, dim}
end


text:
@testset "CUDA" begin
    @testset "DeepONet" begin
        batch_size = 2
        a = [0.83541104, 0.83479851, 0.83404712, 0.83315711, 0.83212979, 0.83096755,
             0.82967374, 0.82825263, 0.82670928, 0.82504949, 0.82327962, 0.82140651,
             0.81943734, 0.81737952, 0.8152405, 0.81302771]
        a = repeat(a, outer=(1, batch_size)) |> gpu
        sensors = collect(range(0, 1, length=16)')
        sensors = repeat(sensors, outer=(batch_size, 1)) |> gpu
        model = DeepONet((16, 22, 30), (2, 16, 24, 30), σ, tanh;
            init_branch=Flux.glorot_normal, bias_trunk=false) |> gpu
        y = model(a, sensors)
        @test size(y) == (batch_size, 16)

        mgrad = Flux.Zygote.gradient(() -> sum(model(a, sensors)), Flux.params(model))
        @test length(mgrad.grads) == 9
    end
end

ext: jl


text:
#Computes the LJ potential energy of the system using linked lists

function u_total_ll(p,Data,first_atom,next_atom,nt)
  ut = 0.
  for iat in 1:Data.N
    icell = trunc(Int64,p[iat][1]/Data.cutoff) + 1
    jcell = trunc(Int64,p[iat][2]/Data.cutoff) + 1
    for i in icell-1:icell+1, j in jcell-1:jcell+1
      iw, jw = wrap_cell(nt,i,j)
      jat = first_atom[iw,jw]
      while jat > 0
        if jat > iat
          dist = distpbc(p[jat],p[iat],Data.side)
          if dist <= Data.cutoff
            ut += upair(dist,Data)
          end
        end
       jat = next_atom[jat]
      end
    end
  end
 return ut
end


text:
# unit test he SVRG solver
include("../src/DMD_Util.jl")
include("../src/DMD_Type.jl")
include("../src/DMD_Func.jl")
include("../src/DMD_SVRG.jl")
include("../src/DMD_BFGS.jl")
include("../src/DMD_Loss.jl")

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 500;    # temporal dimension
n = 100;    # spatial dimension
k = 3;      # number of modes
T = Float64;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);

# additive noise
sigma = 1.0e-4 # background noise level
lds = 0.1 # larger deviation size
ldf = 0.05 # frequency of larger deviations
X1 = X .+ randn(real(eltype(X)),size(X))*sigma .+ (rand(real(eltype(X)),size(X)) .< ldf)*lds
X2 = copy(X1)
# specify the loss
kappa = lds/5.0
lossFunc1(r) = huber_func(r,kappa)
lossGrad1(r) = huber_grad!(r,kappa)
lossFunc2(r) = l2_func(r)
lossGrad2(r) = l2_grad!(r)

#
# Create object
#------------------------------------------------------------------------------
params1 = DMDParams(k, X1, t, lossFunc1, lossGrad1);
params2 = DMDParams(k, X2, t, lossFunc2, lossGrad2);
params2B = DMDParams(k, X2, t, lossFunc2, lossGrad2);
Random.randn!(params1.ar);
copyto!(params2.ar,params1.ar);

#
# Apply Solver
#------------------------------------------------------------------------------
tau = 10;
eta = 1e1;
itm = 2000;
tol = 1e-6;
ptf = 1000;

function prox_stab(ar)
    # project for numerical stability
    k = length(ar) >> 1;
    @inbounds for i = 1:k
        ar[i<<1-1] = min(0.0, ar[i<<1-1]);
    end
end

opts1 = DMD_SVRG_Options(tau, eta, params1, itm=itm, tol=tol, ptf=ptf,
                         prox=prox_stab);
opts2 = DMD_SVRG_Options(tau, eta, params2, itm=itm, tol=tol, ptf=ptf,
                         prox=prox_stab);

sigma   = 1e-3;
itm = 1000;
tol = 1e-6;
ptf = 10;
opts2B = DMD_BFGS_Options(sigma, params2B, itm=itm, tol=tol, ptf=ptf);

@time obj_his1, err_his1 = solveDMD_withSVRG(params1, opts1);
@time obj_his2, err_his2 = solveDMD_withSVRG(params2, opts2);

@time obj_his2B, err_his2B = solveDMD_withBFGS(params2B, opts2B);

@printf "%7.4e abs err evals w/ Huber\n" besterrperm(params1.a,at)
@printf "%7.4e abs err evals w/ L2\n" besterrperm(params2.a,at)
@printf "%7.4e abs err evals w/ L2 (BFGS)\n" besterrperm(params2.a,at)


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(AcuteML.render2file)),NamedTuple{(:id, :age, :field, :GPA, :courses), Tuple{Int64, Int64, String, Float64, Vector{String}}},typeof(render2file),String,Bool})
    Base.precompile(Tuple{typeof(AcuteML.aml_create),Expr,Vector{Union{Expr, Symbol}},Vector{Any},Vector{Union{Expr, Symbol, Type}},Vector{Union{Expr, Symbol}},Vector{Union{Missing, String}},Vector{Union{Missing, Function, Symbol}},Vector{Union{Missing, Type}},String,Type,Array{Union{Missing, Function, Symbol}, 0},Bool,Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Expr})
    Base.precompile(Tuple{typeof(AcuteML.aml_create),Expr,Vector{Union{Expr, Symbol}},Vector{Any},Vector{Union{Expr, Symbol, Type}},Vector{Union{Expr, Symbol}},Vector{Union{Missing, String}},Vector{Union{Missing, Function, Symbol}},Vector{Union{Missing, Type}},String,Type,Array{Union{Missing, Function, Symbol}, 0},Bool,Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Symbol})
    Base.precompile(Tuple{typeof(AcuteML.aml_parse),Expr})
    Base.precompile(Tuple{typeof(AcuteML.findalltext),UnitRange{Int64},Node})
    Base.precompile(Tuple{typeof(AcuteML.findfirsttext),Int64,Node})
    Base.precompile(Tuple{typeof(AcuteML.multiString),Vector{String}})
    Base.precompile(Tuple{typeof(addnode!),Document,String,Vector{String},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Dict{String, String},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Vector{Any},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Vector{Float64},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Vector{Int64},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Vector{String},Type{AbsText}})
    Base.precompile(Tuple{typeof(createnode),Type{AbsHtml},String})
    Base.precompile(Tuple{typeof(findcontent),String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Vector{Float64}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Vector{Int64}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Vector{String}},String,Node,Type{AbsText}})
    Base.precompile(Tuple{typeof(get),Node,String,Nothing})
    Base.precompile(Tuple{typeof(newTemplate),String,Symbol})
    Base.precompile(Tuple{typeof(updatecontent!),Dict{String, String},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Node,Type{AbsAttribute}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),String,String,Node,Type{AbsText}})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{Any},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{Any},Vector{Node},Node,Type})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{Float64},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{Float64},Vector{Node},Node,Type})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{Int64},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{Int64},Vector{Node},Node,Type})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{String},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{String},String,Node,Type{AbsText}})
end

max_stars_repo_name: aminya/AML

max_forks_repo_name: aminya/AML


text:
# ---------------------------- Convenience functions ---------------------------------------------

#=
bounds and inequality constraints (dl <= d <= du)

m - size of c
p - size of d

increases the number of variables to (n + p) and the number of equality constraints
to (m + p) via the introduction of slack variables

=#
function optimize(f, c!, d!, dl, du, x0::Vector{Float64}, xl, xu, m::Int64, p::Int64, param::LFPSQPParams=LFPSQPParams())
	# if no inequality constraints
	if isnothing(d!) || p == 0
		return optimize(f, c!, x0, xl, xu, m, param)
	end

	if !(length(dl) == length(du) == p)
		error("Bound vectors dl and du must be of size p")
	end

	# fill auxiliary x0 with x0 and the initial values of the slack variables
	n = length(x0)

	x0_aux = similar(x0, n + p)
	x0_aux[1:n] .= x0
	d!(view(x0_aux, n+1:n+p), x0)

	xl_aux = similar(xl, n + p)
	xl_aux[1:n] .= xl
	xl_aux[n+1:n+p] .= dl

	xu_aux = similar(xu, n + p)
	xu_aux[1:n] .= xu
	xu_aux[n+1:n+p] .= du

	function f_aux(x)
		return f(view(x, 1:n))
	end

	function c_aux!(cval, x)
		if m > 0
			c!(view(cval, 1:m), view(x, 1:n))
		end

		d!(view(cval, m+1:m+p), view(x, 1:n))
		cval[m+1:m+p] .-= view(x, n+1:n+p)

		return cval
	end

	# generate first-order functions
	grad! = generate_gradient(f_aux, x0_aux)
	jac! = generate_jacobian(c_aux!, x0_aux, m + p)

	# generate Lagrangian Hessian function
	x0_aux_dual = zeros(Dual{nothing,Float64,1}, length(x0_aux))
	x0_aux_dual .= x0_aux

	grad_dual! = generate_gradient(f_aux, x0_aux_dual)
	jac_dual! = generate_jacobian(c_aux!, x0_aux_dual, m + p)

	hess_lag_vec! = generate_hess_lag_vec(grad_dual!, jac_dual!, x0_aux, m + p)

	# truncate output
	x, obj_values, λ_kkt, term_info = optimize(f_aux, grad!, c_aux!, jac!, hess_lag_vec!, x0_aux, xl_aux, xu_aux, m + p, param)
	trunc_x = x[1:n]

	return trunc_x, obj_values, λ_kkt, term_info
end

#=
bounds and inequality constraints (d <= 0)

m - size of c
p - size of d

increases the number of variables to (n + p) and the number of equality constraints
to (m + p) via the introduction of slack variables

=#
function optimize(f, c!, d!, x0::Vector{Float64}, xl, xu, m::Int64, p::Int64, param::LFPSQPParams=LFPSQPParams())
	return optimize(f, c!, d!, -Inf.*ones(p), zeros(p), x0, xl, xu, m, p, param)
end

# bounds and equality constraints
function optimize(f, c!, x0::Vector{Float64}, xl, xu, m::Int64, param::LFPSQPParams=LFPSQPParams())

	# generate first-order functions
	grad! = generate_gradient(f, x0)
	jac! = (m == 0) ? nothing : generate_jacobian(c!, x0, m)

	# generate Lagrangian Hessian function
	x0dual = zeros(Dual{nothing,Float64,1}, length(x0))
	x0dual .= x0

	grad_dual! = generate_gradient(f, x0dual)
	jac_dual! = (m == 0) ? nothing : generate_jacobian(c!, x0dual, m)

	hess_lag_vec! = generate_hess_lag_vec(grad_dual!, jac_dual!, x0, m)

	optimize(f, grad!, c!, jac!, hess_lag_vec!, x0, xl, xu, m, param)
end

# no bounds on x
function optimize(f, c!, x0::Vector{Float64}, m::Int64, param::LFPSQPParams=LFPSQPParams())
	optimize(f, c!, x0, nothing, nothing, m, param)
end

# unconstrained
function optimize(f, x0::Vector{Float64}, param::LFPSQPParams=LFPSQPParams())
	optimize(f, nothing, x0, nothing, nothing, 0, param)
end

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

function optimize(f, grad!, c!, jac!, hess_lag_vec!, x0::Vector{Float64}, xl, xu, m::Int64, param::LFPSQPParams)
	#= perform constrained optimization of the Helfrich energy

	INPUT
	f - objective function
	grad! - gradient of objective function, called as grad!(g, x)
	c! - constraint function, called as c!(cval, x)
	jac! - jacobian of constraint function, called as jac!(Jc, cval, x)
	hess_lag_vec! - function to calculate the action of the Lagrangian Hessian
	x0 - initial guess
	m - number of constraints
	param - algorithmic parameters in a struct LFPSQPParams

	OUTPUT
	x - the point to which the algorithm converged
	obj_values - an array of f evaluated at all steps taken
	kkt_values - an array containing estimates of the kkt multipliers (if they exist)
	term_info - a TerminationInfo struct containing termination criteria

	=#

	n = length(x0)

	# --------------------- check input data --------------------------------

	if !isnothing(xl) && !isnothing(xu)
		if !(length(xl) == length(xu) == length(x0))
			error("xl, xu, and x0 must all be the same length")
		end
	end

	# set up inequality storage and check bounds
	if (isnothing(xl) && isnothing(xu)) || (all(xl .== -Inf) && all(xu .== Inf))
		ineq = false

		f_aug = f

		ineqdata = InequalityData()
	else
		ineq = true

		if any(xl .> xu)
			error("Infeasible: lower bounds cannot be greater than upper bounds")
		end

		f_aug = x -> f(view(x, 1:n))

		PJct = Array{Float64}(undef, 2*n, m)
		λy_kkt = zeros(n)

		ineqdata = InequalityData(xl, xu)
	end

	n_ineq = ineq ? 2*n : n
	m_ineq = ineq ? m+n : m

	# output arrays
	x = Array{Float64}(undef, n_ineq)
	x[1:n] .= x0

	# find initial values for y
	if ineq
		generate_initial_y!(x, ineqdata)
	end

	obj_values = zeros(0)

	# ---------- storage for steps -----------------
	xtilde = similar(x)					# temporary step
	xnew = similar(x)
	Jc = Array{Float64}(undef, m, n) 	# Jacobian of constraints
	Jct = Array{Float64}(undef, n, m)	# gradient of constraints (Jacobian transpose)
	g = zeros(n_ineq)					# gradient of objective function (0 w.r.t. y if inequalities present)
	d = Array{Float64}(undef, n_ineq)	# step to take

	tmp_n = Array{Float64}(undef, n_ineq)	# temp vectors
	tmp_m = Array{Float64}(undef, m_ineq)

	cval = zeros(m)						# to store constraint values
	λ_kkt = zeros(m)
	term_cond::TerminationCondition = f_tol

	# storage and work vectors for SVD
	U = Array{Float64}(undef, n_ineq, m)
	Σ = Vector{Float64}(undef, m)
	Vt = Array{Float64}(undef, m, m)
	svdwork = Vector{Float64}(undef, 1)

	if m > 0
		ksvd!(ineq ? PJct : Jct, U, Σ, Vt, svdwork, true)		# calculate work vector
	end

	newton_d = Array{Float64}(undef, n_ineq)				# truncated Newton steps
	newton_Δλ = Array{Float64}(undef, m_ineq)
	newton_b2 = zeros(m_ineq)
	projcgwork = ProjCGWork(n_ineq, m_ineq)
	prev_grad_norm = 0.0					# forces tol to be κ*grad_norm initially
	grad_norm = Inf

	# inequality decomposition
	if ineq
		ineqdecomp = InequalityDecomp(U, Σ, Vt, [Array{Float64}(undef, n) for i in 1:3]..., Jct, m)
	else
		ineqdecomp = InequalityDecomp(U, Σ, Vt, [zeros(0) for i in 1:3]..., Jct, m)
	end
	ineqproject = InequalityDecompProject(ineqdecomp)

	# construct Linear map to be used for Newton steps
	if ineq
		newton_map = LinearMap{Float64}((dest, src) -> augmented_hess_lag_vec!(dest, src, hess_lag_vec!, x, λ_kkt, λy_kkt, ineqdata), 2*n, ismutating=true, issymmetric=true)
	else
		newton_map = LinearMap{Float64}((dest, src) -> hess_lag_vec!(dest, src, x, λ_kkt), n, ismutating=true, issymmetric=true)
	end

	# retraction methods
	nr = NR(U, Σ, Vt, param.ϵ_c, param.maxiter_retract, NRWork(m), ineq, ineqdata)
	pp = ProjPenalty(jac!, U, Σ, Vt, m, param.μ0, param.ϵ_c, param.maxiter_retract, param.maxiter_pcg, ProjPenaltyWork(m, n, m_ineq, n_ineq), ineq, ineqdecomp, ineqdata)
	euc = Euclidean()
	yr = YRetract(ineqdata)

	# linesearch work structs
	armijo_work = ArmijoWork(n_ineq)
	exact_work = ExactLinesearchWork(n_ineq)

	# ----------------------- Calculate gradient -------------------------
	i = 0
	f_diff = Inf
	step_diff = Inf
	kkt_diff = Inf

	fval = f_aug(x)
	append!(obj_values, fval)

	m > 0 && c!(cval, view(x, 1:n))

	param.disp == iter && print_iter_header()
	param.disp == iter && print_first_line(fval, norm(cval, Inf))

	while true
		# calculate gradient at current point
		grad!(view(g, 1:n), view(x, 1:n))

		# calculate step and add random noise (if set)
		d .= -1.0 .* g

		if param.β > 0
			randn!(tmp_n)

			if param.t_β > 0
				# if there should be some ramp down of random noise
				axpy!(param.β*max(1 - i/param.t_β, 0.0), tmp_n, d)
			else
				axpy!(param.β, tmp_n, d)
			end
		end

		if ineq
			# calculate gradients of inequality constraints
			inequality_gradient!(ineqdecomp, x, ineqdata)
		end

		rank = m
		if m > 0
			# calculate Jacobian at current point and find SVD
			jac!(Jc, cval, view(x, 1:n))
			transpose!(Jct, Jc)

			if ineq
				# project Jacobian orthogonal to inequality constraints
				PJct[1:n, :] .= (1.0 .- ineqdecomp.Dx.*ineqdecomp.Dx) .* Jct
				PJct[n+1:2*n, :] .= -1.0 .* ineqdecomp.Dy.*ineqdecomp.Dx .* Jct

				ksvd!(PJct, U, Σ, Vt, svdwork)		# O(Nm^2)
			else
				ksvd!(Jct, U, Σ, Vt, svdwork)		# O(Nm^2)
			end

			# find rank of constraint Jacobian
			for (j, a) in enumerate(Σ)
				if a < param.ϵ_rank
					rank = j - 1
					break
				end
			end

			# project step onto tangent space if no inequalities present
			if !ineq
				kgemv!('T', rank, 1.0, U, d, 0.0, tmp_m)		# tmp_m = U' d
				kgemv!('N', rank, -1.0, U, tmp_m, 1.0, d)		# d = d - U*tmp_m
			end
		end

		# project step onto tangent space for inequalities
		if ineq
			# update rank in inequality decomposition
			ineqdecomp.rank = rank

			mul!(tmp_m, ineqproject', d)				# tmp_m = Q'd
			mul!(d, ineqproject, tmp_m, -1.0, 1.0)		# d = d - Q*tmp_m
		end

		kkt_diff = norm(d, Inf)

		# update rank in ProjPenalty
		pp.rank = rank

		steptype = 0
		tn_iter = 0
		tn_res = 0.0

		# calculate λ_kkt
		if ineq
			calculate_λ_kkt!(λ_kkt, λy_kkt, tmp_m, ineqdecomp)
		elseif m > 0
			# tmp_m contains U'd = -U'g
			@inbounds @simd for j = 1:rank
				tmp_m[j] /= Σ[j]
			end
			@inbounds @simd for j = rank+1:m
				tmp_m[j] = 0.0
			end

			gemv!('T', 1.0, Vt, tmp_m, 0.0, λ_kkt) # λ = -V Σ^{-1} U' g
		end

		# -------------------- Check for termination conditions --------------

		if f_diff <= param.ϵ_f
			term_cond = f_tol
			break
		elseif step_diff <= param.ϵ_x
			term_cond = x_tol
			break
		elseif i >= param.maxiter
			term_cond = max_iter
			break
		elseif kkt_diff <= param.ϵ_kkt
			term_cond = kkt_tol
			break
		end

		# -------------------- Truncated Newton step -------------------------

		if param.do_newton
			# prepare views
			if ineq
				Qview = ineqproject
				newton_b2_view = view(newton_b2, 1:n+rank)
			else
				Qview = view(U, :, 1:rank)
				newton_b2_view = view(newton_b2, 1:rank)
			end

			# calculate tolerance for truncated newton step
			grad_norm = norm(d)
			tol = param.tn_κ*min(1, (grad_norm/prev_grad_norm))*grad_norm

			prev_grad_norm = grad_norm  	# update prev_grad_norm for next time

			# take truncated newton step using ProjCG
			tn_iter, tn_res = projcg!(newton_d, newton_Δλ, newton_map, Qview, d, newton_b2_view,
				tol=tol, maxit=param.tn_maxiter, work=projcgwork)

			# choose Newton direction if in gradient direction
			if dot(newton_d, d) > 0.0
				d .= newton_d
				steptype = 1
			end
		end

		# ---------------------- Linesearch ------------------------

		# determine retraction method to use
		if m > 0
			if rank == m && !param.do_project_retract
				retract_method = nr
				mtype = 0
			else
				retract_method = pp
				mtype = 1
			end
		else
			if ineq
				retract_method = yr
				mtype = 0
			else
				retract_method = euc
				mtype = 0
			end
		end

		if param.linesearch == armijo || param.disable_linesearch
			# note - don't need to project g onto tangent manifold for dot product since d is guaranteed to live on tangent manifold
			flag, iter1, iter2, newf, f_diff, step_diff, α = armijo!(xnew, x, n, d, g, f_aug, fval, retract_method, cval, c!, param, armijo_work)
		else
			flag, iter1, iter2, newf, f_diff, step_diff, α = exact_linesearch!(xnew, x, n, d, f_aug, fval, retract_method, cval, c!, param, exact_work)
		end

		# ----------------------- update x and function values ---------------------------
		x .= xnew
		fval = newf
		append!(obj_values, fval)
		param.disp == iter && print_iter(i+1, fval, norm(cval,Inf), f_diff, step_diff, steptype, tn_iter, tn_res, mtype, iter1, iter2, α, flag)

		i += 1

		# callback
		if param.callback != nothing && mod(i, param.callback_period) == 0
			param.callback(i, x)
		end
	end

	if i == param.maxiter
		@warn "Maximum # of outer iterations reached"
	end

	return x[1:n], obj_values, λ_kkt, TerminationInfo(term_cond, f_diff, step_diff, kkt_diff, i)
end

@inline function print_iter_header()
	@printf("   step |          f     ||c||      |Δf|    ||Δx||  |   S iter      res  |   M   iter  (pcg)  |        α  flag\n")
	@printf("--------------------------------------------------------------------------------------------------------------\n")
end

@inline function print_first_line(fval, normc)
	@printf("      0 | %10.3e  %8.1e                      |                    |                    |               \n", fval, normc)
end

@inline function print_iter(i, fval, normc, fstep, normx, steptype, tn_iter, tn_res, methodtype, iter1, iter2, α, flag)
	# print out iteration information
	if methodtype == 0
		method = "NR"
	else
		method = "PP"
	end

	if steptype == 0
		stepname = "GD"
	else
		stepname = "TN"
	end

	@printf("%7d | %10.3e  %8.1e  %8.1e  %8.1e  |  %s %4d %8.1e  |  %s %6d %6d  | %8.1e  %4d\n",
		i, fval, normc, fstep, normx, stepname, tn_iter, tn_res, method, iter1, iter2, α, flag)

	flush(Base.stdout)
end


text:
const DEFAULT_N_SITES = Dict{String,String}(
    "e" => "n_e",
    "e1" => "n_e1",
    "e2" => "n_e2",
    "H" => "n_H"
)

"""
    ParamOptions(;kwargs...)

Struct containing all the options related to parameter parsing:
* `userlocations::Vector{String} = String[]`: List of used-defined locations to search.
* `usergrouplocations::Vector{String} = String[]`: List of used-defined group locations to search.
* `asymmetricparams::Vector{String} = String[]`: List of pair or association parameters that follow that `param[i,j] ≠ param[j,i]`
* `ignore_headers::Vector{String} =  ["dipprnumber", "smiles"]`: List of ignored headers.
* `ignore_missing_singleparams::Bool=false`: List of parameters where checking for missing single parameter values are ignored.
* `verbose::Bool = false`: If `true`, show all operations done by `getparams` displayed in the terminal. this includes the warnings emmited by `CSV.jl`
* `species_columnreference::String ="species"`: column name to check for components. in pair and association params, it will check for `#species#1` and `#species#2`, where `#species#` is the value of this option.
* `site_columnreference::String ="site"`: column name to check for sites in association params, it will check for `#site#1` and `#site#2`, where `#site#` is the value of this option.
* `group_columnreference::String = "groups"`: column name to check for groups in group data.
* `normalisecomponents::Bool = true`: If `true`, performs normalization of strings, on the CSV and input components
* `n_sites_columns::Dict{String,String} = Dict( "e" => "n_e","e1" => "n_e1","e2" => "n_e2","H" => "n_H")`: dictionary to look number of sites. the number of sites is stored as columns in a single parameter csv file. for example, the number of sites of name `e` will be looked on the column `n_e`
* `return_sites::Bool = true`: If set to false, association params will be ignored and sites will not be created, even if they exist in the list of locations.
* `component_delimiter::String = "~|~"`: When there are multiple component names to match, seperate them by this delimiter.
"""
Base.@kwdef struct ParamOptions
    userlocations::Vector{String} = String[]
    usergrouplocations::Vector{String} = String[]
    asymmetricparams::Vector{String}= String[]
    ignore_missing_singleparams::Vector{String} = String[]
    ignore_headers::Vector{String} = ["dipprnumber", "smiles"]
    verbose::Bool = false
    species_columnreference::String ="species"
    source_columnreference::String = "source"
    site_columnreference::String = "site"
    group_columnreference::String = "groups"
    normalisecomponents::Bool = true
    n_sites_columns::Dict{String,String} = DEFAULT_N_SITES
    return_sites::Bool = true
    component_delimiter::String = "~|~"
end

export ParamOptions


text:
using DSP

"""
    y = undbp(x)

Convert `x` from dB to a power
"""
function undbp(x)
    return db2pow.(x)
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
# AUTO GENERATED FILE - DO NOT EDIT

export anchor

"""
    anchor(;kwargs...)
    anchor(children::Any;kwargs...)
    anchor(children_maker::Function;kwargs...)

An Anchor component.
Display links with theme styles. For more information, see: https://mantine.dev/core/anchor/
Keyword arguments:
- `children` (String; optional): Primary content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `align` (a value equal to: "left", "right", "center"; optional): Sets text-align css property
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Text color from theme
- `gradient` (optional): Controls gradient settings in gradient variant only. gradient has the following type: lists containing elements 'from', 'to', 'deg'.
Those elements have the following types:
  - `from` (String; required)
  - `to` (String; required)
  - `deg` (Real; optional)
- `href` (String; optional): href
- `inherit` (Bool; optional): Inherit font properties from parent element
- `inline` (Bool; optional): Sets line-height to 1 for centering
- `lineClamp` (Real; optional): CSS -webkit-line-clamp property
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined font-size from theme.fontSizes
- `style` (Dict; optional): Inline style override
- `target` (a value equal to: "_blank", "_self"; optional): Target
- `transform` (a value equal to: "capitalize", "uppercase", "lowercase"; optional): Sets text-transform css property
- `variant` (a value equal to: "link", "gradient", "text"; optional): Link or text variant
- `weight` (a value equal to: "normal", "bold", "bolder", "lighter", "initial", "inherit" | Real; optional): Sets font-weight css property
"""
function anchor(; kwargs...)
        available_props = Symbol[:children, :id, :align, :className, :color, :gradient, :href, :inherit, :inline, :lineClamp, :size, :style, :target, :transform, :variant, :weight]
        wild_props = Symbol[]
        return Component("anchor", "Anchor", "dash_mantine_components", available_props, wild_props; kwargs...)
end

anchor(children::Any; kwargs...) = anchor(;kwargs..., children = children)
anchor(children_maker::Function; kwargs...) = anchor(children_maker(); kwargs...)


text:
using Gnuplot, RDatasets
cars = dataset("datasets", "mtcars")
byCat = cars.Cyl
categ = unique(byCat)
# aesthetics
@gp linetypes(:Set1_5) :-
@gp :- xlab ="Horse Power" ylab = "Miles/(US) gallon" "set auto fix"
@gp :- title = "Categorical Data" :-
@gp :- "set offsets graph .05, graph .05, graph .05, graph .05" :-
@gp :- "set key t r font ',12' tit 'Cylinders' box" :-
# the actual plot
for c in categ
    indc = findall(x->x == c, byCat)
    @gp :-  cars.HP[indc] cars.MPG[indc]  "w p tit '$(c)' pt 7 ps 1.4 "
end
save(term ="pngcairo size 600,400", output ="Categorical_cars.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
import Test: @test, @testset, @test_throws
import Luna: Output, Processing
using EllipsisNotation

@testset "HDF5" begin
    import HDF5
    import Luna: Utils
    fpath = joinpath(tempname(), "test.h5")
    isfile(fpath) && rm(fpath)
    shape = (1024, 4, 2)
    n = 11
    stat = randn()
    statsfun(y, t, dt) = Dict("stat" => stat, "stat2d" => [stat, stat])
    t0 = 0
    t1 = 10
    t = collect(range(t0, stop=t1, length=n))
    ω = randn((1024,))
    wd = dirname(@__FILE__)
    gitc = Utils.git_commit()
    o = Output.HDF5Output(fpath, t0, t1, n, statsfun; yname="y", tname="t")
    extra = Dict()
    extra["ω"] = ω
    extra["git_commit"] = gitc
    o(extra)
    meta = Dict()
    meta["meta1"] = 100
    meta["meta2"] = "src"
    o(meta, meta=true)
    y0 = randn(ComplexF64, shape)
    y(t) = y0
    for (ii, ti) in enumerate(t)
        o(y0, ti, 0, y)
    end
    @test o(extra, force=true) === nothing
    HDF5.h5open(fpath, "r") do file
        @test all(read(file["t"]) == t)
        global yr = read(file["y"])
        @test all([all(yr[:, :, :, ii] == y0) for ii=1:n])
        @test all(ω == read(file["ω"]))
        @test gitc == read(file["git_commit"])
        @test all(read(file["stats"]["stat"]) .== stat)
        @test all(read(file["stats"]["stat2d"]) .== stat)
        @test Utils.git_commit() == read(file["meta"]["git_commit"])
        # Need to strip out date from sourcecode to compare
        src = read(file["meta"]["sourcecode"])
        @test split(Utils.sourcecode(), '\n')[2:end] == split(src, '\n')[2:end]
        @test 100 == read(file["meta"]["meta1"])
        @test "src" == read(file["meta"]["meta2"])
    end
    @test all(yr .== o["y"])
    @test 100 == o["meta"]["meta1"]
    rm(fpath)
    rm(splitdir(fpath)[1], force=true)
end

@testset "Memory" begin
    import Luna: Utils
    shape = (1024, 4, 2)
    n = 11
    stat = randn()
    statsfun(y, t, dt) = Dict("stat" => stat, "stat2d" => [stat, stat])
    t0 = 0
    t1 = 10
    t = collect(range(t0, stop=t1, length=n))
    ω = randn((1024,))
    wd = dirname(@__FILE__)
    gitc = Utils.git_commit()
    o = Output.MemoryOutput(t0, t1, n, statsfun, yname="y", tname="t")
    extra = Dict()
    extra["ω"] = ω
    extra["git_commit"] = gitc
    o(extra)
    meta = Dict()
    meta["meta1"] = 100
    meta["meta2"] = "src"
    o(meta, meta=true)
    y0 = randn(ComplexF64, shape)
    y(t) = y0
    for (ii, ti) in enumerate(t)
        o(y0, ti, 0, y)
    end
    @test all(o.data["t"] == t)
    @test all([all(o.data["y"][:, :, :, ii] == y0) for ii=1:n])
    @test all(ω == o.data["ω"])
    @test gitc == o.data["git_commit"]
    @test_throws ErrorException o(extra)
    @test o(extra, force=true) === nothing
    @test all(o.data["stats"]["stat"] .== stat)
    @test all(o.data["stats"]["stat2d"] .== stat)
    @test Utils.git_commit() == o.data["meta"]["git_commit"]
    # Need to strip out date from sourcecode to compare
    src = o.data["meta"]["sourcecode"]
    @test split(Utils.sourcecode(), '\n')[2:end] == split(src, '\n')[2:end]
    @test 100 == o.data["meta"]["meta1"]
    @test "src" == o.data["meta"]["meta2"]
end

dirpath = tempname()
fpath = joinpath(dirpath, "test.h5")
fpath_comp = joinpath(dirpath, "test_comp.h5")
@testset "HDF5 vs Memory" begin
    using Luna
    import FFTW
    import HDF5
    import LinearAlgebra: norm

    a = 13e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)
    m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
    aeff = let m=m
        z -> Modes.Aeff(m, z=z)
    end
    energyfun, energyfunω = Fields.energyfuncs(grid)
    densityfun = let dens0=PhysData.density(gas, pres)
        z -> dens0
    end
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)

    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)
    statsfun = Stats.collect_stats(grid, Eω,
                                   Stats.ω0(grid),
                                   Stats.energy(grid, energyfunω))
    hdf5 = Output.HDF5Output(fpath, 0, grid.zmax, 51, statsfun)
    hdf5c = Output.HDF5Output(fpath_comp, 0, grid.zmax, 51, statsfun,
                              compression=true)
    mem = Output.MemoryOutput(0, grid.zmax, 51, statsfun)
    function outfun(args...; kwargs...)
        hdf5(args...; kwargs...)
        hdf5c(args...; kwargs...)
        mem(args...; kwargs...)
    end
    for o in (hdf5, hdf5c, mem)
        o(Dict("λ0" => λ0))
        o("τ", τ)
    end
    Luna.run(Eω, grid, linop, transform, FT, outfun, status_period=10)
    HDF5.h5open(hdf5.fpath, "r") do file
        @test read(file["λ0"]) == mem.data["λ0"]
        Eω = reinterpret(ComplexF64, read(file["Eω"]))
        @test Eω == mem.data["Eω"]
        @test read(file["stats"]["ω0"]) == mem.data["stats"]["ω0"]
        @test read(file["z"]) == mem.data["z"]
        @test read(file["grid"]) == Grid.to_dict(grid)
        @test read(file["simulation_type"]["field"]) == "field-resolved"
        @test read(file["simulation_type"]["transform"]) == string(transform)
    end
    HDF5.h5open(hdf5c.fpath, "r") do file
        @test read(file["λ0"]) == mem.data["λ0"]
        Eω = reinterpret(ComplexF64, read(file["Eω"]))
        @test Eω == mem.data["Eω"]
        @test read(file["stats"]["ω0"]) == mem.data["stats"]["ω0"]
        @test read(file["z"]) == mem.data["z"]
        @test read(file["grid"]) == Grid.to_dict(grid)
        @test read(file["simulation_type"]["field"]) == "field-resolved"
        @test read(file["simulation_type"]["transform"]) == string(transform)
    end
    @test stat(hdf5.fpath).size >= stat(hdf5c.fpath).size
    # Test read-only
    o = Output.HDF5Output(fpath)
    HDF5.h5open(o.fpath, "r") do file
        @test read(file["λ0"]) == mem.data["λ0"]
        Eω = reinterpret(ComplexF64, read(file["Eω"]))
        @test Eω == mem.data["Eω"]
        @test read(file["stats"]["ω0"]) == mem.data["stats"]["ω0"]
        @test read(file["z"]) == mem.data["z"]
        @test read(file["grid"]) == Grid.to_dict(grid)
        @test read(file["simulation_type"]["field"]) == "field-resolved"
        @test read(file["simulation_type"]["transform"]) == string(transform)
    end
    # test slice reading
    @test o["Eω", :, 1] == mem["Eω"][:, 1]
    @test o["Eω", 1, :] == mem["Eω"][1, :]
    @test o["Eω", :, 1:5] == mem["Eω"][:, 1:5]
    @test o["Eω", :, [1, 2, 50]] == mem["Eω"][:, [1, 2, 50]]
    @test o["Eω", .., 1] == mem["Eω"][:, 1]
    @test o["Eω", .., [1, 2, 50]] == mem["Eω"][:, [1, 2, 50]]
    @test mem["Eω", :, 1] == mem["Eω"][:, 1]
    @test mem["Eω", 1, :] == mem["Eω"][1, :]
    @test mem["Eω", :, 1:5] == mem["Eω"][:, 1:5]
    @test mem["Eω", :, [1, 2, 50]] == mem["Eω"][:, [1, 2, 50]]
    @test mem["Eω", .., 1] == mem["Eω"][:, 1]
    @test mem["Eω", .., [1, 2, 50]] == mem["Eω"][:, [1, 2, 50]]
    # test slice reading in processing functions
    ω, Eω, zac = Processing.getEω(o, 5e-2)
    @test (ω, Eω[:, 1]) == Processing.getEω(grid, mem["Eω"][:, 51])
    @test zac[1] == 5e-2
end
rm(fpath)
rm(fpath_comp)
rm(splitdir(fpath)[1], force=true)

##
fpath = joinpath(homedir(), ".luna", "output_test", "test.h5")
@testset "Continuing" begin
    using Luna
    import FFTW
    import HDF5

    a = 13e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)
    m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
    aeff(z) = Modes.Aeff(m, z=z)
    energyfun, energyfunω = Fields.energyfuncs(grid)
    dens0 = PhysData.density(gas, pres)
    densityfun(z) = dens0
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)

    # Run with arbitrary error at 3 cm
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)
    statsfun = Stats.collect_stats(grid, Eω,
                                   Stats.ω0(grid),
                                   Stats.energy(grid, energyfunω))
    output = Output.HDF5Output(fpath, 0, grid.zmax, 51, statsfun)
    function stepfun(Eω, z, dz, interpolant)
        output(Eω, z, dz, interpolant)
        if z > 3e-2
            error("Oh no!")
        end
    end
    stepfun(args...; kwargs...) = output(args...; kwargs...)
    try
        Luna.run(Eω, grid, linop, transform, FT, stepfun, status_period=10, z0=0.0)
    catch
    end

    # Run again, starting from 3 cm
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)
    statsfun = Stats.collect_stats(grid, Eω,
                                   Stats.ω0(grid),
                                   Stats.energy(grid, energyfunω))
    output = Output.HDF5Output(fpath, 0, grid.zmax, 51, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)

    # Run from scratch with MemoryOutput
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)
    statsfun = Stats.collect_stats(grid, Eω,
                                   Stats.ω0(grid),
                                   Stats.energy(grid, energyfunω))
    mem = Output.MemoryOutput(0, grid.zmax, 51, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, mem, status_period=5)

    idx1 = findfirst(grid.ωwin .!= 1)
    idx2 = findlast(grid.ωwin .== 1)
    Eωm = mem["Eω"][idx1:idx2, :]
    Eω = output["Eω"][idx1:idx2, :]
    Iω = abs2.(Eω)
    Iωm = abs2.(Eωm)
    @test norm(Iω - Iωm)/norm(Iω) < 1e-7
    @test all(isapprox.(Eωm, Eω, atol=1e-4*maximum(abs.(Eωm))))
    @test all(isapprox.(output["stats"]["ω0"], mem.data["stats"]["ω0"], rtol=1e-6))
    @test all(output["stats"]["energy"] .≈ mem.data["stats"]["energy"])
    @test output["z"] == mem.data["z"]
    @test output["grid"] == Grid.to_dict(grid)
    @test output["simulation_type"]["field"] == "field-resolved"
    @test output["simulation_type"]["transform"] == string(transform)
end
rm(fpath, force=true)
rm(splitdir(fpath)[1], force=true)


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#
# ----------------------------------------------------------------------------------- #
# Function: DataDictionary
# Description: Holds simulation and model parameters as key => value pairs in a Julia Dict()
# Generated on: 2016-11-02T07:00:00.307
#
# Input arguments:
# time_start::Float64 => Simulation start time value (scalar)
# time_stop::Float64 => Simulation stop time value (scalar)
# time_step::Float64 => Simulation time step (scalar)
#
# Output arguments:
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model and simulation parameters as key => value pairs
# ----------------------------------------------------------------------------------- #
function DataDictionary(time_start::Float64,time_stop::Float64,time_step_size::Float64)

	# stoichiometric_matrix and dilution_matrix -
	stoichiometric_matrix = readdlm("./Network.dat")
	dilution_matrix = readdlm("./Dilution.dat")
	degradation_matrix = readdlm("./Degradation.dat")

	# array of gene lengths -
	# gene_coding_length_array = [
	# 	15000	;	# 1	gene_AP1
	# 	15000	;	# 2	gene_AhR
	# 	15000	;	# 3	gene_CD11b
	# 	15000	;	# 4	gene_CD14
	# 	15000	;	# 5	gene_CD38
	# 	15000	;	# 6	gene_CEBPa
	# 	15000	;	# 7	gene_E2F
	# 	15000	;	# 8	gene_EGR1
	# 	15000	;	# 9	gene_GFI1
	# 	15000	;	# 10	gene_IRF1
	# 	15000	;	# 11	gene_OCT1
	# 	15000	;	# 12	gene_OCT4
	# 	15000	;	# 13	gene_P21
	# 	15000	;	# 14	gene_P47Phox
	# 	15000	;	# 15	gene_PPARg
	# 	15000	;	# 16	gene_PU1
	# 	15000	;	# 17	gene_Trigger
	# 	15000	;	# 18	gene_cRAF
	# ]

	gene_coding_length_array = [
		10323	;	# 1	gene_AP1
		47530	;	# 2	gene_AhR
		72925	;	# 3	gene_CD11b
		8974	;	# 4	gene_CD14
		74978	;	# 5	gene_CD38
		2630	;	# 6	gene_CEBPa
		17919	;	# 7	gene_E2F
		10824	;	# 8	gene_EGR1
		13833	;	# 9	gene_GFI1
		16165	;	# 10	gene_IRF1
		206516 	;	# 11	gene_OCT1
		6356	;	# 12	gene_OCT4
		15651	;	# 13	gene_P21
		3074	;	# 14	gene_P47Phox https://www.ncbi.nlm.nih.gov/nuccore/AF003533.1   #1475	;	# 14	gene_P47Phox -- https://www.ncbi.nlm.nih.gov/nuccore/NM_000265.5
		153507	;	# 15	gene_PPARg
		40782	;	# 16	gene_PU1
		5908	;	# 17	gene_Trigger
		87571; # 18	gene_cRAF 	https://www.ncbi.nlm.nih.gov/nuccore/NG_007467.1 #3291	;	# 18	gene_cRAF --- NCBI homo sapien: https://www.ncbi.nlm.nih.gov/nuccore/NM_002880.3
	];

	# array of mRNA coding lengths -
	mRNA_coding_length_array = [
		gene_coding_length_array[1]	;	# 19	1	mRNA_gene_AP1
		gene_coding_length_array[2]	;	# 20	2	mRNA_gene_AhR
		gene_coding_length_array[3]	;	# 21	3	mRNA_gene_CD11b
		gene_coding_length_array[4]	;	# 22	4	mRNA_gene_CD14
		gene_coding_length_array[5]	;	# 23	5	mRNA_gene_CD38
		gene_coding_length_array[6]	;	# 24	6	mRNA_gene_CEBPa
		gene_coding_length_array[7]	;	# 25	7	mRNA_gene_E2F
		gene_coding_length_array[8]	;	# 26	8	mRNA_gene_EGR1
		gene_coding_length_array[9]	;	# 27	9	mRNA_gene_GFI1
		gene_coding_length_array[10]	;	# 28	10	mRNA_gene_IRF1
		gene_coding_length_array[11]	;	# 29	11	mRNA_gene_OCT1
		gene_coding_length_array[12]	;	# 30	12	mRNA_gene_OCT4
		gene_coding_length_array[13]	;	# 31	13	mRNA_gene_P21
		gene_coding_length_array[14]	;	# 32	14	mRNA_gene_P47Phox
		gene_coding_length_array[15]	;	# 33	15	mRNA_gene_PPARg
		gene_coding_length_array[16]	;	# 34	16	mRNA_gene_PU1
		gene_coding_length_array[17]	;	# 35	17	mRNA_gene_Trigger
		gene_coding_length_array[18]	;	# 36	18	mRNA_gene_cRAF
	]

	# array of mRNA coding lengths -
	# protein_coding_length_array = [
	# 	round((0.33)*mRNA_coding_length_array[1])	;	# 37	1	protein_gene_AP1
	# 	round((0.33)*mRNA_coding_length_array[2])	;	# 38	2	protein_gene_AhR
	# 	round((0.33)*mRNA_coding_length_array[3])	;	# 39	3	protein_gene_CD11b
	# 	round((0.33)*mRNA_coding_length_array[4])	;	# 40	4	protein_gene_CD14
	# 	round((0.33)*mRNA_coding_length_array[5])	;	# 41	5	protein_gene_CD38
	# 	round((0.33)*mRNA_coding_length_array[6])	;	# 42	6	protein_gene_CEBPa
	# 	round((0.33)*mRNA_coding_length_array[7])	;	# 43	7	protein_gene_E2F
	# 	round((0.33)*mRNA_coding_length_array[8])	;	# 44	8	protein_gene_EGR1
	# 	round((0.33)*mRNA_coding_length_array[9])	;	# 45	9	protein_gene_GFI1
	# 	round((0.33)*mRNA_coding_length_array[10])	;	# 46	10	protein_gene_IRF1
	# 	round((0.33)*mRNA_coding_length_array[11])	;	# 47	11	protein_gene_OCT1
	# 	round((0.33)*mRNA_coding_length_array[12])	;	# 48	12	protein_gene_OCT4
	# 	round((0.33)*mRNA_coding_length_array[13])	;	# 49	13	protein_gene_P21
	# 	round((0.33)*mRNA_coding_length_array[14])	;	# 50	14	protein_gene_P47Phox
	# 	round((0.33)*mRNA_coding_length_array[15])	;	# 51	15	protein_gene_PPARg
	# 	round((0.33)*mRNA_coding_length_array[16])	;	# 52	16	protein_gene_PU1
	# 	round((0.33)*mRNA_coding_length_array[17])	;	# 53	17	protein_gene_Trigger
	# 	round((0.33)*mRNA_coding_length_array[18])	;	# 54	18	protein_gene_cRAF
	# ]

	# array of mRNA coding lengths -
	protein_coding_length_array = [
		331	;	# 37	1	protein_gene_AP1
		848	;	# 38	2	protein_gene_AhR
		1153	;	# 39	3	protein_gene_CD11b
		375	;	# 40	4	protein_gene_CD14
		300	;	# 41	5	protein_gene_CD38
		393	;	# 42	6	protein_gene_CEBPa
		437	;	# 43	7	protein_gene_E2F
		543	;	# 44	8	protein_gene_EGR1
		422	;	# 45	9	protein_gene_GFI1
		325	;	# 46	10	protein_gene_IRF1
		741.3333333333	;	# 47	11	protein_gene_OCT1
		206.3333333	;	# 48	12	protein_gene_OCT4
		198	;	# 49	13	protein_gene_P21
		390	;	# 50	14	protein_gene_P47Phox -- https://www.ncbi.nlm.nih.gov/protein/NP_000256.4
		250	;	# 51	15	protein_gene_PPARg
		270.5	;	# 52	16	protein_gene_PU1
		420.6666667		;	# 53	17	protein_gene_Trigger
		648	;	# 54	18	protein_gene_cRAF --https://www.ncbi.nlm.nih.gov/protein/NP_002871.1
	]

	# ------------------------------------------------------------------------------------------#
	# constants (from bionumbers)       units
	# ------------------------------------------------------------------------------------------#
	cell_diameter = 12                  # mum
	number_of_rnapII = 75000            # copies/cells
	number_of_ribosome = 1e6            # copies/cells
	mRNA_half_life_TF = 2               # hrs
	protein_half_life = 10              # hrs
	doubling_time_cell = 19.5           # hrs
	max_translation_rate = 5            # aa/sec
	max_transcription_rate = 6.0        # nt/sec
	average_transcript_length = 15000   # nt
	average_protein_length = 5000       # aa
	fraction_nucleus = 0.49             # dimensionless
	av_number = 6.02e23                 # number/mol
	avg_gene_number = 2                 # number of copies of a gene
	# ------------------------------------------------------------------------------------------#
	#
	# ------------------------------------------------------------------------------------------#
	# Calculate constants using bionumber values
	# ------------------------------------------------------------------------------------------#
	# Calculate the volume (convert to L)
	V = ((1-fraction_nucleus)*(1/6)*(3.14159)*(cell_diameter)^3)*(1e-15)

	# Calculate the rnapII_concentration and ribosome_concentration
	rnapII_concentration = number_of_rnapII*(1/av_number)*(1/V)*1e9                   # nM
	ribosome_concentration = number_of_ribosome*(1/av_number)*(1/V)*1e9               # nM

	# degrdation rate constants -
	degradation_constant_mRNA = -(1/mRNA_half_life_TF)*log(0.5)                       # hr^-1
	degradation_constant_protein = -(1/protein_half_life)*log(0.5)                    # hr^-1

	# kcats for transcription and translation -
	kcat_transcription = max_transcription_rate*(3600/average_transcript_length)      # hr^-1
	kcat_translation = max_translation_rate*(3600/average_protein_length)             # hr^-1

	# Maximum specific growth rate -
	maximum_specific_growth_rate = (1/doubling_time_cell)*log(2)                      # hr^-1

	# What is the average gene concentration -
	avg_gene_concentration = avg_gene_number*(1/av_number)*(1/V)*1e9                  # nM

	# How fast do my cells die?
	death_rate_constant = 0.2*maximum_specific_growth_rate                            # hr^-1

	# Saturation constants for translation and trascription -
	saturation_transcription = 4600*(1/av_number)*(1/V)*1e9                           # nM
	saturation_translation = 100000*(1/av_number)*(1/V)*1e9                           # nM
	# -------------------------------------------------------------------------------------------#

	# initial condition array -
	initial_condition_array = [
		avg_gene_concentration	;	# 1	gene_AP1
		avg_gene_concentration	;	# 2	gene_AhR
		avg_gene_concentration	;	# 3	gene_CD11b
		avg_gene_concentration	;	# 4	gene_CD14
		avg_gene_concentration	;	# 5	gene_CD38
		avg_gene_concentration	;	# 6	gene_CEBPa
		avg_gene_concentration	;	# 7	gene_E2F
		avg_gene_concentration	;	# 8	gene_EGR1
		avg_gene_concentration	;	# 9	gene_GFI1
		avg_gene_concentration	;	# 10	gene_IRF1
		avg_gene_concentration	;	# 11	gene_OCT1
		avg_gene_concentration	;	# 12	gene_OCT4
		avg_gene_concentration	;	# 13	gene_P21
		avg_gene_concentration	;	# 14	gene_P47Phox
		avg_gene_concentration	;	# 15	gene_PPARg
		avg_gene_concentration	;	# 16	gene_PU1
		avg_gene_concentration	;	# 17	gene_Trigger
		avg_gene_concentration	;	# 18	gene_cRAF
		0.0	;	# 19	mRNA_gene_AP1
		0.0	;	# 20	mRNA_gene_AhR
		0.0	;	# 21	mRNA_gene_CD11b
		0.0	;	# 22	mRNA_gene_CD14
		0.0	;	# 23	mRNA_gene_CD38
		0.0	;	# 24	mRNA_gene_CEBPa
		0.0	;	# 25	mRNA_gene_E2F
		0.0	;	# 26	mRNA_gene_EGR1
		0.0	;	# 27	mRNA_gene_GFI1
		0.0	;	# 28	mRNA_gene_IRF1
		0.0	;	# 29	mRNA_gene_OCT1
		0.0	;	# 30	mRNA_gene_OCT4
		0.0	;	# 31	mRNA_gene_P21
		0.0	;	# 32	mRNA_gene_P47Phox
		0.0	;	# 33	mRNA_gene_PPARg
		0.0	;	# 34	mRNA_gene_PU1
		0.0	;	# 35	mRNA_gene_Trigger
		0.0	;	# 36	mRNA_gene_cRAF
		0.0	;	# 37	protein_gene_AP1
		0.0	;	# 38	protein_gene_AhR
		0.0	;	# 39	protein_gene_CD11b
		0.0	;	# 40	protein_gene_CD14
		0.0	;	# 41	protein_gene_CD38
		0.0	;	# 42	protein_gene_CEBPa
		0.0	;	# 43	protein_gene_E2F
		0.0	;	# 44	protein_gene_EGR1
		0.0	;	# 45	protein_gene_GFI1
		0.0	;	# 46	protein_gene_IRF1
		0.0	;	# 47	protein_gene_OCT1
		0.0	;	# 48	protein_gene_OCT4
		0.0	;	# 49	protein_gene_P21
		0.0	;	# 50	protein_gene_P47Phox
		0.0	;	# 51	protein_gene_PPARg
		0.0	;	# 52	protein_gene_PU1
		0.0	;	# 53	protein_gene_Trigger
		0.0	;	# 54	protein_gene_cRAF
	]

	binding_parameter_dictionary = Dict{AbstractString,Float64}()
	binding_parameter_dictionary["n_gene_AP1_gene_AhR"] = 1.0
	binding_parameter_dictionary["K_gene_AP1_gene_AhR"] = 120.0
	binding_parameter_dictionary["n_gene_AP1_gene_PU1"] = 1.0
	binding_parameter_dictionary["K_gene_AP1_gene_PU1"] = 120.0
	binding_parameter_dictionary["n_gene_AP1_gene_PPARg"] = 1.0
	binding_parameter_dictionary["K_gene_AP1_gene_PPARg"] = 120.0
	binding_parameter_dictionary["n_gene_AhR_gene_Trigger"] = 1.0
	binding_parameter_dictionary["K_gene_AhR_gene_Trigger"] = 120.0
	binding_parameter_dictionary["n_gene_CD11b_gene_PU1_gene_cRAF"] = 1.0
	binding_parameter_dictionary["K_gene_CD11b_gene_PU1_gene_cRAF"] = 120.0
	binding_parameter_dictionary["n_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"] = 1.0
	binding_parameter_dictionary["K_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"] = 120.0
	binding_parameter_dictionary["n_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"] = 1.0
	binding_parameter_dictionary["K_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"] = 120.0
	binding_parameter_dictionary["n_gene_CEBPa_gene_Trigger"] = 1.0
	binding_parameter_dictionary["K_gene_CEBPa_gene_Trigger"] = 120.0
	binding_parameter_dictionary["n_gene_CEBPa_gene_PPARg"] = 1.0
	binding_parameter_dictionary["K_gene_CEBPa_gene_PPARg"] = 120.0
	binding_parameter_dictionary["n_gene_CEBPa_gene_CEBPa"] = 1.0
	binding_parameter_dictionary["K_gene_CEBPa_gene_CEBPa"] = 120.0
	binding_parameter_dictionary["n_gene_CEBPa_gene_GFI1"] = 1.0
	binding_parameter_dictionary["K_gene_CEBPa_gene_GFI1"] = 120.0
	binding_parameter_dictionary["n_gene_E2F_gene_E2F"] = 1.0
	binding_parameter_dictionary["K_gene_E2F_gene_E2F"] = 120.0
	binding_parameter_dictionary["n_gene_E2F_gene_PPARg"] = 1.0
	binding_parameter_dictionary["K_gene_E2F_gene_PPARg"] = 120.0
	binding_parameter_dictionary["n_gene_E2F_gene_CEBPa"] = 1.0
	binding_parameter_dictionary["K_gene_E2F_gene_CEBPa"] = 120.0
	binding_parameter_dictionary["n_gene_E2F_gene_GFI1"] = 1.0
	binding_parameter_dictionary["K_gene_E2F_gene_GFI1"] = 120.0
	binding_parameter_dictionary["n_gene_E2F_gene_cRAF"] = 1.0
	binding_parameter_dictionary["K_gene_E2F_gene_cRAF"] = 120.0
	binding_parameter_dictionary["n_gene_EGR1_gene_Trigger"] = 1.0
	binding_parameter_dictionary["K_gene_EGR1_gene_Trigger"] = 120.0
	binding_parameter_dictionary["n_gene_EGR1_gene_PU1"] = 1.0
	binding_parameter_dictionary["K_gene_EGR1_gene_PU1"] = 120.0
	binding_parameter_dictionary["n_gene_EGR1_gene_PPARg"] = 1.0
	binding_parameter_dictionary["K_gene_EGR1_gene_PPARg"] = 120.0
	binding_parameter_dictionary["n_gene_EGR1_gene_GFI1"] = 1.0
	binding_parameter_dictionary["K_gene_EGR1_gene_GFI1"] = 120.0
	binding_parameter_dictionary["n_gene_GFI1_gene_CEBPa"] = 1.0
	binding_parameter_dictionary["K_gene_GFI1_gene_CEBPa"] = 120.0
	binding_parameter_dictionary["n_gene_GFI1_gene_EGR1"] = 1.0
	binding_parameter_dictionary["K_gene_GFI1_gene_EGR1"] = 120.0
	binding_parameter_dictionary["n_gene_IRF1_gene_Trigger"] = 1.0
	binding_parameter_dictionary["K_gene_IRF1_gene_Trigger"] = 120.0
	binding_parameter_dictionary["n_gene_IRF1_gene_AhR"] = 1.0
	binding_parameter_dictionary["K_gene_IRF1_gene_AhR"] = 120.0
	binding_parameter_dictionary["n_gene_IRF1_gene_PPARg"] = 1.0
	binding_parameter_dictionary["K_gene_IRF1_gene_PPARg"] = 120.0
	binding_parameter_dictionary["n_gene_OCT1_gene_PPARg"] = 1.0
	binding_parameter_dictionary["K_gene_OCT1_gene_PPARg"] = 120.0
	binding_parameter_dictionary["n_gene_OCT4_gene_Trigger"] = 1.0
	binding_parameter_dictionary["K_gene_OCT4_gene_Trigger"] = 120.0
	binding_parameter_dictionary["n_gene_OCT4_gene_AhR"] = 1.0
	binding_parameter_dictionary["K_gene_OCT4_gene_AhR"] = 120.0
	binding_parameter_dictionary["n_gene_OCT4_gene_cRAF"] = 1.0
	binding_parameter_dictionary["K_gene_OCT4_gene_cRAF"] = 120.0
	binding_parameter_dictionary["n_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"] = 1.0
	binding_parameter_dictionary["K_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"] = 120.0
	binding_parameter_dictionary["n_gene_P21_gene_GFI1"] = 1.0
	binding_parameter_dictionary["K_gene_P21_gene_GFI1"] = 120.0
	binding_parameter_dictionary["n_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"] = 1.0
	binding_parameter_dictionary["K_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"] = 120.0
	binding_parameter_dictionary["n_gene_P47Phox_gene_PPARg"] = 1.0
	binding_parameter_dictionary["K_gene_P47Phox_gene_PPARg"] = 120.0
	binding_parameter_dictionary["n_gene_PPARg_gene_Trigger"] = 1.0
	binding_parameter_dictionary["K_gene_PPARg_gene_Trigger"] = 120.0
	binding_parameter_dictionary["n_gene_PPARg_gene_CEBPa"] = 1.0
	binding_parameter_dictionary["K_gene_PPARg_gene_CEBPa"] = 120.0
	binding_parameter_dictionary["n_gene_PPARg_gene_EGR1"] = 1.0
	binding_parameter_dictionary["K_gene_PPARg_gene_EGR1"] = 120.0
	binding_parameter_dictionary["n_gene_PPARg_gene_PU1"] = 1.0
	binding_parameter_dictionary["K_gene_PPARg_gene_PU1"] = 120.0
	binding_parameter_dictionary["n_gene_PPARg_gene_AP1"] = 1.0
	binding_parameter_dictionary["K_gene_PPARg_gene_AP1"] = 120.0
	binding_parameter_dictionary["n_gene_PU1_gene_Trigger"] = 1.0
	binding_parameter_dictionary["K_gene_PU1_gene_Trigger"] = 120.0
	binding_parameter_dictionary["n_gene_PU1_gene_CEBPa"] = 1.0
	binding_parameter_dictionary["K_gene_PU1_gene_CEBPa"] = 120.0
	binding_parameter_dictionary["n_gene_PU1_gene_PU1"] = 1.0
	binding_parameter_dictionary["K_gene_PU1_gene_PU1"] = 120.0
	binding_parameter_dictionary["n_gene_PU1_gene_AP1"] = 1.0
	binding_parameter_dictionary["K_gene_PU1_gene_AP1"] = 120.0
	binding_parameter_dictionary["n_gene_PU1_gene_OCT1"] = 1.0
	binding_parameter_dictionary["K_gene_PU1_gene_OCT1"] = 120.0
	binding_parameter_dictionary["n_gene_PU1_gene_AhR"] = 1.0
	binding_parameter_dictionary["K_gene_PU1_gene_AhR"] = 120.0
	binding_parameter_dictionary["n_gene_PU1_gene_GFI1"] = 1.0
	binding_parameter_dictionary["K_gene_PU1_gene_GFI1"] = 120.0

	# Alias the control function parameters -
	control_parameter_dictionary = Dict{AbstractString,Float64}()
	control_parameter_dictionary["W_gene_AP1_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_AP1_gene_AhR"] = 1.0
	control_parameter_dictionary["W_gene_AP1_gene_PU1"] = 1.0
	control_parameter_dictionary["W_gene_AP1_gene_PPARg"] = 1.0
	control_parameter_dictionary["W_gene_AhR_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_AhR_gene_Trigger"] = 1.0
	control_parameter_dictionary["W_gene_CD11b_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_CD11b_gene_PU1_gene_cRAF"] = 1.0
	control_parameter_dictionary["W_gene_CD14_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"] = 0.0
	control_parameter_dictionary["W_gene_CD38_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"] = 1.0
	control_parameter_dictionary["W_gene_CEBPa_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_CEBPa_gene_Trigger"] = 1.0
	control_parameter_dictionary["W_gene_CEBPa_gene_PPARg"] = 1.0
	control_parameter_dictionary["W_gene_CEBPa_gene_CEBPa"] = 1.0
	control_parameter_dictionary["W_gene_CEBPa_gene_GFI1"] = 1.0
	control_parameter_dictionary["W_gene_E2F_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_E2F_gene_E2F"] = 1.0
	control_parameter_dictionary["W_gene_E2F_gene_PPARg"] = 1.0
	control_parameter_dictionary["W_gene_E2F_gene_CEBPa"] = 1.0
	control_parameter_dictionary["W_gene_E2F_gene_GFI1"] = 1.0
	control_parameter_dictionary["W_gene_E2F_gene_cRAF"] = 1.0
	control_parameter_dictionary["W_gene_EGR1_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_EGR1_gene_Trigger"] = 1.0
	control_parameter_dictionary["W_gene_EGR1_gene_PU1"] = 1.0
	control_parameter_dictionary["W_gene_EGR1_gene_PPARg"] = 1.0
	control_parameter_dictionary["W_gene_EGR1_gene_GFI1"] = 1.0
	control_parameter_dictionary["W_gene_GFI1_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_GFI1_gene_CEBPa"] = 1.0
	control_parameter_dictionary["W_gene_GFI1_gene_EGR1"] = 1.0
	control_parameter_dictionary["W_gene_IRF1_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_IRF1_gene_Trigger"] = 1.0
	control_parameter_dictionary["W_gene_IRF1_gene_AhR"] = 1.0
	control_parameter_dictionary["W_gene_IRF1_gene_PPARg"] = 1.0
	control_parameter_dictionary["W_gene_OCT1_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_OCT1_gene_PPARg"] = 1.0
	control_parameter_dictionary["W_gene_OCT4_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_OCT4_gene_Trigger"] = 1.0
	control_parameter_dictionary["W_gene_OCT4_gene_AhR"] = 1.0
	control_parameter_dictionary["W_gene_OCT4_gene_cRAF"] = 1.0
	control_parameter_dictionary["W_gene_P21_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"] = 1.0
	control_parameter_dictionary["W_gene_P21_gene_GFI1"] = 1.0
	control_parameter_dictionary["W_gene_P47Phox_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"] = 1.0
	control_parameter_dictionary["W_gene_P47Phox_gene_PPARg"] = 1.0
	control_parameter_dictionary["W_gene_PPARg_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_PPARg_gene_Trigger"] = 1.0
	control_parameter_dictionary["W_gene_PPARg_gene_CEBPa"] = 1.0
	control_parameter_dictionary["W_gene_PPARg_gene_EGR1"] = 1.0
	control_parameter_dictionary["W_gene_PPARg_gene_PU1"] = 1.0
	control_parameter_dictionary["W_gene_PPARg_gene_AP1"] = 1.0
	control_parameter_dictionary["W_gene_PU1_RNAP"] = 0.001
	control_parameter_dictionary["W_gene_PU1_gene_Trigger"] = 1.0
	control_parameter_dictionary["W_gene_PU1_gene_CEBPa"] = 1.0
	control_parameter_dictionary["W_gene_PU1_gene_PU1"] = 1.0
	control_parameter_dictionary["W_gene_PU1_gene_AP1"] = 1.0
	control_parameter_dictionary["W_gene_PU1_gene_OCT1"] = 1.0
	control_parameter_dictionary["W_gene_PU1_gene_AhR"] = 1.0
	control_parameter_dictionary["W_gene_PU1_gene_GFI1"] = 1.0

	# These need to be zero => no background expression -
	control_parameter_dictionary["W_gene_Trigger_RNAP"] = 0.0;
	control_parameter_dictionary["W_gene_cRAF_RNAP"] = 0.0;

	# Parameter name index array -
	parameter_name_mapping_array = [
		"n_gene_AP1_gene_AhR"	;	# 1
		"K_gene_AP1_gene_AhR"	;	# 2
		"n_gene_AP1_gene_PU1"	;	# 3
		"K_gene_AP1_gene_PU1"	;	# 4
		"n_gene_AP1_gene_PPARg"	;	# 5
		"K_gene_AP1_gene_PPARg"	;	# 6
		"n_gene_AhR_gene_Trigger"	;	# 7
		"K_gene_AhR_gene_Trigger"	;	# 8
		"n_gene_CD11b_gene_PU1_gene_cRAF"	;	# 9
		"K_gene_CD11b_gene_PU1_gene_cRAF"	;	# 10
		"n_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 11
		"K_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 12
		"n_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 13
		"K_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 14
		"n_gene_CEBPa_gene_Trigger"	;	# 15
		"K_gene_CEBPa_gene_Trigger"	;	# 16
		"n_gene_CEBPa_gene_PPARg"	;	# 17
		"K_gene_CEBPa_gene_PPARg"	;	# 18
		"n_gene_CEBPa_gene_CEBPa"	;	# 19
		"K_gene_CEBPa_gene_CEBPa"	;	# 20
		"n_gene_CEBPa_gene_GFI1"	;	# 21
		"K_gene_CEBPa_gene_GFI1"	;	# 22
		"n_gene_E2F_gene_E2F"	;	# 23
		"K_gene_E2F_gene_E2F"	;	# 24
		"n_gene_E2F_gene_PPARg"	;	# 25
		"K_gene_E2F_gene_PPARg"	;	# 26
		"n_gene_E2F_gene_CEBPa"	;	# 27
		"K_gene_E2F_gene_CEBPa"	;	# 28
		"n_gene_E2F_gene_GFI1"	;	# 29
		"K_gene_E2F_gene_GFI1"	;	# 30
		"n_gene_E2F_gene_cRAF"	;	# 31
		"K_gene_E2F_gene_cRAF"	;	# 32
		"n_gene_EGR1_gene_Trigger"	;	# 33
		"K_gene_EGR1_gene_Trigger"	;	# 34
		"n_gene_EGR1_gene_PU1"	;	# 35
		"K_gene_EGR1_gene_PU1"	;	# 36
		"n_gene_EGR1_gene_PPARg"	;	# 37
		"K_gene_EGR1_gene_PPARg"	;	# 38
		"n_gene_EGR1_gene_GFI1"	;	# 39
		"K_gene_EGR1_gene_GFI1"	;	# 40
		"n_gene_GFI1_gene_CEBPa"	;	# 41
		"K_gene_GFI1_gene_CEBPa"	;	# 42
		"n_gene_GFI1_gene_EGR1"	;	# 43
		"K_gene_GFI1_gene_EGR1"	;	# 44
		"n_gene_IRF1_gene_Trigger"	;	# 45
		"K_gene_IRF1_gene_Trigger"	;	# 46
		"n_gene_IRF1_gene_AhR"	;	# 47
		"K_gene_IRF1_gene_AhR"	;	# 48
		"n_gene_IRF1_gene_PPARg"	;	# 49
		"K_gene_IRF1_gene_PPARg"	;	# 50
		"n_gene_OCT1_gene_PPARg"	;	# 51
		"K_gene_OCT1_gene_PPARg"	;	# 52
		"n_gene_OCT4_gene_Trigger"	;	# 53
		"K_gene_OCT4_gene_Trigger"	;	# 54
		"n_gene_OCT4_gene_AhR"	;	# 55
		"K_gene_OCT4_gene_AhR"	;	# 56
		"n_gene_OCT4_gene_cRAF"	;	# 57
		"K_gene_OCT4_gene_cRAF"	;	# 58
		"n_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 59
		"K_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 60
		"n_gene_P21_gene_GFI1"	;	# 61
		"K_gene_P21_gene_GFI1"	;	# 62
		"n_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 63
		"K_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 64
		"n_gene_P47Phox_gene_PPARg"	;	# 65
		"K_gene_P47Phox_gene_PPARg"	;	# 66
		"n_gene_PPARg_gene_Trigger"	;	# 67
		"K_gene_PPARg_gene_Trigger"	;	# 68
		"n_gene_PPARg_gene_CEBPa"	;	# 69
		"K_gene_PPARg_gene_CEBPa"	;	# 70
		"n_gene_PPARg_gene_EGR1"	;	# 71
		"K_gene_PPARg_gene_EGR1"	;	# 72
		"n_gene_PPARg_gene_PU1"	;	# 73
		"K_gene_PPARg_gene_PU1"	;	# 74
		"n_gene_PPARg_gene_AP1"	;	# 75
		"K_gene_PPARg_gene_AP1"	;	# 76
		"n_gene_PU1_gene_Trigger"	;	# 77
		"K_gene_PU1_gene_Trigger"	;	# 78
		"n_gene_PU1_gene_CEBPa"	;	# 79
		"K_gene_PU1_gene_CEBPa"	;	# 80
		"n_gene_PU1_gene_PU1"	;	# 81
		"K_gene_PU1_gene_PU1"	;	# 82
		"n_gene_PU1_gene_AP1"	;	# 83
		"K_gene_PU1_gene_AP1"	;	# 84
		"n_gene_PU1_gene_OCT1"	;	# 85
		"K_gene_PU1_gene_OCT1"	;	# 86
		"n_gene_PU1_gene_AhR"	;	# 87
		"K_gene_PU1_gene_AhR"	;	# 88
		"n_gene_PU1_gene_GFI1"	;	# 89
		"K_gene_PU1_gene_GFI1"	;	# 90
		"W_gene_AP1_RNAP"	;	# 91
		"W_gene_AP1_gene_AhR"	;	# 92
		"W_gene_AP1_gene_PU1"	;	# 93
		"W_gene_AP1_gene_PPARg"	;	# 94
		"W_gene_AhR_RNAP"	;	# 95
		"W_gene_AhR_gene_Trigger"	;	# 96
		"W_gene_CD11b_RNAP"	;	# 97
		"W_gene_CD11b_gene_PU1_gene_cRAF"	;	# 98
		"W_gene_CD14_RNAP"	;	# 99
		"W_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 100
		"W_gene_CD38_RNAP"	;	# 101
		"W_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 102
		"W_gene_CEBPa_RNAP"	;	# 103
		"W_gene_CEBPa_gene_Trigger"	;	# 104
		"W_gene_CEBPa_gene_PPARg"	;	# 105
		"W_gene_CEBPa_gene_CEBPa"	;	# 106
		"W_gene_CEBPa_gene_GFI1"	;	# 107
		"W_gene_E2F_RNAP"	;	# 108
		"W_gene_E2F_gene_E2F"	;	# 109
		"W_gene_E2F_gene_PPARg"	;	# 110
		"W_gene_E2F_gene_CEBPa"	;	# 111
		"W_gene_E2F_gene_GFI1"	;	# 112
		"W_gene_E2F_gene_cRAF"	;	# 113
		"W_gene_EGR1_RNAP"	;	# 114
		"W_gene_EGR1_gene_Trigger"	;	# 115
		"W_gene_EGR1_gene_PU1"	;	# 116
		"W_gene_EGR1_gene_PPARg"	;	# 117
		"W_gene_EGR1_gene_GFI1"	;	# 118
		"W_gene_GFI1_RNAP"	;	# 119
		"W_gene_GFI1_gene_CEBPa"	;	# 120
		"W_gene_GFI1_gene_EGR1"	;	# 121
		"W_gene_IRF1_RNAP"	;	# 122
		"W_gene_IRF1_gene_Trigger"	;	# 123
		"W_gene_IRF1_gene_AhR"	;	# 124
		"W_gene_IRF1_gene_PPARg"	;	# 125
		"W_gene_OCT1_RNAP"	;	# 126
		"W_gene_OCT1_gene_PPARg"	;	# 127
		"W_gene_OCT4_RNAP"	;	# 128
		"W_gene_OCT4_gene_Trigger"	;	# 129
		"W_gene_OCT4_gene_AhR"	;	# 130
		"W_gene_OCT4_gene_cRAF"	;	# 131
		"W_gene_P21_RNAP"	;	# 132
		"W_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 133
		"W_gene_P21_gene_GFI1"	;	# 134
		"W_gene_P47Phox_RNAP"	;	# 135
		"W_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 136
		"W_gene_P47Phox_gene_PPARg"	;	# 137
		"W_gene_PPARg_RNAP"	;	# 138
		"W_gene_PPARg_gene_Trigger"	;	# 139
		"W_gene_PPARg_gene_CEBPa"	;	# 140
		"W_gene_PPARg_gene_EGR1"	;	# 141
		"W_gene_PPARg_gene_PU1"	;	# 142
		"W_gene_PPARg_gene_AP1"	;	# 143
		"W_gene_PU1_RNAP"	;	# 144
		"W_gene_PU1_gene_Trigger"	;	# 145
		"W_gene_PU1_gene_CEBPa"	;	# 146
		"W_gene_PU1_gene_PU1"	;	# 147
		"W_gene_PU1_gene_AP1"	;	# 148
		"W_gene_PU1_gene_OCT1"	;	# 149
		"W_gene_PU1_gene_AhR"	;	# 150
		"W_gene_PU1_gene_GFI1"	;	# 151
		"W_gene_Trigger_RNAP"	;	# 152
		"W_gene_cRAF_RNAP"	;	# 153
	]

	# Background copy number - (added by JV)
	background_copy_number_dictionary = Dict{AbstractString,Float64}()
	background_copy_number_dictionary["protein_gene_AP1"] = 100.0
	background_copy_number_dictionary["protein_gene_AhR"] = 100.0
	background_copy_number_dictionary["protein_gene_CD11b"] = 100.0
	background_copy_number_dictionary["protein_gene_CD14"] = 10.0
	background_copy_number_dictionary["protein_gene_CD38"] = 100.0
	background_copy_number_dictionary["protein_gene_CEBPa"] = 100.0
	background_copy_number_dictionary["protein_gene_E2F"] = 10000.0
	background_copy_number_dictionary["protein_gene_OCT4"] = 10000.0
	background_copy_number_dictionary["protein_gene_EGR1"] = 100.0
	background_copy_number_dictionary["protein_gene_GFI1"] = 100.0
	background_copy_number_dictionary["protein_gene_OCT1"] = 100.0
	background_copy_number_dictionary["protein_gene_P21"] = 1000.0
	background_copy_number_dictionary["protein_gene_P47Phox"] = 100.0
	background_copy_number_dictionary["protein_gene_PPARg"] = 100.0
	background_copy_number_dictionary["protein_gene_PU1"] = 1000.0
	background_copy_number_dictionary["protein_gene_IRF1"] = 100.0

	# =============================== DO NOT EDIT BELOW THIS LINE ============================== #
	data_dictionary = Dict{AbstractString,Any}()
	data_dictionary["initial_condition_array"] = initial_condition_array
	data_dictionary["average_transcript_length"] = average_transcript_length
	data_dictionary["average_protein_length"] = average_protein_length
	data_dictionary["gene_coding_length_array"] = gene_coding_length_array
	data_dictionary["mRNA_coding_length_array"] = mRNA_coding_length_array
	data_dictionary["protein_coding_length_array"] = protein_coding_length_array
	data_dictionary["rnapII_concentration"] = rnapII_concentration  # muM
	data_dictionary["ribosome_concentration"] = ribosome_concentration # muM
	data_dictionary["degradation_constant_mRNA"] = degradation_constant_mRNA  # hr^-1
	data_dictionary["degradation_constant_protein"] = degradation_constant_protein  # hr^-1
	data_dictionary["kcat_transcription"] = kcat_transcription  # hr^-1
	data_dictionary["kcat_translation"] = kcat_translation  # hr^-1
	data_dictionary["maximum_specific_growth_rate"] = maximum_specific_growth_rate  # hr^-1
	data_dictionary["death_rate_constant"] = death_rate_constant
	data_dictionary["avg_gene_concentration"] = avg_gene_concentration
	data_dictionary["saturation_constant_transcription"] = saturation_transcription
	data_dictionary["saturation_constant_translation"] = saturation_translation

	data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
	data_dictionary["dilution_matrix"] = dilution_matrix
	data_dictionary["degradation_matrix"] = degradation_matrix

	data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary
	data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary
	data_dictionary["parameter_name_mapping_array"] = parameter_name_mapping_array

	# Added by JV -
	data_dictionary["background_copy_number_dictionary"] = background_copy_number_dictionary;
	data_dictionary["cell_volume"] = V;
	data_dictionary["av_number"] = av_number;
	data_dictionary["number_of_binding"] = 90;

	# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #
	return data_dictionary
end


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

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

# ╔═╡ f6d07dee-fdb9-11ea-004a-f1283db60877
begin
	using Pkg
	Pkg.add.(["Plots", "StatsBase", "PlutoUI", "DataStructures", "GR"])

	using DataStructures
	using Plots
	using StatsBase
	using PlutoUI
end

# ╔═╡ c3dc6d9c-fdb6-11ea-3b74-e1ecfa6c6f49
md"# Probability via computation"

# ╔═╡ f07d314c-fdb6-11ea-0e6e-173e625133cf
md"""
We would like to model how an epidemic spreads in a population. There are many approaches to do so.

One approach would be to think about individual people, or **agents**. They move around in an environment and interact with other agents when they are close by. If I spend time near a person who is infectious, I might catch the infection. It's too hard to model the details of the physical process by which virus particles get transmitted from one person to another, so we will model that as saying that I have a certain **probability** (chance) $p_I$ to get infected whenever I meet somebody infectious.
"""

# ╔═╡ 8fbf80ac-fe0c-11ea-1848-9fae8515757f
md"![](https://i.imgur.com/OzlPEjn.gif)"

# ╔═╡ 89cbbbe8-fe0c-11ea-1cf8-e33f6bd7cee2
md"""
Once I'm infected, a complicated process takes place in my immune system which is also difficult to model. So we could say that each day I have a certain **probability** $p_R$ to recover.

In order to implement an **agent-based model** (or "individual-based model") like this, we first need to understand what it means for something to happen with probability $p$, and how to implement that on a computer.

We will take a practical and pragmatic approach to probability *via computation*, and leave discussions of the philosophical underpinnings of probability to other venues!
"""

# ╔═╡ e906450c-fdb6-11ea-13c3-d1348fdae587
md"## Random sampling"

# ╔═╡ cd32d0c6-fdb7-11ea-3d3f-d9981951293f
md"We have already generated random objects using the `rand` function. Let's look at that function in a bit more detail. In this notebook we won't discuss how randomness is generated on a computer; we will just assume that these methods are already available.

[Note that computers cannot generate true randomness -- they simulate it using complicated deterministic processes. There are [ways](https://www.random.org/) to obtain true randomness from physical processes, such as atmospheric or electronic noisde.]"

# ╔═╡ ebccc984-fdb9-11ea-2f94-997184ccc66d
md"Calling `rand` on a collection like an array or a tuple returns one of the elements in that collection with equal (**uniform**) probability, i.e. each element of the collection has the same chance of being returned.

For example, we can choose a random friend to call:
"

# ╔═╡ c55235d4-fdbc-11ea-0942-45ac804aa51d
friends = ["Alicia", "Jaime", "Elena", "Roberto"]

# ╔═╡ cf2a7ec2-fdbc-11ea-1010-190de6ff85d1
rand(friends)

# ╔═╡ 1ef8cf70-fe13-11ea-257b-cbc9abd8ab3b
rand( (1, 2) )

# ╔═╡ f317779a-fdbc-11ea-1afa-2dc25c9a16f8
md"Choosing a random element like this is called (uniform) **sampling**.

The `rand` function has many methods that sample from different types of objects."

# ╔═╡ e31461dc-fdbc-11ea-1398-89b7ba513f5c
md"## Rolling a die"

# ╔═╡ b6aec1dc-fdbc-11ea-2345-831cff1d5e7c
 md"""Let's roll an (unbiased) die. ["Die" is the singular of "dice".]"""

# ╔═╡ 2fea1cb6-fdba-11ea-2744-f5dd5fe0fcd7
begin
	num_sides = 6
	sides = 1:num_sides
end

# ╔═╡ 74430946-fe95-11ea-31d6-31bde3a357ba
sides

# ╔═╡ 75e017d0-fe95-11ea-37b1-357616582a9e
typeof(sides)

# ╔═╡ 5a5d74d4-fdba-11ea-16a3-ed2c24f47821
rand(sides)

# ╔═╡ 61f7ab10-fdba-11ea-271c-d17d8d4ebc8e
md"
We can repeat this sampling $N$ times by adding a second argument to `rand`:
"

# ╔═╡ a36f553e-fdba-11ea-2a6a-b3ba05d11acd
md"What does it mean for the sampling to be uniform? It means that we expect each element to occur with the same **frequency** (number of occurrences).

We can compute the frequencies using the `countmap` function from the `StatsBase.jl` package. But you should think about how we could write this function ourselves!"

# ╔═╡ a03b164e-fdba-11ea-0c7f-a59d54b2f813
sample = rand(sides, 10)

# ╔═╡ 10b83e00-fdbe-11ea-260d-65cf5e080534
frequencies = StatsBase.countmap(sample)

# ╔═╡ a6da8b1a-fe89-11ea-364c-958bddfa4612
bar(frequencies, alpha=0.5,
	size=(400, 300), leg=false,
	xlabel="outcome", ylabel="number of rolls", xlim=(0, num_sides+1))

# ╔═╡ 2d0be59e-fdbf-11ea-0b86-a50db0aaccaa
md"We expect the **relative frequency** or **proportion** to be near $1/6$:"

# ╔═╡ db102cb6-fdbf-11ea-0408-61334cb89f6d
begin
	max_rolls = 10000
	rolls2 = rand(1:num_sides, max_rolls)
end

# ╔═╡ b596f406-fdbf-11ea-0ed7-ab15d2da9fd2
@bind num_rolls2 Slider(1:max_rolls, show_value=true)

# ╔═╡ 46ae7d02-fdbf-11ea-1f19-4f94f9bbf851
begin
	freqs2 = SortedDict(StatsBase.countmap(rolls2[1:num_rolls2]))

	ks = collect(keys(freqs2))
	vs = collect(values(freqs2)) ./ num_rolls2

	bar(ks, vs, leg=false, alpha=0.5, xlims=(0, num_sides+1),
		size=(400,300),
		xlabel="value", ylabel="relative frequency")

	hline!([1 / num_sides], ls=:dash, lw=3, c=:red)

	 ylims!(0, 0.3)
end

# ╔═╡ f97c2aa6-fdbf-11ea-11bf-2b44cd0a90bc
md"Note that we **pre-generated** the data, to avoid generating new random samples each time."

# ╔═╡ 24674186-fdc5-11ea-2a67-3bd71199ae6e
md"## Random variables"

# ╔═╡ 2b954fca-fdc5-11ea-149f-e5f7f6d49407
md"Let's call $X$ the outcome of flipping one coin. Each time we run the experiment, $X$ will take a different value. This makes $X$ an example of a **random variable**.

[Giving a proper mathematical description of this is rather tricky; see an advanced probability course.]

We write

$$\mathbb{P}(X=1) = \textstyle \frac{1}{6}$$

to say that the probability that $X$ takes the value 1 is $\frac{1}{6}$.
"

# ╔═╡ 112c47a4-fe0f-11ea-04a4-adc7e339a352
md"## Uniform random numbers"

# ╔═╡ 69c2c1e4-fdc0-11ea-0e76-335b6e06057d
md"From the above we can see how to do something with probability $1/6$:

> roll a 6-sided die, and return `true` if the outcome is 1 (for example).

Let's extend this to events which occur with a given probability $p$ with $0 \le p \le 1$: How could we sample an event with probability $p$?"

# ╔═╡ cbe5bbc6-fe0e-11ea-17c7-095c0d485bdd
md"""Calling the `rand()` function with *no* arguments returns a random `Float64` between $0$ and $1$. This approximates uniform sampling of a *real* number between $0$ and $1$."""

# ╔═╡ 6d7f1ca4-fe8a-11ea-3b20-b3df98a4201b
rand(Float64)

# ╔═╡ d5743c58-fe0e-11ea-0391-c7591b944bdd
rand()

# ╔═╡ bf084b00-fdc1-11ea-04ca-1757996125fd
md"""
Let's see what this looks like. We'll sample uniform random variates in $[0, 1]$ and plot them as a function of "time" on the $x$ axis (i.e. the number of the trial in which they appeared). We also plot the outcomes along the $y$ axis so that we can see how the interval $[0, 1]$ gets covered:
"""

# ╔═╡ 019bd65c-fdc3-11ea-01ba-51d497f875d2
max_samples = 500

# ╔═╡ 2dd45d84-fdc3-11ea-1b61-6bc6467a5013
sample2 = rand(max_samples)

# ╔═╡ 08af602e-fdc3-11ea-0e47-af1047c7bd13
@bind num_samples Slider(1:max_samples, show_value=true)

# ╔═╡ da913e08-fdc2-11ea-14c1-6bb4fdaf9353
begin
	r = sample2[1:num_samples]
	scatter(r, 1:length(r), alpha=0.5, leg=false,
			size=(400, 300))
	scatter!(r, zeros(length(r)), alpha=0.5)

	ylims!(-10, max_samples)
	xlims!(0, 1)

	xlabel!("random variate")
	ylabel!("position in sequence (time)")
end

# ╔═╡ e5994394-fe0e-11ea-1350-51442f863799
md"How uniform is the sampling? Let's plot a **histogram**. This splits the interval $[0, 1]$ up into boxes and counts the number (frequency) of data that fall into each box:"

# ╔═╡ 3620150e-fe0f-11ea-0eff-91a85ccf3864
@bind bins Slider(10:1000, show_value=true)

# ╔═╡ 479e0232-fe0f-11ea-308e-2928fe213807
@bind random_samples Slider(1:max_samples, show_value=true)

# ╔═╡ 022adcee-fde1-11ea-1ae5-8bcf6b129404
histogram(sample2[1:random_samples], bins=bins, alpha=0.5, leg=false,
		size=(400, 300), xlims=(0, 1))

# ╔═╡ 5d4102dc-fdc0-11ea-049f-7bc11f5b728d
md"## Sampling events with a given probability: Bernoulli trials"

# ╔═╡ 86731f7a-fdc3-11ea-2955-6b938f3346a5
md"""
How can we use *uniform* sampling to do *non-uniform* sampling? For example, to simulate an event that occurs with probability 0.25? This means that when we sample many points, the proportion that fall into $X$ should be approximately 0.25.

Since the sampling is uniform, we can take $X = [0, 0.25]$. The definition of "uniform" is that the probability is proportional to the length of the interval.

The recipe for sampling an event with probability $p$ is then as follows:

> 1. Generate a uniform random variate $r \in [0, 1]$, using `rand()`.
>  2. If $r < p$ then return `true`, else return `false`.
"""

# ╔═╡ ba9939f0-fe0f-11ea-0328-6780c29cc01c
md"To get a feeling for why this is a good definition, let's go back to our picture of uniform sampling and add in the cutoff at height $p$:"

# ╔═╡ ca62d4a4-fe0f-11ea-3623-73a11752dc8c
@bind num_samples2 Slider(1:max_samples, show_value=true)

# ╔═╡ 8fa5e992-fde1-11ea-0c65-3d83d42729ea
@bind pp Slider(0:0.01:1, show_value=true)

# ╔═╡ 39b9d2a0-fde1-11ea-2189-97a95849152f
begin
	r2 = sample2[1:num_samples2]

	scatter(r2, 1:length(r2), alpha=0.5, leg=false,
			size=(400, 300))
	scatter!(r2, zeros(length(r2)), alpha=0.5)

	ylims!(-10, max_samples)
	xlims!(0, 1)

	vline!([pp], ls=:dash, lw=3, c=:green)
	annotate!(170, pp+0.06, text("p", color=:green))

	which = findall(x -> x .< pp, r2)
	scatter!(r2[which], which, m=:square, alpha=0.5)
	scatter!(r2[which], zeros(length(r2[which])), m=:square, alpha=0.5)

end

# ╔═╡ f80a700a-fe15-11ea-2c48-9f499842e45d
md"We expect that out of $N$ trials, about $pN$ will be below $p$. So the proportion of the total number is $(pN) / N \simeq p$.
"

# ╔═╡ 1fa92ed0-fdc4-11ea-0f6c-9178b21b259c
md"The above procedure is called a [**Bernoulli trial**](https://en.wikipedia.org/wiki/Bernoulli_trial) with probability $p$.
We can think of it as flipping a coin that is **biased** to fall on heads a higher proportion of the time (e.g. by making the coin out of two different metals, one of which is more dense than the other).

Julia allows us to code this in the following very compact way:"

# ╔═╡ 4bf82e84-fdc4-11ea-2abf-19174c74f726
bernoulli(p) = rand() < p

# ╔═╡ bc549f38-fe15-11ea-10dd-db165347709d
3 < 4

# ╔═╡ 9ab21d04-fe97-11ea-0f88-11a86a4974c6
function bernoulli2(p)
	if rand() < p
		return true
	else
		return false
	end
end

# ╔═╡ 4f9ebeea-fdc4-11ea-052d-3b9e7dfbc02b
bernoulli(0.25)

# ╔═╡ 5640f132-fdc4-11ea-372e-799a17716334
md"We are using the short-form version of a function definition, together with the fact that a comparison like `x < y` returns a Boolean value."

# ╔═╡ db63e988-fde1-11ea-152e-8b194c107e78
md"Note that the value returned, $B$, is a random variable! It has

$$\mathbb{P}(B = 1) = p$$
$$\mathbb{P}(B = 0) = 1 - p$$
"

# ╔═╡ 11e7a0b2-fde2-11ea-0c37-6be01573167b
md"This collection of the probabilities of all the possible outcomes is the **probability distribution** of the Bernoulli random variable."

# ╔═╡ c9b287a2-fdc4-11ea-04c6-6b0d1a6c5321
md"## Flipping many biased coins"

# ╔═╡ ca7b9150-fdc5-11ea-2281-b3b026825ea5
md"
Let's suppose we flip $n$ coins, each with the same probability $p$. We will suppose that the coins are **independent**, i.e. that flipping one of them has no influence on the outcome of the others. (This would not be the case if the coins were metallic and attracted each other, for example.)

There are many quantities of interest, or **observables**, that we might want to measure. One obvious one is the number of coins that come out as heads:
"

# ╔═╡ 60949cce-fde2-11ea-1946-31d1b1bce37a
[bernoulli(0.3) for i = 1:10]  # array comprehension

# ╔═╡ 727b0e56-fdd7-11ea-28d1-5b859a5f8a29
flips(n, p) = count( bernoulli(p) for i = 1:n )

# ╔═╡ dabc00d6-fe97-11ea-0180-cd138faaab5d
( bernoulli(0.3) for i = 1:10 )

# ╔═╡ d450742c-fdd7-11ea-2477-b37441b6ddfa
flips(20, 0.3)

# ╔═╡ dca0b5a6-fdd7-11ea-15e4-0320a1ceebff
md"The function `count` counts the number of `true`s in its argument. The argument itself looks like an array comprehension, but with no square brackets; this is called a **generator expression**. Since `count` does not need to **materialize** the array, i.e. create it in memory, it can be faster."

# ╔═╡ 2c2479a0-fdd8-11ea-19a2-db21219dd7c3
md"""The outcome of flipping the coin is another random variable, $H_{n}$. Clearly its possible values are the integers between $0$ and $n$. But intuitively we expect that it's very *unlikely* to have either exactly $0$ or exactly $n$ heads. This is another example of a *non-uniform* random variable.

What again need to calculate the **probability distribution** of the random variable $H_{n,p}$, i.e. the probability that it takes on each of its possible values.

We can calculate this numerically by thinking of "flipping $n$ coins" as a single experiment, and running that experiment many times:
"""

# ╔═╡ aa63b25e-fdd8-11ea-2b30-ed13d2072488
run_experiment(n, p, num_times) = [flips(n, p) for i in 1:num_times]

# ╔═╡ 48acfd94-fde3-11ea-3ae8-0d8ba164a68c
num_experiments = 1000

# ╔═╡ 521c85c0-fde3-11ea-1040-09a6e13a4328
num_coins = 100

# ╔═╡ c55013e6-fdd8-11ea-3784-cd27c2e3d0d8
run_experiment(num_coins, 0.3, num_experiments)

# ╔═╡ caf791ca-fdd8-11ea-2234-bdc14d214cac
md"As expected, we get back a dataset containing integers between 0 and 20, but actually 1 and 20 occur few times.

How can we calculate the probability distribution? We can use `countmap` again to obtain the frequencies:"

# ╔═╡ 15231f76-fdd9-11ea-361e-0f85e1c3ace2
begin
	max_times = 10^6
	data = run_experiment(num_coins, 0.3, max_times)
end

# ╔═╡ 29edd4e6-fdd9-11ea-016b-9142be69bf68
@bind num_times Slider(1:max_times, show_value=true)

# ╔═╡ 3398d3a6-fdd9-11ea-0845-f54dab63b2dc
begin
	freqs3 = SortedDict(countmap(data[1:num_times]))
	bar(collect(keys(freqs3)), collect(values(freqs3)) ./ num_times,
		alpha=0.5,
		size=(400, 300), leg=false)

	xlims!(0, num_coins)
	ylims!(0, 0.2)

	xlabel!("number of heads")
	ylabel!("proportion")
end

# ╔═╡ adf1c626-fdd9-11ea-3535-2d2219feda56
md"As we take more samples, the **empirical distribution** (calculated from the data) converges to the true underlying **population distribution** (the theoretical distribution for this random variable)."

# ╔═╡ 92efff82-fde3-11ea-010d-a188a852b79b
md"If we take a large number of coins then we see the well-known bell-shaped distribution, called a  **Gaussian distribution** or **normal distribution**.

If we call the $i$th coin flip result $B_i$, which takes the value 0 or 1, then we have

$$H_{n} = \sum_{i=1}^n B_i.$$

We are summing up independent random variables which take the value 0 or 1 (the number of heads in each coin flip), and the Central Limit Theorem says that the sum of independent random variables usually converges to a normal distribution.

Note that we have not said what we mean by **converges** here. Making that precise turns out to be quite difficult. Intuitively it means that the function describing the probability distribution converges at each point to a limiting curve, but unfortunately this intuition is not always correct.
"

# ╔═╡ Cell order:
# ╠═f6d07dee-fdb9-11ea-004a-f1283db60877
# ╟─c3dc6d9c-fdb6-11ea-3b74-e1ecfa6c6f49
# ╟─f07d314c-fdb6-11ea-0e6e-173e625133cf
# ╟─8fbf80ac-fe0c-11ea-1848-9fae8515757f
# ╟─89cbbbe8-fe0c-11ea-1cf8-e33f6bd7cee2
# ╟─e906450c-fdb6-11ea-13c3-d1348fdae587
# ╟─cd32d0c6-fdb7-11ea-3d3f-d9981951293f
# ╟─ebccc984-fdb9-11ea-2f94-997184ccc66d
# ╠═c55235d4-fdbc-11ea-0942-45ac804aa51d
# ╠═cf2a7ec2-fdbc-11ea-1010-190de6ff85d1
# ╠═1ef8cf70-fe13-11ea-257b-cbc9abd8ab3b
# ╟─f317779a-fdbc-11ea-1afa-2dc25c9a16f8
# ╟─e31461dc-fdbc-11ea-1398-89b7ba513f5c
# ╟─b6aec1dc-fdbc-11ea-2345-831cff1d5e7c
# ╠═2fea1cb6-fdba-11ea-2744-f5dd5fe0fcd7
# ╠═74430946-fe95-11ea-31d6-31bde3a357ba
# ╠═75e017d0-fe95-11ea-37b1-357616582a9e
# ╠═5a5d74d4-fdba-11ea-16a3-ed2c24f47821
# ╟─61f7ab10-fdba-11ea-271c-d17d8d4ebc8e
# ╟─a36f553e-fdba-11ea-2a6a-b3ba05d11acd
# ╠═a03b164e-fdba-11ea-0c7f-a59d54b2f813
# ╠═10b83e00-fdbe-11ea-260d-65cf5e080534
# ╟─a6da8b1a-fe89-11ea-364c-958bddfa4612
# ╟─2d0be59e-fdbf-11ea-0b86-a50db0aaccaa
# ╠═db102cb6-fdbf-11ea-0408-61334cb89f6d
# ╠═b596f406-fdbf-11ea-0ed7-ab15d2da9fd2
# ╟─46ae7d02-fdbf-11ea-1f19-4f94f9bbf851
# ╟─f97c2aa6-fdbf-11ea-11bf-2b44cd0a90bc
# ╟─24674186-fdc5-11ea-2a67-3bd71199ae6e
# ╟─2b954fca-fdc5-11ea-149f-e5f7f6d49407
# ╟─112c47a4-fe0f-11ea-04a4-adc7e339a352
# ╟─69c2c1e4-fdc0-11ea-0e76-335b6e06057d
# ╟─cbe5bbc6-fe0e-11ea-17c7-095c0d485bdd
# ╠═6d7f1ca4-fe8a-11ea-3b20-b3df98a4201b
# ╠═d5743c58-fe0e-11ea-0391-c7591b944bdd
# ╟─bf084b00-fdc1-11ea-04ca-1757996125fd
# ╠═019bd65c-fdc3-11ea-01ba-51d497f875d2
# ╠═2dd45d84-fdc3-11ea-1b61-6bc6467a5013
# ╟─08af602e-fdc3-11ea-0e47-af1047c7bd13
# ╟─da913e08-fdc2-11ea-14c1-6bb4fdaf9353
# ╟─e5994394-fe0e-11ea-1350-51442f863799
# ╠═3620150e-fe0f-11ea-0eff-91a85ccf3864
# ╠═479e0232-fe0f-11ea-308e-2928fe213807
# ╠═022adcee-fde1-11ea-1ae5-8bcf6b129404
# ╟─5d4102dc-fdc0-11ea-049f-7bc11f5b728d
# ╟─86731f7a-fdc3-11ea-2955-6b938f3346a5
# ╟─ba9939f0-fe0f-11ea-0328-6780c29cc01c
# ╠═ca62d4a4-fe0f-11ea-3623-73a11752dc8c
# ╠═8fa5e992-fde1-11ea-0c65-3d83d42729ea
# ╟─39b9d2a0-fde1-11ea-2189-97a95849152f
# ╟─f80a700a-fe15-11ea-2c48-9f499842e45d
# ╟─1fa92ed0-fdc4-11ea-0f6c-9178b21b259c
# ╠═4bf82e84-fdc4-11ea-2abf-19174c74f726
# ╠═bc549f38-fe15-11ea-10dd-db165347709d
# ╠═9ab21d04-fe97-11ea-0f88-11a86a4974c6
# ╠═4f9ebeea-fdc4-11ea-052d-3b9e7dfbc02b
# ╟─5640f132-fdc4-11ea-372e-799a17716334
# ╟─db63e988-fde1-11ea-152e-8b194c107e78
# ╟─11e7a0b2-fde2-11ea-0c37-6be01573167b
# ╟─c9b287a2-fdc4-11ea-04c6-6b0d1a6c5321
# ╟─ca7b9150-fdc5-11ea-2281-b3b026825ea5
# ╠═60949cce-fde2-11ea-1946-31d1b1bce37a
# ╠═727b0e56-fdd7-11ea-28d1-5b859a5f8a29
# ╠═dabc00d6-fe97-11ea-0180-cd138faaab5d
# ╠═d450742c-fdd7-11ea-2477-b37441b6ddfa
# ╟─dca0b5a6-fdd7-11ea-15e4-0320a1ceebff
# ╟─2c2479a0-fdd8-11ea-19a2-db21219dd7c3
# ╠═aa63b25e-fdd8-11ea-2b30-ed13d2072488
# ╠═48acfd94-fde3-11ea-3ae8-0d8ba164a68c
# ╠═521c85c0-fde3-11ea-1040-09a6e13a4328
# ╠═c55013e6-fdd8-11ea-3784-cd27c2e3d0d8
# ╟─caf791ca-fdd8-11ea-2234-bdc14d214cac
# ╠═15231f76-fdd9-11ea-361e-0f85e1c3ace2
# ╠═29edd4e6-fdd9-11ea-016b-9142be69bf68
# ╟─3398d3a6-fdd9-11ea-0845-f54dab63b2dc
# ╟─adf1c626-fdd9-11ea-3535-2d2219feda56
# ╟─92efff82-fde3-11ea-010d-a188a852b79b

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
for i = 1:1000
    for j = 1:(1000-i)
        k = (1000-i-j)
        if i^2+j^2 == k^2
            println(i*j*k)
            quit()
        end
    end
end

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving