text:
using REPLCompletions
using Base.Test

module CompletionFoo
    module CompletionFoo2

    end
    const bar = 1
    foo() = bar
end

test_complete(s) = completions(s,endof(s))
test_scomplete(s) = shell_completions(s,endof(s))

s = ""
c,r = test_complete(s)
@test in("CompletionFoo",c)
@test r == 0:-1
@test s[r] == ""

s = "Comp"
c,r = test_complete(s)
@test in("CompletionFoo",c)
@test r == 1:4
@test s[r] == "Comp"

s = "Main.Comp"
c,r = test_complete(s)
@test in("CompletionFoo",c)
@test r == 6:9
@test s[r] == "Comp"

s = "Main.CompletionFoo."
c,r = test_complete(s)
@test in("bar",c)
@test r == 20:19
@test s[r] == ""

s = "Main.CompletionFoo.f"
c,r = test_complete(s)
@test in("foo",c)
@test r == 20:20
@test s[r] == "f"

# Test completion of packages
mkp(p) = ((@assert !isdir(p)); mkdir(p))
mkp(Pkg.dir("MyAwesomePackage"))
mkp(Pkg.dir("CompletionFooPackage"))

s = "using MyAwesome"
c,r = test_complete(s)
@test in("MyAwesomePackage",c)
@test s[r] == "MyAwesome"

s = "using Completion"
c,r = test_complete(s)
@test in("CompletionFoo",c) #The module
@test in("CompletionFooPackage",c) #The package
@test s[r] == "Completion"

s = "using CompletionFoo.Completion"
c,r = test_complete(s)
@test in("CompletionFoo2",c)
@test s[r] == "Completion"

rmdir(Pkg.dir("MyAwesomePackage"))
rmdir(Pkg.dir("CompletionFooPackage"))

@unix_only begin
    #Assume that we can rely on the existence and accessibility of /tmp
    s = "/t"
    c,r = test_scomplete("/t")
    @test in("tmp/",c)
    @test r == 2:2
    @test s[r] == "t"

    s = "/tmp"
    c,r = test_scomplete(s)
    @test in("tmp/",c)
    @test r == 2:4
    @test s[r] == "tmp"

    # This should match things that are inside the tmp directory
    if !isdir("/tmp/tmp")
        s = "/tmp/"
        c,r = test_scomplete(s)
        @test !in("tmp/",c)
        @test r == 6:5
        @test s[r] == ""
    end

    s = "cd \$(Pk"
    c,r = test_scomplete(s)
    @test in("Pkg",c)
    @test r == 6:7
    @test s[r] == "Pk"
end


text:
# This script should eventually become a ultimate copy of the all encompasing
# textbook on multivariate time-series: Lütkepohl (2007). Some of the methods will be added
# from other papers, such as generalised VAR structure from Pesaran, Shin or computing
# all reorderings in the standard identification structure.

#
###### Definition of the specialised type
#

# The definition of the type rests on several specific neccessities
# There are things that are neccessary for proper definition of the VAR model and as
# such have to be included in the definition of the type. That is for example number
# of lags, the original data, the type.
# Other things can be foregone, but are often desirable to be stored, because they are
# computationally intensive and computing them only once is typically much better than
# recomputing them every time. However, in many applications, one might not even need them.
#
# It is neccessary to balance the needs and strike the balance between initial speed of
# estimation and future speed. For example, I probably wouldn't like to have bootstrapped
# confidence intervals for some things, but I would surely keep the possibility to have them.
#
# Moreover, the current estimation technique will probably be replaced by the method from
# function restrictVAR2, because it is more general.

import Base: show

function estimateVAR(data::Matrix{Float64}, lags::Int, typ::String)
	X = dataMatrix(data, lags, typ)
	Y = data[(1+lags):end,:]
	obs, vars = size(data)
	C = (X'*X) \ X'*Y
	obs = obs - lags
	# Get numerical type
	typ=="None" && (ntyp = 0)
	typ=="Const" && (ntyp = 1)
	typ=="Const and trend" && (ntyp = 2)
	typ=="Trend" && (ntyp = 1)
	u = (Y-X*C)
	Σ = Hermitian(1/(obs-lags*vars-ntyp) * u'*u)

	return (X, Y, obs, vars, C, Σ, ntyp)
end

function estimateVARShrinkage(data::Matrix{Float64}, lags::Int, typ::String, λ::Float64)
	X = dataMatrix(data, lags, typ)
	Y = data[(1+lags):end,:]
	obs, vars = size(data)
	C = (X'*X + λ * obs * eye(size(X)[2])) \ X'*Y
	obs = obs - lags
	# Get numerical type
	typ=="None" && (ntyp = 0)
	typ=="Const" && (ntyp = 1)
	typ=="Const and trend" && (ntyp = 2)
	typ=="Trend" && (ntyp = 1)
	u = (Y-X*C)
	Σ = 1/(obs-lags*vars-ntyp) * u'*u

	return (X, Y, obs, vars, C, Σ, ntyp)
end

type varEstimate <: VARRepresentation
	lags::Int
	typ::String
	ntyp::Int
	vars::Int
	obs::Int
	X::Matrix{Float64}
	Y::Matrix{Float64}
	C::Matrix{Float64}
	Σ::Matrix{Float64}
	HPsi::Int
	HPhi::Int
	Phi::Array{Float64, 3}
	Psi::Array{Float64, 3}
	seriesNames::Vector{String}

	function varEstimate(data::Matrix{Float64},lags::Int, typ::String)
		(X, Y, obs, vars, C, Σ, ntyp) = estimateVAR(data::Matrix{Float64},lags::Int,typ::String)
		new(lags, typ, ntyp, vars, obs, X, Y, C, Σ, 0, 0, zeros(vars,vars,1), zeros(vars,vars,1), ["" for i=1:vars])
	end
end

type varEstimateShrink <: VARRepresentation
	lags::Int
	typ::String
	ntyp::Int
	vars::Int
	obs::Int
	X::Matrix{Float64}
	Y::Matrix{Float64}
	C::Matrix{Float64}
	Σ::Matrix{Float64}
	HPsi::Int
	HPhi::Int
	Phi::Array{Float64, 3}
	Psi::Array{Float64, 3}
	seriesNames::Vector{String}
	λ

	function varEstimateShrink(data::Matrix{Float64},lags::Int, typ::String, λ)
		(X, Y, obs, vars, C, Σ, ntyp) = estimateVARShrinkage(data::Matrix{Float64},lags::Int,typ::String, λ)
		new(lags, typ, ntyp, vars, obs, X, Y, C, Σ, 0, 0, zeros(vars,vars,1), zeros(vars,vars,1), ["" for i=1:vars], λ)
	end
end

function Clag(est::VARRepresentation, lag)
	C = est.C
	ntyp = est.ntyp
	(coefs, vars) = size(C)
	lags = (coefs-ntyp)/vars
	@assert lags >= lag "The demanded lag is higher than existing"
	return C[1+ntyp+(((lag-1)*vars):((lag)*vars-1)),:]
end

function λ(est::VARRepresentation)
	# using Lutkepohl pg.139
	T = est.obs - est.lags
	K = est.vars
	β = vec(est.C)
	y = vec(est.Y)
	Z = est.X'
	Σ = est.Σ

	λv = -(K*T)/2 * log(2π) - T/2 * log(det(Σ)) - 1/2 * (y-kron(eye(K),Z')*β)'*(kron(inv(Σ), eye(T)))*(y-kron(eye(K),Z')*β)
	return λv[1]
end

function getR(a::Matrix{Bool})
	b = prod(size(a))
	c = sum(a)
	mat = fill(0.0, b, c)
	a = vec(a)
	d=1
	for i=1:b
		if a[i]
			mat[i,:] = [fill(0.0, d-1), 1, fill(0.0, c-d)]
			d += 1
		else
			mat[i,:] = fill(0.0, c)
		end
	end

	return mat
end

function restrictVAR(est::varEstimate, restrictions::Expr, egls::Bool=false)
	# TODO: Add small r from Lutkepohl 5.2
	# The original data are not changed, due to very different structure in the estimation (vectorisation)
	# The matrix R as the argument has to be of the size of the coefficients
	# Y stay the same

	# Naming conventions, this does not eat up memory, it is just renaming for more readability,
	# it is parsed out
	Z = est.X'
	Y = est.Y
	Σ = est.Σ
	L = est.lags
	K = est.vars
	k = est.obs-est.lags
	R = restrictions.R
	r = restrictions.r

	z = vec(Y) - kron(eye(K), Z')*r

	if egls
		Σ1 = inv(Σ)
		γ = R'*kron(Σ1, Z*Z')*R \ R'*kron(Σ1, Z)*z
	else
		γ = (R'*kron(eye(K), Z*Z')*R) \ R'*kron(eye(K), Z)*z
	end
	β = R*γ + r
	est.C = reshape(β, L*K+1, K)

	est.Σ = (est.X*est.C-est.Y)'*(est.X*est.C-est.Y) / k
	return est
end

function restrictVAR2(est::varEstimate, R::Matrix{Bool}, egls::Bool=false)
	# TODO: Add small r from Lutkepohl 5.2
	# The original data are not changed, due to very different structure in the estimation (vectorisation)
	# The matrix R as the argument has to be of the size of the coefficients
	# Y stay the same

	# Naming conventions, this does not eat up memory, it is just renaming for more readability,
	# it is parsed out
	Z = est.X'
	Y = est.Y
	Σ = est.Σ
	L = est.lags
	K = est.vars
	k = est.obs-est.lags

	# Temporary restriction to restrict betas only to 0
	r = fill(0.0, prod(size(R)))

	z = vec(Y) - kron(eye(K), Z')*r
	R = getR(R')

	if egls
		Σ1 = inv(Σ)
		γ = R'*kron(Σ1, Z*Z')*R \ R'*kron(Σ1, Z)*z
	else
		γ = (R'*kron(eye(K), Z*Z')*R) \ R'*kron(eye(K), Z)*z
	end
	β = R*γ + r
	est.C = reshape(β, L*K+1, K)

	est.Σ = (est.X*est.C-est.Y)'*(est.X*est.C-est.Y) / k
	return est
end

# Nejde rychlejc, zkusil jsem spoustu možností včetně prealokací
function Phi(estimate::VARRepresentation, H)
	phi = zeros(estimate.vars, estimate.vars, H)
	phi[:,:,1] = eye(estimate.vars)

	for j=2:H
		k = j > estimate.lags ? estimate.lags + 1 : j
		for i=1:(k-1)
			phi[:,:,j] += phi[:,:,j-i]*Clag(estimate, i)
		end
	end
	estimate.Phi = phi
	estimate.HPhi = H
	return phi
end

function Psi(estimate::VARRepresentation, H)
	Phi(estimate, H)
	psi = deepcopy(estimate.Phi)
	P = chol(estimate.Σ)
	for i=1:H
		psi[:,:,i] = P*psi[:,:,i]
	end
	estimate.Psi = psi
	estimate.HPsi = H
	return psi
end

function show(io::IO, a::VARRepresentation)
	line = "--------------------------------------------------------"
	println(io, line)
	println(io, "                    VAR ESTIMATE")
	println(io, line)
	print(io, "Series: ")
	i=1
	while i < (length(a.seriesNames))
		print(io, (a.seriesNames)[i])
		print(io, ", ")
		i += 1
	end
	print(io, a.seriesNames[end])
	print(io, "\n")

	print(io, "Lags: ")
	println(io, a.lags)

	print(io, "Type: ")
	println(io, a.typ)

	println(io, line)

	if (a.typ == "Const") || (a.typ == "Trend")
		println(io, a.typ)
		println(io, mapslices((x)->join(x, "\n"), mapslices((x)->join(x, "  "), map((x)->format(x, precision = 2, signed = true), a.C[1:1,:]), [2]), [1])[1,1])
	elseif a.typ == "Const and trend"
		println(io, a.typ)
		println(io, mapslices((x)->join(x, "\n"), mapslices((x)->join(x, "  "), map((x)->format(x, precision = 2, signed = true), a.C[e:2,:]), [2]), [1])[1,1])
	end

	for i=1:a.lags
		println(io, line)
		print(io,"Lag ")
		println(io, i)
		println(io, line)
		println(io, mapslices((x)->join(x, "\n"), mapslices((x)->join(x, "  "), map((x)->format(x, precision = 2, signed = true), Clag(a,i)), [2]), [1])[1,1])
	end
end

function simulateVAR(estimate::VARRepresentation, N, burnout=1000)
	exogen = zeros((N+burnout, estimate.vars))
	for i=1:(N + burnout)
		if estimate.typ == "None"
			exogen[i, :] = zeros(size(estimate.C)[2])
		elseif estimate.typ == "Const"
			exogen[i, :] = estimate.C[1,:]
		elseif estimate.typ == "Trend"
			exogen[i, :] = estimate.C[1,:]*i
		elseif estimate.typ == "Const and trend"
			exogen[i, :] = estimate.C[1,:] + estimate.C[2,:]*i
		end
	end
	return simulateVAR(estimate.Σ, [Clag(estimate, i) for i = 1:estimate.lags], exogen, N, burnout)
end

function simulateVAR(Σ, Π, exogen, N, burnout=1000)
	d = MvNormal(Σ)
	errors = rand(d, N+burnout)
	simulation = zeros((N+burnout, size(Σ)[1]))
	for i=(size(Π)[1]+1):(N + burnout)
		simulation[i, :] = exogen[i, :]
		for j=1:(size(Π)[1])
			simulation[i, :] += Π'[j] * simulation[i-j,:]
		end
		simulation[i, :] += errors[:, i]
	end
	return simulation[(burnout+1):(size(simulation)[1]), :]
end

ext: jl


text:
# This file exemplifies the workflow from data input to optimization result generation
using CapacityExpansion
using Clp
using Plots

## LOAD DATA ##
state="GER_18" # or "GER_18" or "CA_1" or "TX_1"
# laod ts-data
ts_input_data = load_timeseries_data_provided(state; T=24) #CEP
# load cep-data
cep_data = load_cep_data_provided(state)

## CLUSTERING ##
# run aggregation with kmeans
ts_clust_data = run_clust(ts_input_data;method="kmeans",representation="centroid",n_init=100,n_clust=5) # Increase n_init to 10000 for a "real" run
## OPTIMIZATION EXAMPLES##
# select optimizer
optimizer=Clp.Optimizer

# Optimize the capacity expansion problem with a co2_limit of 1000
cep = run_opt(ts_clust_data.best_results,cep_data,optimizer;co2_limit=1000)

# Extract the CAP-Variable
cap=cep.variables["CAP"]

# Prepare the data that shall be plotted
dat=sum(cap[:, "new", :], dims=2)
# Prepare the xticks
xticks=(1:length(axes(cap,"tech")),axes(cap,"tech"))

# Plot as a bar-plot
bar(dat,title="Cap", xticks=xticks ,legend=false)


text:
export
ruleSPBernoulliOutNP,
ruleSPBernoulliIn1PN,
ruleSPBernoulliOutNB,
ruleVBBernoulliOut,
ruleVBBernoulliIn1

ruleSPBernoulliOutNP(msg_out::Nothing, msg_p::Message{PointMass, Univariate}) = Message(Univariate, Bernoulli, p=deepcopy(msg_p.dist.params[:m]))

ruleSPBernoulliIn1PN(msg_out::Message{PointMass, Univariate}, msg_p::Nothing) = Message(Univariate, Beta, a=msg_out.dist.params[:m]+1, b=2-msg_out.dist.params[:m])

ruleSPBernoulliOutNB(msg_out::Nothing, msg_p::Message{Beta, Univariate}) = Message(Univariate, Bernoulli, p=msg_p.dist.params[:a]/(msg_p.dist.params[:a] + msg_p.dist.params[:b]))

function ruleVBBernoulliOut(marg_out::Any, marg_p::ProbabilityDistribution{Univariate})
    rho_1 = clamp(exp(unsafeLogMean(marg_p)), tiny, huge)
    rho_2 = clamp(exp(unsafeMirroredLogMean(marg_p)), tiny, huge)

    Message(Univariate, Bernoulli, p=rho_1/(rho_1 + rho_2))
end

ruleVBBernoulliIn1(marg_out::ProbabilityDistribution, marg_p::Any) = Message(Univariate, Beta, a=unsafeMean(marg_out) + 1.0, b= 2.0 - unsafeMean(marg_out))


text:
using UnPack

@testset "Parameters" begin
    # First constructor
    p = Parameters(0.4, 1.4737, 1e-5, Float32)
    @unpack ϕ, kT, N, τ, seed = p
    @test ϕ == 0.4f0
    @test kT == 1.4737f0
    @test N == 512
    @test τ == 0.00001f0
    @test seed == 0

    # Second constructor, define number of particles
    num_particles = 256
    p1 = Parameters(0.4, 1.4737, num_particles, 1e-5, Float32)
    @unpack ϕ, kT, N, τ, seed = p1
    @test ϕ == 0.4f0
    @test kT == 1.4737f0
    @test N == num_particles
    @test τ == 0.00001f0
    @test seed == 0

    # Third constructor, define number of particles and random seed
    num_particles = 256
    new_seed = 804533
    p2 = Parameters{Float32}(0.4, 1.4737, num_particles, 1e-5, new_seed)
    @unpack ϕ, kT, N, τ, seed = p2
    @test ϕ == 0.4f0
    @test kT == 1.4737f0
    @test N == num_particles
    @test τ == 0.00001f0
    @test seed == new_seed
end
@testset "SimulationSystem" begin
    # Create the parameters and simulation system from it
    p = Parameters(0.4, 1.4737, 1e-5, Float64)
    # 3 is the dimension of the system, a 3D simulation system
    dims = 3
    simul_sys = SimulationSystem(p, dims, Float64)
    @unpack params, ρ, boxl, rc, system, energy = simul_sys
    # Check the parameters
    @unpack ϕ, kT, N, τ, seed = p
    @test ϕ == 0.4
    @test kT == 1.4737
    @test N == 512
    @test τ == 0.00001
    @test seed == 0
    # Now check the system parameters
    density = 6.0 * params.ϕ / π
    @test density ≈ ρ
    @test boxl == cbrt(params.N / density)
    @test rc == boxl * 0.5
    @test energy == 0.0
    # Check that the positions and forces are correctly initialized
    total_size = (512, dims)
    @test size(system.positions) == total_size
    @test size(system.forces) == total_size
    @test all(system.positions .== 0.0)
    @test all(system.forces .== 0.0)
end


text:
function write_varattr(io, attr::Dict, gms::oneProblem, attrkey::AbstractString, attrval::AbstractString, funcstrs::Any)

    transpose_attr = Dict(k=>attr[k] for k in keys(attr) if attr[k] == attrkey)
    write_varattr(io, transpose_attr, gms, funcstrs, attrval)
    return
end

function write_varattr(io, attr::Dict, gms::oneProblem, funcstrs::Any, useval=nothing)

    vals = [i[2] for i in attr]
    freqtable = Dict(v => [] for v in unique(vals))
    isa(funcstrs, Vector) ? funcstrs = funcstrs : funcstrs = [funcstrs]

    for i in attr
        push!(freqtable[i[2]], i[1])
    end

    for val in keys(freqtable)
        idxseq, nameseq = separate_varseq(freqtable[val])
        useval == nothing ? valstring = val : valstring = useval
        for i in nameseq
            for fs in funcstrs
                write(io, "$(fs)($(i), $(valstring))\n")
            end
        end
        for i in idxseq
            if length(i) > 3
                write(io, "for i in $(i[1][1]):$(i[end][1])\n")
                for fs in funcstrs
                    write(io, "\t$(fs)($(i[1][2])[i], $(valstring))\n")
                end
                write(io, "end\n")
            else
                for j in i
                    for fs in funcstrs
                        write(io, "$(fs)($(j[2])[$(j[1])], $(valstring))\n")
                    end
                end
            end
        end
    end

    return
end

function write_varattr(io, gms::oneProblem, mode::AbstractString)

    for col in gms.cols

        if mode == "raw"
            colName = col
        elseif mode == "index"
            colName = gms.cols2vars[col]
        end

        if haskey(gms.lb, col)
            write(f, "setlowerbound($(colName), $(gms.lb[col]))\n")
        end
        if haskey(gms.ub, col)
            write(f, "setupperbound($(colName), $(gms.ub[col]))\n")
        end
        if haskey(gms.fx, col)
            write(f, "setlowerbound($(colName), $(gms.fx[col]))\n")
            write(f, "setupperbound($(colName), $(gms.fx[col]))\n")
        end
        if haskey(gms.l, col)
            write(f, "setvalue($(colName), $(gms.l[col]))\n")
        end
        if haskey(gms.prior, col)
            warning("Branching priority indicated in gms file. This behavior is solver dependent in JuMP. Consider implement them for exact gms instruction.")
        end
        if haskey(gms.scale, col)
            error("variable scaling applied in gms file. Don't know how to handle this one.")
        end
        if haskey(gms.stage, col)
            warning("Stage variable attribute indicated in gms file. Generic JuMP doesn't support his. Consider StructJuMP.jl for this implementation.")
        end
    end

    return
end

function separate_varseq(varseq::Vector)

    nameseq = []
    idxseq = []
    for i in varseq
        if match(r"\d+", i) == nothing
            push!(nameseq, i)
        else
            idx = convert(Int, parse(match(r"\d+", i).match))
            vname = match(r"[a-zA-Z]+", i).match
            push!(idxseq, (idx, vname))
        end
    end

    return chop_idx_seq(sort(idxseq)), nameseq
end

function chop_idx_seq(xidxs::Vector)

    if length(xidxs) == 1
        return [xidxs]
    end

    chop = []

    N = length(xidxs)
    num = 1
    cnt = 0
    for i in 2:N
        # Sequence break or variable name can cause a partition
        if xidxs[i][1] - xidxs[i-1][1] > 1 || xidxs[i][2] != xidxs[i-1][2]
            subseq = [xidxs[num:(i-1)];]
            push!(chop, subseq)
            num = i
        end
    end

    push!(chop, [xidxs[num:N];])

    return chop
end

function try_iflinear(c::AbstractString, quadNL::Bool=false)

	contains(c, "sqrt") && return false
	contains(c, "^") && return false
	contains(c, "abs") && return false

    linear = true

    try
        con = eval(parse(c))
        if quadNL && isa(con, JuMP.ConstraintRef{JuMP.Model,JuMP.GenericQuadConstraint{JuMP.GenericQuadExpr{Float64,JuMP.Variable}}})
            linear = false
        end
    catch e
        linear = false
    end

    return linear
end

function try_addvar(m_t, var::AbstractString, var_idx::AbstractString)
    vs = "@variable(m_t, $(var)[$(var_idx)])"
    @eval $(parse(vs))
    return m
end

function try_addvar(m_t, var::AbstractString)
    vs = "@variable(m_t, $(var))"
    @eval $(parse(vs))
    return m
end

function get_one_line(file::IOStream; one_line=" ")
    while one_line[end] != ';'
        one_line = string(one_line, lstrip(readline(file),' '))
    end
    one_line = split(one_line, ';', keep=false)
    @assert length(one_line) == 1
    return one_line[1]
end

function replace_vars(gms::oneProblem)
    for row in gms.rows
		gms.rowsLHS[row] = replace(gms.rowsLHS[row], r"x(\d+)", s"x[\1]")
		gms.rowsLHS[row] = replace(gms.rowsLHS[row], r"i(\d+)", s"i[\1]")
		gms.rowsLHS[row] = replace(gms.rowsLHS[row], r"b(\d+)", s"b[\1]")
    end
    gms.objective = replace(gms.objective, r"x(\d+)", s"x[\1]")
    gms.objective = replace(gms.objective, r"i(\d+)", s"i[\1]")
    gms.objective = replace(gms.objective, r"b(\d+)", s"b[\1]")
end

function replace_oprs(line::AbstractString; kwargs...)

    line = replace(line, "\n", "")
    line = replace(line, " ", "")
    contains(line, "sqr") && (line = _replace_sqr(line))
    contains(line, "POWER") && (line = _replace_POWER(line))
	contains(line, "power") && (line = _replace_POWER(line,powerstring="power("))
    line = replace(line, "**", "^")

	unsupported_opr = ["arctan", "ARCTAN",
					   "ceil", "CEIL",
					   "errorof", "ERROROF",
					   "floor", "FLOOR",
					   "mapVal", "MAPVAL",
					   "max", "MAX",
					   "min", "MIN",
					   "mod", "MOD",
					   "normal", "NORMAL",
					   "round", "ROUND",
					   "sign", "SIGN",
					   "trunc", "TRUNC",
					   "uniform", "UNIFORM"]
	for o in unsupported_opr
		contains(line, o) && error("Unsupported operator $(o) detected")
	end

    return line
end

function _replace_sqr(str::AbstractString)
    # @show "Begin ",str
    for i in 1:(length(str)-4)
        if str[i:(i+3)] == "sqr("
            warp = 0
            for j in (i+4):length(str)
                if str[j] == '('
                    warp += 1
                elseif str[j] == ')'
                    warp -= 1
                end
                if warp == -1
                    substart = i+4
                    subclose = j-1
                    subStr = str[substart:subclose]
                    # @show "Before replace ", subStr
                    str = replace(str, "sqr($subStr)"," ($subStr)^2")
                    # @show "After replace ", str
                    subStr = str[substart-2:subclose-2]
                    subStr = _replace_sqr(subStr)
                    # @show "Processed subStr = ", subStr
                    str = string(str[1:substart-3], subStr, str[subclose-1:end])
                    # @show "After recomposing", length(str)
                    break
                end
            end
        end
    end
    # @show "exiting this level with ", str
    return str
end

function _replace_POWER(str::AbstractString; powerstring="POWER(")
    pVal = 0.0
    for i in 1:length(str)-6
        if str[i:(i+5)] == powerstring
            warp = 1
            pValStart = 0
            for j in (i+6):length(str)
                if str[j] == '('
                    warp += 1
                elseif str[j] == ')'
                    warp -= 1
                elseif str[j] == ','
                    pValStart = j
                end
                if warp == 0
                    @assert pValStart > 0   # By the time a function is closed, pVal segment must be detected
                    substart = i+6
                    subclose = j-1
                    subStr = str[substart:subclose]
                    argStr = str[substart:(pValStart-1)]
                    pVal = parse(str[pValStart+1:subclose])   # Don't allow complicated expression here
					str = replace(str, "$(powerstring)$(subStr))"," ($argStr)^$(pVal)")
                    subStr = str[substart-5:subclose-5]
                    subStr = _replace_POWER(subStr)
                    str = string(str[1:substart-6], subStr, str[subclose-4:end])
                    i = substart + 1 # Dail back the looper due to the reduced length of string
                    break
                end
            end
        end
        i > length(str)-6 && break
    end

    return str
end

function read_command(file::IOStream, gms::oneProblem, lInit::AbstractString; kwargs...)
    warn("I am not that smart to parse command in gms files, YET")
    return 0
end

function convert_equality(probname="")

    println("This function handles problem with too many equality constraints.", prefix="POD Experiment: ")
    if !isfile("$(Pkg.dir())/toJuMP/instances/$(probname).jl")
        println("NO problem file detected in $(Pkg.dir())/toJuMP/instances/$(probname).jl", prefix="POD Experiment: ")
        return
    end

    f = open("$(Pkg.dir())/toJuMP/instances/$(probname).jl", "r")
    outf = open("$(Pkg.dir())/toJuMP/instances/$(probname)_gl.jl", "w")

    for l in readlines(f)
        if ismatch(r"==", l)
            geq = replace(l, "==", ">=")
            leq = replace(l, "==", "<=")
            write(outf, geq)
            write(outf, "\n")
            write(outf, leq)
            write(outf, "\n")
        elseif ismatch(Regex(probname), l)
            write(outf, "function $(probname)_gl(;options=Dict())\n")
        else
            write(outf, l)
            write(outf, "\n")
        end
    end

    close(f)
    close(outf)

    return
end

function parse_varname(gms::oneProblem)

    if isempty(gms.vars)

        for varString in gms.cols
            varName = split(varString, r"\d+")[1]
            if varName != varString
                varSplit = split(varString, varName)
                varIndex = parse(varSplit[2])
                @assert isa(varIndex, Int)
                if !haskey(gms.vars, varName)
                    gms.vars[varName] = []
                end
                if varIndex in gms.vars[varName]
                    error("ERROR|gms2jump.jl|parse_varname()|Conflicting indice variable names")
                end
                push!(gms.vars[varName], varIndex)
                gms.cols2vars[varString] = parse(string(varName, "[",varIndex,"]"))
                gms.vars2cols[parse(string(varName,"[",varIndex,"]"))] = varString
            else
                if !haskey(gms.vars, varName)
                    gms.vars[varName] = 0
                else
                    error("ERROR|gms2jump.jl|parse_varname()|Conflicting symbolic variable names.")
                end
                gms.cols2vars[varString] = parse(varName)
                gms.vars2cols[parse(varName)] = varString
            end

        end
    else
        return 0
    end

    return
end

function clear_m_tester()
    cleanstring = "m_tester = Model()"
    eval(parse(cleanstring))
    return
end


text:
using Pkg
Pkg.activate(".")
using Psychometrics
using Distributions
using LinearAlgebra
using Dates
using Random
#using StatsPlots

function est()
    I_total = 40
    N = 500

    # ITEM PARAMETERS AND LATENTS
    items = [Item(i, string("item_",i), ["math"], Parameters2PL(Product([LogNormal(0.3, 0.2), Normal(0,1)]), [1e-5,5.0], [-6.0, 6.0])) for i = 1 : I_total];
    examinees = [Examinee(e, string("examinee_",e), Latent1D(Normal(0,1), [-6.0, 6.0])) for e = 1 : N];

    # RESPONSES

    responses = answer(examinees, items);

    # Get response matrix (zeros are for wrong answers and for missing answers)
    response_matrix = get_response_matrix(responses, I_total, N);

    #Set Seeds for Random Generation
    #Random.seed!()

    # function tr_norm_gen(number, mu, v)
    #     temp_u = Distributions.rand(Distributions.Uniform(zero(Float64),one(Float64)), number)
    #     norm_dist = Distributions.Normal(mu, sqrt(v))
    #     norm_cdf = Distributions.cdf(norm_dist, zero(Float64))
    #     tr_norm_value = temp_u .- (temp_u .* norm_cdf) .+ norm_cdf
    #     return Distributions.quantile(norm_dist, tr_norm_value)
    # end

    # using StatsPlots
    # A short example of sampling truncated standard normal distribution
    # histogram(tr_norm_gen(800,0,1))
    # #histogram(truncate_rand(Normal(0.0,1.0), [0.0, Inf], n=800))
    # histogram(rand(Distributions.TruncatedNormal(0.0,1.0, 0.0, Inf), 800))

    # Simulation code by Jiang and Templin
    # Gibbs Samplers for Logistic Item Response Models via the Pólya–GammaDistribution:
    # A Computationally Efficient Data-Augmentation StrategyArticleinPsychometrika
    # October 2018 DOI: 10.1007/s11336-018-9641-xCITATIONS2READS692

    ##################################################################################################
    ###########################      Data and Settings    ############################################
    ##################################################################################################

    #Specified by users:Response matrix X
    X = permutedims(response_matrix);
    ##################################################################################################
    ###########################      PolyGamma MCMC sampler   ########################################
    ##################################################################################################

    #K matrix for the coming PolyaGamma samplings
    K = X .- 0.5;

    #Initial Values, need to make sure that all Variance needs to be positive
    ini_a = ones(I_total);
    ini_b = zeros(I_total);
    ini_theta = zeros(N);
    ini_kernel = zeros(N, I_total);

    ini_w = [Distributions.rand(Psychometrics.PolyaGamma(1.0, ini_a[i] * (ini_theta[p] - ini_b[i]))) for p = 1:N, i = 1:I_total] ;

    function theta_sampler(
        vec_a, # I_total x 1
        vec_b, # I_total x 1
        mat_w_id, # N x 1
        mat_k_id, # N x 1
        theta_prior_mu, # 1 X 1
        theta_prior_var,# 1 X 1
    )
        omega_theta = LinearAlgebra.Diagonal(mat_w_id) # I_total x I_total
        theta_variance = 1 / ((vec_a' * omega_theta * vec_a) + (1 / theta_prior_var)) # 1 x 1
        z_theta = vec_a .* vec_b .* mat_w_id + mat_k_id # I_total x 1 # times w_e
        theta_mu = theta_variance * ((vec_a' * z_theta) + (theta_prior_mu / theta_prior_var))
        rand(Distributions.Normal(theta_mu, sqrt(theta_variance)))
    end
    function a_sampler(
        vec_theta, # N x 1
        vec_b_id, # I_total x 1
        mat_w_id, # N x I_total
        mat_k_id, # N x I_total
        a_prior_mu, # 1 x 1
        a_prior_var, # 1 x 1
    )
        omega_a = LinearAlgebra.Diagonal(mat_w_id) # N x N
        b_omega_a = vec_theta .- vec_b_id
        a_variance =
            1 / (
                (b_omega_a' * omega_a * b_omega_a) +
                (1 / a_prior_var)
            ) # 1 x 1
        #z_a = mat_k_id ./ mat_w_id # don't need it
        a_mu = a_variance * (b_omega_a' * mat_k_id + (a_prior_mu / a_prior_var)) # ok
        #return tr_norm_gen(1, a_mu, a_variance)[1]
        return  rand(Distributions.TruncatedNormal(a_mu, sqrt(a_variance), 0.0, Inf))
    end
    #  ( T   -   b   ) ^2 *  w
    # N x 1   I_total x 1        N x I_total

    function b_sampler(vec_theta, vec_a_id, mat_w_id, mat_k_id, b_prior_mu, b_prior_var, N)
        omega_b = LinearAlgebra.Diagonal(mat_w_id) # N x N
        vec_a_id_vec = fill(vec_a_id, N) #N
        b_variance = 1 / (vec_a_id_vec' * omega_b * vec_a_id_vec + (1 / b_prior_var))  # 1 x 1 ok, vec_a_id_vec can be positive because it is squared anyway
        z_b = mat_k_id .- (vec_a_id .* vec_theta .* mat_w_id) # N x 1 # times w_i
        b_mu = b_variance * (((.-vec_a_id_vec)' * z_b) + (b_prior_mu / b_prior_var)) # ok
        return rand(Distributions.Normal(b_mu, sqrt(b_variance)))
    end
    function i_sampler(vec_theta, vec_a_id, vec_b_id, mat_w_id, mat_k_id, a_prior_mu, a_prior_var, b_prior_mu, b_prior_var)
        sigma2 = mapreduce(
            (e, w) -> [(e - vec_b_id)^2, vec_a_id^2 ] .* w,
            +,
            vec_theta,
            mat_w_id,
        )
        sigma2 = 1 ./ (sigma2 + (1 ./ [a_prior_var, b_prior_var]))
        mu = mapreduce(
            (e, w, r) -> [
                (e - vec_b_id) *
                (r),
                -vec_a_id * (
                    #(get_responses_by_examinee_id(e.id, responses)[1].val - 0.5) -
                    (r) -
                    (vec_a_id * e * w)
                ),
            ],
            +,
            vec_theta,
            mat_w_id,
            mat_k_id
        )
        mu =
            sigma2 .* (
                mu + (
                    [a_prior_mu, b_prior_mu] ./
                    [a_prior_var, b_prior_var]
                )
            )
        return rand(Distributions.Product([
            Distributions.TruncatedNormal(mu[1], sqrt(sigma2[1]), 0.0, Inf),
            Distributions.Normal(mu[2], sqrt(sigma2[2])),
        ]))
    end
    # using RCall
    # R"""
    # rpg_sp <- function(x){
    # BayesLogit::rpg(1,1,x)
    # }
    # """
    function w_sampler(vec_theta, vec_a, vec_b, N, I_total)
        #w_temp = [vec_a[i] * (vec_theta[p] - vec_b[i]) for p = 1:N, i = 1:I_total]
        #pg = rcopy(R"apply($w_temp, c(1,2), rpg_sp)")
        pg =[rand(PolyaGamma(1.0, vec_a[i] * (vec_theta[p] - vec_b[i]))) for p = 1:N, i = 1:I_total]
        return pg
    end

    ##################################################################################################
    ###########################   Estimation Starts from here  #######################################
    ##################################################################################################

    #Specified by users:Construct 2PL IRT Estimator

    Iter = 4000
    theta_prior_mu = 0
    theta_prior_var = 1
    a_prior_mu = 1
    a_prior_var = 5
    b_prior_mu = 0
    b_prior_var = 5

    sav_a = zeros(Float64, Iter, I_total);
    sav_b = zeros(Float64, Iter, I_total);
    sav_theta = zeros(Float64, Iter, N);

    sav_a[1, :] = ini_a
    sav_b[1, :] = ini_b
    sav_theta[1, :] = ini_theta
    #sav.iter.hist<-matrix(0,N,I_total)

    for iter = 2:Iter
        sav_w = w_sampler(sav_theta[iter-1, :], sav_a[iter-1, :], sav_b[iter-1, :], N, I_total)
        if mod(iter,100)==0
            println(iter)
        end
        sav_theta_iter = map(
            p -> begin
                no_missing_i = findall(.!ismissing.(K[p, :]))
                theta_sampler(
                sav_a[iter-1, no_missing_i],
                sav_b[iter-1, no_missing_i],
                sav_w[p, no_missing_i],
                K[p, no_missing_i],
                theta_prior_mu,
                theta_prior_var,
            )
        end
        ,
        1:N
        )
        res_iter = map(
            i -> begin
                no_missing_e = findall(.!ismissing.(K[:, i]))
                i_sampler(
                    sav_theta_iter,
                    sav_a[iter-1, i],
                    sav_b[iter-1, i],
                    sav_w[no_missing_e, i],
                    K[no_missing_e, i],
                    a_prior_mu,
                    a_prior_var,
                    b_prior_mu,
                    b_prior_var
                )
            end,
            1:I_total,
        )
        sav_a_iter = map( r -> r[1], res_iter)
        sav_b_iter = map( r -> r[2], res_iter)

        # sav_a_iter = map(
        #     i -> a_sampler(
        #         sav_theta_iter,
        #         sav_b[iter-1, i],
        #         sav_w[:, i],
        #         K[:, i],
        #         a_prior_mu,
        #         a_prior_var,
        #     ),
        #     1:I_total,
        # )
        # sav_b_iter = map(
        #     i -> b_sampler(
        #         sav_theta_iter,
        #         sav_a_iter[i],
        #         sav_w[:, i],
        #         K[:, i],
        #         b_prior_mu,
        #         b_prior_var,
        #         N,
        #     ),
        #     1:I_total,
        # )
        sav_theta[iter, :] .= copy(sav_theta_iter)
        sav_a[iter, :] .= copy(sav_a_iter)
        sav_b[iter, :] .= copy(sav_b_iter)
    end

    sampled = (mcmc_iter-1000):mcmc_iter
    mean_a = [mean(sav_a[sampled, i]) for i = 1:I_total];
    mean_b = [mean(sav_b[sampled, i]) for i = 1:I_total];
    mean_theta = [mean(sav_theta[sampled, i]) for i = 1:N];

    # hcat(map(i -> i.parameters.a, items), mean_a)
    # hcat(map(i -> i.parameters.b, items), mean_b)

    # RMSEs
    println(sqrt(sum((map(i -> i.parameters.a, items) .- mean_a).^2)/I_total))
    println(sqrt(sum((map(i -> i.parameters.b, items) .- mean_b).^2)/I_total))
    println(sqrt(sum((map(e -> e.latent.val, examinees) .- mean_theta).^2)/N))

    return examinees, examinees_est, items, items_est, responses
end
examinees, examinees_est, items, items_est, responses = est();

# using StatsPlots
# plot(sav_a[:,1])

# RMSE_a_500 = zeros(Iter-500);
# RMSE_b_500 = zeros(Iter-500);
# RMSE_t_500 = zeros(Iter-500);
# for j in 1:(Iter-500)
# mean_a = [mean(sav_a[j:(j+499), i]) for i = 1:I_total];
# mean_b = [mean(sav_b[j:(j+499), i]) for i = 1:I_total];
# mean_theta = [mean(sav_theta[j:(j+499), i]) for i = 1:N];
# # RMSEs
# RMSE_a_500[j]=sqrt(sum((map(i -> i.parameters.a, items) .- mean_a).^2)/I_total)
# RMSE_b_500[j]=sqrt(sum((map(i -> i.parameters.b, items) .- mean_b).^2)/I_total)
# RMSE_t_500[j]=sqrt(sum((map(e -> e.latent.val, examinees) .- mean_theta).^2)/N)
# end

# plot(hcat(RMSE_a_500, RMSE_b_500, RMSE_t_500))

# RMSE_a_100 = zeros(Iter-100);
# RMSE_b_100 = zeros(Iter-100);
# RMSE_t_100 = zeros(Iter-100);
# for j in 1:(Iter-100)
# mean_a = [mean(sav_a[j:(j+99), i]) for i = 1:I_total];
# mean_b = [mean(sav_b[j:(j+99), i]) for i = 1:I_total];
# mean_theta = [mean(sav_theta[j:(j+99), i]) for i = 1:N];
# # RMSEs
# RMSE_a_100[j]=sqrt(sum((map(i -> i.parameters.a, items) .- mean_a).^2)/I_total)
# RMSE_b_100[j]=sqrt(sum((map(i -> i.parameters.b, items) .- mean_b).^2)/I_total)
# RMSE_t_100[j]=sqrt(sum((map(e -> e.latent.val, examinees) .- mean_theta).^2)/N)
# end

# plot(hcat(RMSE_a_100, RMSE_b_100, RMSE_t_100))

max_stars_repo_path: examples/Jiang Templin/MCMC sampler.jl

max_issues_repo_path: examples/Jiang Templin/MCMC sampler.jl

max_forks_repo_path: examples/Jiang Templin/MCMC sampler.jl


text:
(
    runtime = now(),
    database = (
        connection = (
            url = "https://secureserver.me/staging",
            port = 3601,
        ),
        credentials = (
            username = "stagingadmin",
            password = "",
            empty = ()
        ),
    ),
    empty = Dict(),
)

max_stars_repo_name: citkane/Config

max_issues_repo_name: citkane/Config

max_forks_repo_name: citkane/Config


text:
module layout

import Term: int, get_last_valid_str_idx
import ..renderables: RenderablesUnion, Renderable, AbstractRenderable
import ..measure: Measure
import ..segment: Segment
using ..box
import ..consoles: console_width, console_height

export Padding, vstack, hstack
export Spacer, vLine, hLine

# ---------------------------------------------------------------------------- #
#                                    PADDING                                   #
# ---------------------------------------------------------------------------- #
"""
    Padding

Stores empty string to pad a string to the desired with.
"""
struct Padding
    left::AbstractString
    right::AbstractString
end

"""
    Padding(text::AbstractString, target_width::Int, method::Symbol)::Padding

Create `Padding` for a string given a justification `method` ∈ (:left, :center, :right).
"""
function Padding(text::AbstractString, target_width::Int, method::Symbol)::Padding
    # get total padding size
    lw = Measure(text).w
    if lw == target_width
        return Padding("", "")
    else
        if lw < target_width
            @debug "\e[31mTarget width is $target_width but the text has width $lw: \e[0m\n   $text\n     Cleaned: '$(remove_ansi(remove_markup(text)))'"
        end
    end

    if target_width < lw
        @debug "While creating padding, target width < lw" target_width lw method text
        return Padding("", "")
    end
    padding = " "^(target_width - lw - 1)
    pad_width = Measure(padding).w

    # split left/right padding for left/right justify
    if pad_width > 0
        cut = pad_width % 2 == 0 ? Int(pad_width / 2) : (Int ∘ floor)(pad_width / 2)
        left, right = padding[1:cut], padding[(cut + 1):end]
        @assert length(left * right) <= length(padding) "\e[31m$(length(left * right)) instead of $(length(padding))"
    else
        left, right = "", ""
    end

    # craete padding
    if method == :center
        padding = Padding(" " * left, right)
    elseif method == :left
        padding = Padding(" ", padding)
    elseif method == :right
        padding = Padding(padding, " ")
    end

    # @info "made padding" text target_width method padding lw pad_width left right
    @assert Measure(padding.left * text * padding.right).w <= target_width "\e[31mPadded width $(Measure(padding.left * text * padding.right).w), target: $target_width $padding"

    return padding
end

function Base.show(io::IO, padding::Padding)
    return print(
        io,
        "$(typeof(padding))  \e[2m(left: $(length(padding.left)), right: $(length(padding.right)))\e[0m",
    )
end

# ---------------------------------------------------------------------------- #
#                                   STACKING                                   #
# ---------------------------------------------------------------------------- #
"""
    vstack(r1::RenderablesUnion, r2::RenderablesUnion)

Vertically stack two renderables to give a new renderable.
"""
function vstack(r1::RenderablesUnion, r2::RenderablesUnion)
    r1 = Renderable(r1)
    r2 = Renderable(r2)

    # get dimensions of final renderable
    w1 = r1.measure.w
    w2 = r2.measure.w

    # create segments stack
    segments::Vector{Segment} = vcat(r1.segments, r2.segments)
    measure = Measure(max(w1, w2), length(segments))

    return Renderable(segments, measure)
end

"""
    vstack(renderables...)

Vertically stack a variable number of renderables
"""
function vstack(renderables...)
    renderable = Renderable()

    for ren in renderables
        renderable = vstack(renderable, ren)
    end
    return renderable
end

"""
    hstack(r1::RenderablesUnion, r2::RenderablesUnion)

Horizontally stack two renderables to give a new renderable.
"""
function hstack(r1::RenderablesUnion, r2::RenderablesUnion)
    r1 = Renderable(r1)
    r2 = Renderable(r2)

    # get dimensions of final renderable
    h1 = r1.measure.h
    h2 = r2.measure.h

    # make sure both renderables have the same number of segments
    Δh = abs(h2 - h1)
    if h1 > h2
        r2.segments = vcat(r2.segments, [Segment(" "^r2.measure.w) for i in 1:Δh])
    elseif h1 < h2
        r1.segments = vcat(r1.segments, [Segment(" "^r1.measure.w) for i in 1:Δh])
    end

    # combine segments
    segments = [Segment(s1.text * s2.text) for (s1, s2) in zip(r1.segments, r2.segments)]

    return Renderable(segments, Measure(segments))
end

"""
    hstack(renderables...)

Horizonatlly stack a variable number of renderables
"""
function hstack(renderables...)
    renderable = Renderable()

    for ren in renderables
        renderable = hstack(renderable, ren)
    end
    return renderable
end

# --------------------------------- operators -------------------------------- #

"""
    Operators for more coincise layout of renderables
"""

Base.:/(r1::RenderablesUnion, r2::RenderablesUnion) = vstack(r1, r2)

Base.:*(r1::AbstractRenderable, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractString, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractRenderable, r2::AbstractString) = hstack(r1, r2)

# ---------------------------------------------------------------------------- #
#                                LINES & SPACER                                #
# ---------------------------------------------------------------------------- #
abstract type AbstractLayoutElement <: AbstractRenderable end

# ---------------------------------- spacer ---------------------------------- #
"""
        Spacer

A box of empty text with given width and height.
"""
mutable struct Spacer <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

function Spacer(width::Number, height::Number; char::Char = ' ')
    width = int(width)
    height = int(height)

    line = char^width
    segments = [Segment(line) for i in 1:height]
    return Spacer(segments, Measure(segments))
end

# ----------------------------------- vline ---------------------------------- #
"""
    vLine

A multi-line renderable with each line made of a | to create a vertical line
"""
mutable struct vLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
    height::Int
end

"""
    vLine(height::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a `vLine` given a height and, optionally, style information.
"""
function vLine(
    height::Number; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED
)
    height = int(height)
    char = string(eval(box).head.left)
    segments = [Segment(char, style) for i in 1:height]
    return vLine(segments, Measure(segments), height)
end

"""
    vLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Create a `vLine` as tall as the `stdout` console
"""
function vLine(; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED)
    return vLine(console_height(); style = style, box = box)
end

"""
    hLine

A 1-line renderable made of repeated character from a Box.
"""
mutable struct hLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
    width::Int
end

"""
    hLine(width::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a styled `hLine` of given width.
"""
function hLine(
    width::Number; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED
)
    width = int(width)
    char = eval(box).row.mid
    segments = [Segment(char^width, style)]
    return hLine(segments, Measure(segments), width)
end

"""
    hLine(width::Number, text::String; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Creates an hLine object with texte centered horizontally.
"""
function hLine(
    width::Number,
    text::String;
    style::Union{String,Nothing} = nothing,
    box::Symbol = :ROUNDED,
)
    box = eval(box)

    initial_line = box.top.mid^width

    cutval = int(ncodeunits(initial_line) / 2 - ncodeunits(text) - 5)
    cut_start = get_last_valid_str_idx(initial_line, cutval)

    pre = Segment(
        Segment(initial_line[1:cut_start], style) *
        "\e[0m" *
        " " *
        Segment(text, style) *
        " ",
    )

    post = Segment(box.top.mid^(length(initial_line) - pre.measure.w - 1), style)

    segments = [Segment(pre * (post), style)]
    return hLine(segments, Measure(segments), width)
end

"""
    hLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout`
"""
function hLine(; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED)
    return hLine(console_width(); style = style, box = box)
end

"""
    hLine(text::AbstractString; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout` with centered text.
"""
function hLine(
    text::AbstractString; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED
)
    return hLine(console_width(), text; style = style, box = box)
end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 32be2dac-7618-11eb-2ea5-9f165c3c91e1
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 334f2774-7610-11eb-197e-eb4c20530a0c
md"""# Fruitful Functions

Many of the Julia functions we have used, such as the math functions, produce return values. But the functions we’ve written are all void: they have an effect, like printing a value or moving a turtle, but they return nothing. In this chapter you will learn to write fruitful functions."""

# ╔═╡ b1a8c332-7610-11eb-1b18-f7bf7f28ea27
md"""## Return Values

Calling a function generates a return value, which we usually assign to a variable or use as part of an expression:

```julia
e = exp(1.0)
height = radius * sin(radians)
```

!!! tip
    Julia provides the Euler's number ``\mathcal e`` by the symbol `ℯ` (**`\euler TAB`**).

The functions we have written so far are void. Speaking casually, they have no return value; more precisely, their return value is nothing. In this chapter, we are (finally) going to write fruitful functions. The first example is area, which returns the area of a circle with the given radius:

```julia
function area(radius)
	a = π * radius^2
	return a
end
```

We have seen the return statement before, but in a fruitful function the return state‐ ment includes an expression. This statement means: “Return immediately from this function and use the following expression as a return value.” The expression can be arbitrarily complicated, so we could have written this function more concisely:

```julia
function area(radius)
	π * radius^2
end
```

where the last expression is automatically returned; or

```julia
area(radius) = π * radius^2
```

but this is only usefull for a function consisting of a oneline expression.

However, *temporary variables* like `a` and explicit return statements can make debugging easier.

The value returned by a function is the value of the last expression evaluated, which, by default, is the last expression in the body of the function definition.

Sometimes it is useful to have multiple return statements, one in each branch of a conditional:

```julia
function absvalue(x)
	if x < 0
		return -x
	else
		return x
	end
end
```

or

```julia
function absvalue(x)
	if x < 0
		-x
	else
		x
	end
end
```

because the last expression evaluated is automatically returned.

Since these return statements are in an alternative conditional, only one runs.

As soon as a return statement runs, the function terminates without executing any subsequent statements. Code that appears after a return statement, or any other place the flow of execution can never reach, is called *dead code*.

In a fruitful function, it is a good idea to ensure that every possible path through the program hits a return statement. Consider this example:

```julia
function absvalue(x)
	if x < 0
		return -x
	end
	if x > 0
		return x
	end
end
```

This function is incorrect because if `x` happens to be `0`, neither condition is true, and the function ends without hitting a return statement. If the flow of execution gets to the end of a function, the return value is `nothing`, which is not the absolute value of `0`.

```julia
julia> show(absvalue(0))
nothing
```

!!! tip
    Julia provides a built-in function called `abs` that computes absolute values.

"""

# ╔═╡ 4219f3ac-7612-11eb-3167-f55e643f9235
md"""#### Exercise 6-1

Write a compare function that takes two values, `x` and `y`, and returns `1` if `x > y`, `0` if `x == y`, and `-1` if `x < y`.
"""

# ╔═╡ 8b79bba4-7682-11eb-2ea1-0986a9911b40
md"""!!! languages
    In Python and the C programming language, `return` cannot be omitted, neither is the last evaluated expression returned automatically.

    ```python
    def absvalue(x):
        if x < 0:
            return -x
        else:
            return x

    ```

    The return type has to be provided in C:

    ```c
    double absvalue(double x) {
        if (x < 0) {
            return -x;
        } else {
            return x;
        }
    }
    ```

    Matlab doesn't use the `return` keyword but the return values have to be specified in the header:

    ```matlab
    function [ret] = absvalue(x)
        if x < 0
            ret = -x;
        else
            ret = x;
        end
    end
    ```

    A function in Matlab has also to be defined in a file having the same name as the function, in this case *`absvalue.m`*.
"""

# ╔═╡ bfecb80e-7612-11eb-0633-035cbcd3e426
md"""## Incremental Development

As you write larger functions, you might find yourself spending more time debugging.

To deal with increasingly complex programs, you might want to try a process called *incremental development*. The goal of incremental development is to avoid long debugging sessions by adding and testing only a small amount of code at a time.

As an example, suppose you want to find the distance between two points, given by the coordinates ``(x_1, y_1)`` and ``(x_2, y_2)``. By the Pythagorean theorem, the distance is:

```math
d=\sqrt{(x_2−x_1)^2+(y_2−y_1)^2}
```

The first step is to consider what a distance function should look like in Julia. In other words, what are the inputs (parameters) and what is the output (return value)?

In this case, the inputs are two points, which you can represent using four numbers. The return value is the distance represented by a floating-point value.

Immediately you can write an outline of the function:

```julia
function distance(x₁, y₁, x₂, y₂)
	0.0
end
```

Obviously, this version doesn’t compute distances; it always returns zero. But it is syntactically correct, and it runs, which means that you can test it before you make it more complicated. The subscript numbers are available in the Unicode character encoding (**`\_1 TAB`**, **`\_2 TAB`**, etc.).

To test the new function, call it with sample arguments:

```julia
julia> distance(1.0, 2.0, 4.0, 6.0)
0.0
```

I chose these values so that the horizontal distance is ``3`` and the vertical distance is ``4``; that way, the result is ``5``, the hypotenuse of a 3-4-5 triangle. When testing a function, it is useful to know the right answer.

At this point we have confirmed that the function is syntactically correct, and we can start adding code to the body. A reasonable next step is to find the differences ``x_2 − x_1`` and ``y_2 − y_1``. The next version stores those values in temporary variables and prints them with the `@show` macro:

```julia
function distance(x₁, y₁, x₂, y₂)
	dx = x₂ - x₁
	dy = y₂ - y₁
	@show dx dy
	0.0
end
```

If the function is working, it should display `dx = 3` and `dy = 4`. If so, we know that the function is getting the right arguments and performing the first computation correctly. If not, there are only a few lines to check.

Next, we compute the sum of squares of `dx` and `dy`:

```julia
function distance(x₁, y₁, x₂, y₂)
	dx = x₂ - x₁
	dy = y₂ - y₁
	d² = dx^2 + dy^2
	@show d²
	0.0
end
```

Again, you would run the program at this stage and check the output (which should be `25`). Superscript numbers are also available (**`\^2 TAB`**). Finally, you can use `sqrt` or `√`(**`\sqrt TAB`**) to compute and return the result:

```julia
function distance(x₁, y₁, x₂, y₂)
	dx = x₂ - x₁
	dy = y₂ - y₁
	d² = dx^2 + dy^2
	√d²
end
```

If that works correctly, you are done. Otherwise, you might want to print the value of `√d²` before the return statement.

The final version of the function doesn’t display anything when it runs; it only returns a value. The print statements we wrote are useful for debugging, but once you get the function working, you should remove them. Code like that is called *scaffolding* because it is helpful for building the program but is not part of the final product.

When you start out, you should add only a line or two of code at a time. As you gain more experience, you might find yourself writing and debugging bigger chunks. Either way, incremental development can save you a lot of debugging time.

The key aspects of the process are:

1. Start with a working program and make small, incremental changes. At any point, if there is an error, you should have a good idea where it is.
2. Use variables to hold intermediate values so you can display and check them.
3. Once the program is working, you might want to remove some of the scaffolding or consolidate multiple statements into compound expressions, but only if it does not make the program difficult to read.
"""

# ╔═╡ f5edc4c2-7613-11eb-0765-cba4f0286212
md"""#### Exercise 6-2

Use incremental development to write a function called hypotenuse that returns the length of the hypotenuse of a right triangle given the lengths of the other two legs as arguments. Record each stage of the development process as you go."""

# ╔═╡ 47e45c90-7615-11eb-0fe7-358b9d219497
md""" ## Composition

As you should expect by now, you can call one function from within another. As an example, we’ll write a function that takes two points, the center of the circle and a point on the perimeter, and computes the area of the circle.

Assume that the center point is stored in the variables `xc` and `yc`, and the perimeter point is in `xp` and `yp`. The first step is to find the radius of the circle, which is the distance between the two points. We just wrote a function, distance, that does that:

```julia
radius = distance(xc, yc, xp, yp)
```

The next step is to find the area of a circle with that radius; we just wrote that, too:

```julia
result = area(radius)
```

Encapsulating these steps in a function, we get:

```julia
function circlearea(xc, yc, xp, yp)
	radius = distance(xc, yc, xp, yp)
	result = area(radius)
	return result
end
```

The temporary variables `radius` and `result` are useful for development and debugging, but once the program is working, we can make it more concise by composing the function calls:

```julia
circlearea(xc, yc, xp, yp) = area(distance(xc, yc, xp, yp))
```
"""

# ╔═╡ a74ddcba-7615-11eb-1ba0-c1d9a4eacfb5
md"""## Boolean Functions

Functions can return Booleans, which is often convenient for hiding complicated tests inside functions. For example:

```julia
function isdivisible(x, y)
	if x % y == 0
		return true
	else
		return false
	end
end
```

It is common to give Boolean functions names that sound like yes/no questions; `isdivisible` returns either `true` or `false` to indicate whether `x` is divisible by `y`.

Here is an example:

```julia
julia> isdivisible(6, 4)
false
julia> isdivisible(6, 3)
true
```

The result of the `==` operator is a Boolean, so we can write the function more con‐ cisely by returning it directly:

```julia
isdivisible(x, y) =  x % y == 0
```

Boolean functions are often used in conditional statements:

```julia
if isdivisible(x, y)
	println("x is divisible by y")
end
```

It might be tempting to write something like:

```julia
if isdivisible(x, y) == true
	println("x is divisible by y")
end
```

But the extra comparison with `true` is unnecessary.
"""

# ╔═╡ 3a20604e-7616-11eb-20c6-1fb57dbbb3c8
md"""#### Exercise 6-3
Write a function `isbetween(x, y, z)` that returns `true` if `x ≤ y ≤ z` or `false` otherwise.
"""

# ╔═╡ 53822bd0-7616-11eb-08df-0b6ef72120a8
md"""## More Recursion

We have only covered a small subset of Julia, but you might be interested to know that this subset is a *complete* programming language, which means that anything that can be computed can be expressed in this language. Any program ever written could be rewritten using only the language features you have learned so far (actually, you would need a few commands to control devices like the mouse, disks, etc., but that’s all).

Proving that claim is a nontrivial exercise first accomplished by Alan Turing, one of the first computer scientists (some would argue that he was a mathematician, but a lot of early computer scientists started as mathematicians). Accordingly, it is known as the Turing thesis. For a more complete (and accurate) discussion of the Turing thesis, I recommend Michael Sipser’s book *Introduction to the Theory of Computation* (Cengage).

To give you an idea of what you can do with the tools you have learned so far, we’ll evaluate a few recursively defined mathematical functions. A recursive definition is similar to a circular definition, in the sense that the definition contains a reference to the thing being defined. A truly circular definition is not very useful:

*vorpal*:
An adjective used to describe something that is vorpal.

If you saw that definition in the dictionary, you might be annoyed. On the other hand, if you looked up the definition of the factorial function, denoted with the symbol ``!``, you might get something like this:

```math
n!=\begin{cases}
1 & \textrm{if }n=0\\
n(n-1)! & \textrm{if }n>0
\end{cases}
```

This definition says that the factorial of ``0`` is ``1``, and the factorial of any other value, ``n``, is ``n`` multiplied by the factorial of ``n − 1``.

So, ``3!`` is ``3`` times ``2!``, which is ``2`` times ``1!``, which is ``1`` times ``0!``. Putting it all together, ``3!`` equals ``3`` times ``2`` times ``1`` times ``1``, which is ``6``.

If you can write a recursive definition of something, you can write a Julia program to evaluate it. The first step is to decide what the parameters should be. In this case it should be clear that factorial takes an integer:

```julia
function fact(n) end
```

If the argument happens to be `0` all we have to do is return `1`:

```julia
function fact(n)
	if n == 0
		1
	end
end
```

Otherwise, and this is the interesting part, we have to make a recursive call to find the factorial of `n-1` and then multiply it by `n`:

```julia
function fact(n)
	if n==0
		return 1
	else
		recurse = fact(n-1)
		result = n * recurse
		return result
	end
end
```

The flow of execution for this program is similar to the flow of `countdown` in “Recursion”. If we call fact with the value `3`:

* Since `3` is not `0`, we take the second branch and calculate the factorial of `n-1`...
  * Since `2` is not `0`, we take the second branch and calculate the factorial of `n-1`...
    * Since `1` is not `0`, we take the second branch and calculate the factorial of `n-1`...
      * Since `0` equals `0`, we take the first branch and return `1` without making an more recursive calls.
    * The return value, `1`, is multiplied by `n`, which is `1`, and the result is returned.
  * The return value, `1`, is multiplied by `n`, which is `2`, and the result is returned.
* The return value `2` is multiplied by `n`, which is `3`, and the result, `6`, becomes the return value of the function call that started the whole process.

Figure 6-1 shows what the stack diagram looks like for this sequence of function calls.
"""

# ╔═╡ 2e8091da-7618-11eb-3181-a5801b7aaacb
Drawing(width=720, height=200) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	text(x=140, y=30, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("Main")
	end
    rect(x=170, y=10, width=400, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=140, y=70, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("fact")
	end
    rect(x=170, y=50, width=400, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=190, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=210, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=240, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("3")
	end
	text(x=300, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("recurse")
	end
	text(x=370, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=400, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	text(x=460, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("result")
	end
	text(x=520, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=550, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("6")
	end
	text(x=140, y=110, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("fact")
	end
    rect(x=170, y=90, width=400, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=190, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=210, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=240, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	text(x=300, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("recurse")
	end
	text(x=370, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=400, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=460, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("result")
	end
	text(x=520, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=550, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	text(x=140, y=150, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("fact")
	end
    rect(x=170, y=130, width=400, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=190, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=210, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=240, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=300, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("recurse")
	end
	text(x=370, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=400, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=460, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("result")
	end
	text(x=520, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=550, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=140, y=190, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("fact")
	end
    rect(x=170, y=170, width=400, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=190, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=210, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=240, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("0")
	end
	path(d="M 570 50 C 595 50 595 40 580 40", stroke="black", fill="transparent", marker_end="url(#arrow)")
	text(x=600, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("6")
	end
	path(d="M 570 90 C 595 90 595 80 580 80", stroke="black", fill="transparent", marker_end="url(#arrow)")
	text(x=600, y=90, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	path(d="M 570 130 C 595 130 595 120 580 120", stroke="black", fill="transparent", marker_end="url(#arrow)")
	text(x=600, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	path(d="M 570 170 C 595 170 595 160 580 160", stroke="black", fill="transparent", marker_end="url(#arrow)")
	text(x=600, y=170, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
end

# ╔═╡ 336130c4-761d-11eb-3574-3b57c12a393a
md"*Figure 6-1. Stack diagram.*"

# ╔═╡ 885ca14c-761f-11eb-072c-c9a4a2421f8f
md"""The return values are shown being passed back up the stack. In each frame, the return value is the value of result, which is the product of n and recurse.

In the last frame, the local variables recurse and result do not exist, because the branch that creates them does not run.

!!! tip
    Julia provides the function factorial to calculate the factorial of an integer number."""

# ╔═╡ a4c4b8a8-761f-11eb-345c-572d653d9b0b
md"""## Leap of Faith

Following the flow of execution is one way to read programs, but it can quickly become overwhelming. An alternative is what I call the “leap of faith.” When you come to a function call, instead of following the flow of execution, you *assume* that the function works correctly and returns the right result.

In fact, you are already practicing this leap of faith when you use built-in functions. When you call `cos` or `exp`, you don’t examine the bodies of those functions. You just assume that they work because the people who wrote the built-in functions were good programmers.

The same is true when you call one of your own functions. For example, in “Boolean Functions”, we wrote a function called `isdivisible` that determines whether one number is divisible by another. Once we have convinced ourselves that this function is correct—by examining the code and testing it—we can use the function without looking at the body again.

The same is true of recursive programs. In our example, when you get to the recursive call, instead of following the flow of execution, you should assume that the recursive call works (returns the correct result) and then ask yourself, “Assuming that I can find the factorial of `n-1`, can I compute the factorial of n?” It is clear that you can, by multiplying by `n`.

Of course, it’s a bit strange to assume that the function works correctly when you haven’t finished writing it, but that’s why it’s called a leap of faith!"""

# ╔═╡ d9ce6ae4-761f-11eb-346a-df17c8a60871
md"""## One More Example

After factorial, the most common example of a recursively defined mathematical function is Fibonacci:

```math
\textit{fib}(n)=\begin{cases}
0&\textrm{if }n=0\\
1&\textrm{if }n=1\\
\textit{fib}(n-1)+\textit{fib}(n-2)&\textrm{if }n>1
\end{cases}
```

Translated into Julia, it looks like this:

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

If you try to follow the flow of execution here, even for fairly small values of `n`, your head explodes. But if you take the leap of faith and you assume that the two recursive calls work correctly, then it is clear that you get the right result by adding them together.
"""

# ╔═╡ 8e7878f2-7620-11eb-058a-7346748f29b0
md"""## Checking Types

What happens if we call `fact` and give it `1.5` as an argument?

```julia
julia> fact(1.5)
ERROR: StackOverflowError:
Stacktrace:
 [1] fact(::Float64) at ./REPL[15]:0
 [2] fact(::Float64) at ./REPL[15]:5 (repeats 79983 times)
```

It looks like an infinite recursion. How can that be? The function has a base case— when `n == 0`. But if `n` is not an integer, we can *miss* the base case and recurse forever.

In the first recursive call, the value of `n` is `0.5`. In the next, it is `-0.5`. From there, it gets smaller (more negative), but it will never be `0`.

We have two choices. We can try to generalize the factorial function to work with floating-point numbers, or we can make `fact` check the type of its argument. The first option is called the gamma function, and it’s a little beyond the scope of this book. So we’ll go for the second.

We can use the built-in operator `isa` to verify the type of the argument. While we’re at it, we can also make sure the argument is positive:

```julia
function fact(n)
	if !(n isa Int64)
		error("Factorial is only defined for integers.")
	else if n < 0
		error("Factorial is not defined for negative integers.")
	elseif n == 0
		return 1
	else
		return n * fact(n-1)
	end
end
```

The first base case handles nonintegers; the second handles negative integers. In both cases, the program prints an error message and returns nothing to indicate that something went wrong:

```julia
julia> fact("fred")
ERROR: Factorial is only defined for integers.
julia> fact(-2)
ERROR: Factorial is not defined for negative integers.
```

If we get past both checks, we know that `n` is positive or `0`, so we can prove that the recursion terminates.

This program demonstrates a pattern sometimes called a guardian. The first two con‐ ditionals act as guardians, protecting the code that follows from values that might cause an error. The guardians make it possible to prove the correctness of the code.

In “Catching Exceptions” we will see a more flexible alternative to printing an error message: raising an exception."""

# ╔═╡ 431b4dc0-7621-11eb-1eb8-59d8c7c61a48
md"""## Debugging

Breaking a large program into smaller functions creates natural checkpoints for debugging. If a function is not working, there are three possibilities to consider:

* There is something wrong with the arguments the function is getting; a precondition is violated.
* There is something wrong with the function; a postcondition is violated.
* There is something wrong with the return value or the way it is being used.

To rule out the first possibility, you can add a print statement at the beginning of the function and display the values of the parameters (and maybe their types). Or you can write code that checks the preconditions explicitly.

If the parameters look good, add a print statement before each return statement and display the return value. If possible, check the result by hand. Consider calling the function with values that make it easy to check the result (as in “Incremental Development”).

If the function seems to be working, look at the function call to make sure the return value is being used correctly (or used at all!).

Adding print statements at the beginning and end of a function can help make the flow of execution more visible. For example, here is a version of `fact` with print statements:

```julia
function fact(n)
	space = " " ^ (4 * n)
	println(space, "factorial ", n)
	if n == 0
		println(space, "returning 1")
		return 1
	else
		recurse = fact(n-1)
		result = n * recurse
		println(space, "returning ", result)
		return result
	end
end
```

`space` is a string of space characters that controls the indentation of the output:

```julia
julia> fact(4)
                factorial 4
            factorial 3
        factorial 2
    factorial 1
factorial 0
returning 1
    returning 1
        returning 2
            returning 6
                returning 24
24
```

If you are confused about the flow of execution, this kind of output can be helpful. It takes some time to develop effective scaffolding, but a little bit of scaffolding can save a lot of debugging."""

# ╔═╡ b23efabc-7621-11eb-1d20-9959e0cda861
md"""## Glossary

*temporary variable*:
A variable used to store an intermediate value in a complex calculation.

*dead code*:
Part of a program that can never run, often because it appears after a return statement.

*incremental development*:
A program development plan intended to avoid debugging by adding and testing only a small amount of code at a time.

*scaffolding*:
Code that is used during program development but is not part of the final version.

*guardian*:
A programming pattern that uses a conditional statement to check for and handle circumstances that might cause an error."""

# ╔═╡ e30a3db6-7621-11eb-0dbd-b962759a4c5a
md"""## Exercises

#### Exercise 6-4

Draw a stack diagram for the following program. What does the program print?

```julia
function b(z)
	prod = a(z, z)
    println(z, " ", prod)
	prod
end

function a(x, y)
	x=x+1
	x*y
end

function c(x, y, z)
	total=x+y+z
	square = b(total)^2
	square
end

x=1
y=x+1
println(c(x, y+3, x+y))
```
"""

# ╔═╡ 1269fe6e-7622-11eb-3da2-df0623082981
md"""#### Exercise 6-5

The Ackermann function, ``A(m, n)``, is defined as:

```math
A(m,n)=\begin{cases}
n+1&\textrm{if }m=0\\
A(m-1,1)&\textrm{if }m>0\textrm{ and }n=0\\
A(m-1,A(m,n-1))&\textrm{if }m>0\textrm{ and }n>0
\end{cases}
```

Write a function named `ack` that evaluates the Ackermann function. Use your function to evaluate `ack(3, 4)`, which should be `125`. What happens for larger values of `m` and `n`?
"""

# ╔═╡ 80d44aee-7622-11eb-0555-e3f485a65fe9
md"""#### Exercise 6-6

A palindrome is a word that is spelled the same backward and forward, like “noon” or “redivider.” Recursively, a word is a palindrome if the first and last letters are the same and the middle is a palindrome.

The following are functions that take a string argument and return the first, last, and middle letters:

```julia
function first(word)
	first = firstindex(word)
	word[first]
end

function last(word)
	last = lastindex(word)
	word[last]
end

function middle(word)
	first = firstindex(word)
	last = lastindex(word)
	word[nextind(word, first) : prevind(word, last)]
end
```

We’ll see how they work in Chapter 8.

1. Test these functions out. What happens if you call `middle` with a string with two letters? One letter? What about the empty string, which is written "" and contains no letters?
2. Write a function called `ispalindrome` that takes a string argument and returns `true` if it is a palindrome and `false` otherwise. Remember that you can use the built-in function length to check the length of a string."""

# ╔═╡ cc009a36-7622-11eb-1c86-1fbb98e9fdb3
md"""#### Exercise 6-7

A number, ``a``, is a power of ``b`` if it isdivisible by ``b`` and ``\frac{a}{b}`` is a power of ``b``. Write a function called `ispower` that takes parameters `a` and `b` and returns `true` if `a` is a power of `b`.

!!! tip
    You will have to think about the base case.
"""

# ╔═╡ 260955f6-7623-11eb-332f-dbba32b4a3d9
md"""#### Exercise 6-8

The greatest common divisor (GCD) of ``a`` and ``b`` is the largest number that divides both of them with no remainder.

One way to find the GCD of two numbers is based on the observation that if ``r`` is the remainder when ``a`` is divided by ``b``, then `gcd(a, b) = gcd(b, r)`. As a base case, we can use `gcd(a, 0) = a`.

Write a function called `gcd` that takes parameters `a` and `b` and returns their greatest common divisor.

*Credit*: This exercise is based on an example from Hal Abelson and Gerald Jay Sussman’s *Structure and Interpretation of Computer Programs* (MIT Press)."""

# ╔═╡ Cell order:
# ╟─32be2dac-7618-11eb-2ea5-9f165c3c91e1
# ╟─334f2774-7610-11eb-197e-eb4c20530a0c
# ╟─b1a8c332-7610-11eb-1b18-f7bf7f28ea27
# ╟─4219f3ac-7612-11eb-3167-f55e643f9235
# ╟─8b79bba4-7682-11eb-2ea1-0986a9911b40
# ╟─bfecb80e-7612-11eb-0633-035cbcd3e426
# ╟─f5edc4c2-7613-11eb-0765-cba4f0286212
# ╟─47e45c90-7615-11eb-0fe7-358b9d219497
# ╟─a74ddcba-7615-11eb-1ba0-c1d9a4eacfb5
# ╟─3a20604e-7616-11eb-20c6-1fb57dbbb3c8
# ╟─53822bd0-7616-11eb-08df-0b6ef72120a8
# ╟─2e8091da-7618-11eb-3181-a5801b7aaacb
# ╟─336130c4-761d-11eb-3574-3b57c12a393a
# ╟─885ca14c-761f-11eb-072c-c9a4a2421f8f
# ╟─a4c4b8a8-761f-11eb-345c-572d653d9b0b
# ╟─d9ce6ae4-761f-11eb-346a-df17c8a60871
# ╟─8e7878f2-7620-11eb-058a-7346748f29b0
# ╟─431b4dc0-7621-11eb-1eb8-59d8c7c61a48
# ╟─b23efabc-7621-11eb-1d20-9959e0cda861
# ╟─e30a3db6-7621-11eb-0dbd-b962759a4c5a
# ╟─1269fe6e-7622-11eb-3da2-df0623082981
# ╟─80d44aee-7622-11eb-0555-e3f485a65fe9
# ╟─cc009a36-7622-11eb-1c86-1fbb98e9fdb3
# ╟─260955f6-7623-11eb-332f-dbba32b4a3d9


text:
using Documenter, Waveforms

makedocs(modules = [Waveforms], sitename = "Waveforms.jl")
deploydocs(repo = "github.com/Paalon/Waveforms.jl.git")


text:
# compare exponential kernel specialized system with
# exponential kernel in unspecialized system

# TODO (ask claudia to do it? :-P)


text:
export DiffEqObjective, build_loss_objective

struct DiffEqObjective{F,F2} <: Function
  cost_function::F
  cost_function2::F2
end

(f::DiffEqObjective)(x) = f.cost_function(x)
(f::DiffEqObjective)(x,y) = f.cost_function2(x,y)

function build_loss_objective(prob::DiffEqBase.DEProblem,alg,loss,regularization=nothing;
                              prior=nothing,mpg_autodiff = false,
                              verbose_opt = false,verbose_steps = 100,
                              prob_generator = STANDARD_PROB_GENERATOR,
                              autodiff_prototype = mpg_autodiff ? zero(prob.p) : nothing,
                              autodiff_chunk = mpg_autodiff ? ForwardDiff.Chunk(autodiff_prototype) : nothing,
                              kwargs...)
  if verbose_opt
    count = 0 # keep track of # function evaluations
  end
  cost_function = function (p)
    tmp_prob = prob_generator(prob,p)
    if typeof(loss) <: Union{L2Loss,LogLikeLoss}
      sol = solve(tmp_prob,alg;saveat=loss.t,save_everystep=false,dense=false,kwargs...)
    else
      sol = solve(tmp_prob,alg;kwargs...)
    end

    loss_val = loss(sol)

    if prior != nothing
      loss_val += prior_loss(prior,p)
    end
    if regularization != nothing
      loss_val += regularization(p)
    end

    if verbose_opt
      count::Int += 1
      if mod(count,verbose_steps) == 0
        println("Iteration: $count")
        println("Current Cost: $loss_val")
        println("Parameters: $p")
      end
    end
    loss_val
  end

  if mpg_autodiff
    gcfg = ForwardDiff.GradientConfig(cost_function, autodiff_prototype, autodiff_chunk)
    g! = (x, out) -> ForwardDiff.gradient!(out, cost_function, x, gcfg)
  else
    g! = (x, out) -> Calculus.finite_difference!(cost_function,x,out,:central)
  end
  cost_function2 = function (p,grad)
    if length(grad)>0
      g!(p,grad)
    end
    cost_function(p)
  end
  DiffEqObjective(cost_function,cost_function2)
end


text:
@testset "Learner" begin
    learn = test_learner()
    add_cb!(learn,TestCallback())

    run_learner(learn)
end


text:
@testset "Partial MD" begin
    st = raw"""
        \newcommand{\com}{HH}
        \newcommand{\comb}[1]{HH#1HH}
        A list
        * \com and \comb{blah}
        * $f$ is a function
        * a last element
        """

    steps = explore_md_steps(st)
    lxdefs, tokens, braces, blocks, lxcoms = steps[:latex]

    @test length(braces) == 1
    @test F.content(braces[1]) == "blah"

    @test length(blocks) == 1
    @test blocks[1].name == :MATH_A
    @test F.content(blocks[1]) == "f"

    b2insert, = steps[:b2insert]

    inter_md, mblocks = F.form_inter_md(st, b2insert, lxdefs)
    @test inter_md // """
            A list
            *  ##FDINSERT## and  ##FDINSERT##
            *  ##FDINSERT## is a function
            * a last element"""
    inter_html = F.md2html(inter_md)
    @test inter_html // """

A list

##FDINSERT## and  ##FDINSERT##

##FDINSERT## is a function

a last element

        """
end

# index arithmetic over a string is a bit trickier when using all symbols
# we can use `prevind` and `nextind` to make sure it works properly
@testset "Inter Md 2" begin
    st = raw"""
        ~~~
        this⊙ then ⊙ ⊙ and
        ~~~
        finally ⊙⊙𝛴⊙ and
        ~~~
        escape ∀⊙∀
        ~~~
        done
        """

    inter_md, = explore_md_steps(st)[:inter_md]
    @test inter_md // """
                     ##FDINSERT##
                    finally ⊙⊙𝛴⊙ and
                     ##FDINSERT##
                    done
                    """
end

@testset "Latex eqa" begin
    st = raw"""
        a\newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}b@@d .@@
        \eqa{\sin^2(x)+\cos^2(x) &=& 1}
        """

    steps = explore_md_steps(st)
    lxdefs, tokens, braces, blocks, lxcoms = steps[:latex]
    b2insert, = steps[:b2insert]
    inter_md, mblocks = steps[:inter_md]
    @test inter_md // """
                ab

                ##FDINSERT##

                 ##FDINSERT##
                """

    inter_html, = steps[:inter_html]

    @test F.convert_block(b2insert[1], lxdefs) == "
.
"
    @test isapproxstr(F.convert_block(b2insert[2], lxdefs), "\\[\\begin{array}{c} \\sin^2(x)+\\cos^2(x) &=& 1\\end{array}\\]")
    hstring = F.convert_inter_html(inter_html, b2insert, lxdefs)
    @test isapproxstr(hstring, raw"""

ab

.

                    \[\begin{array}{c} \sin^2(x)+\cos^2(x) &=& 1\end{array}\]""")
end

@testset "MD>HTML" begin
    st = raw"""
        text A1 \newcommand{\com}{blah}text A2 \com and
        ~~~
        escape B1
        ~~~
        \newcommand{\comb}[ 1]{\mathrm{#1}} text C1 $\comb{b}$ text C2
        \newcommand{\comc}[ 2]{part1:#1 and part2:#2} then \comc{AA}{BB}.
        """

    steps = explore_md_steps(st)
    lxdefs, tokens, braces, blocks, lxcoms = steps[:latex]
    b2insert, = steps[:b2insert]
    inter_md, mblocks = steps[:inter_md]
    inter_html, = steps[:inter_html]

    @test isapproxstr(inter_md, """
            text A1 text A2  ##FDINSERT## and
             ##FDINSERT##
             text C1  ##FDINSERT## text C2
             then  ##FDINSERT##.""")

    @test isapproxstr(inter_html, """
text A1 text A2  ##FDINSERT## and  ##FDINSERT##  text C1  ##FDINSERT## text C2  then  ##FDINSERT##.
""")

    hstring = F.convert_inter_html(inter_html, b2insert, lxdefs)
    @test isapproxstr(hstring, raw"""

text A1 text A2 blah and
                    escape B1
                    text C1 \(\mathrm{ b}\) text C2  then part1: AA and part2: BB.
""")
end

@testset "headers" begin
    set_curpath("index.md")
    h = """
        # Title
        and then
        ## Subtitle cool!
        done
        """ |> seval
    @test h // """

Title

and then

Subtitle cool!

done
"""
end


text:
### Guru Interfaces

# abbreviated name for the type for C pointer used in plan struct
finufft_plan_c = Ptr{Cvoid}

mutable struct finufft_plan{T}
    type       :: Cint
    ntrans     :: Cint
    dim        :: Cint
    ms         :: BIGINT
    mt         :: BIGINT
    mu         :: BIGINT
    nj         :: BIGINT
    nk         :: BIGINT
    plan_ptr   :: finufft_plan_c
end

"""
    finufft_makeplan(type::Integer,
                          n_modes_or_dim::Union{Array{Int64},Integer},
                          iflag::Integer,
                          ntrans::Integer,
                          eps::Real;
                          dtype=Float64,
                          kwargs...)

Creates a `finufft_plan` object for the guru interface to FINUFFT, of
 type 1, 2 or 3, and with given numbers of Fourier modes (unless type 3).

 # Inputs
  - `type`            transform type: 1, 2, or 3
  - `n_modes_or_dim`  if `type` is 1 or 2, the number of Fourier modes in each
                    dimension: `ms` in 1D, `[ms mt]` in 2D, or `[ms mt mu]` in 3D.
                    Its length thus sets the dimension, which must be 1, 2 or 3.
                    If `type` is 3, in contrast, its *value* fixes the dimension.
  -  `iflag`   if >=0, uses + sign in exponential, otherwise - sign.
  -  `ntrans`          number of transforms to compute simultaneously
  -  `eps`     relative precision requested (generally between 1e-15 and 1e-1),
            real, need not match type of `dtype`
  - `dtype`           `Float32` or `Float64`, plan for single precision or double precision
  -  `kwargs`  (optional): for options, see `nufft_opts` and
    https://finufft.readthedocs.io/en/latest/opts.html
 # Returns
  - finufft_plan struct
# Examples
```julia-repl
julia> p = finufft_makeplan(2,10,+1,1,1e-6);
```
creates a plan for a 1D type 2 Float64 transform with 10 Fourier modes
and tolerance 1e-6.
```julia-repl
julia> p = finufft_makeplan(1,[10 20],+1,1,1e-6);
```
creates a plan for a 2D type 1 Float64 transform with 10*20 Fourier modes.
```julia-repl
julia> p = finufft_makeplan(3,2,+1,1,1e-4,dtype=Float32,nthreads=4);
```
creates a plan for a 2D type 3 Float32 transform with tolerance 1e-4,
to use 4 threads.
"""
function finufft_makeplan(type::Integer,
                          n_modes_or_dim::Union{Array{BIGINT},Integer},
                          iflag::Integer,
                          ntrans::Integer,
                          eps::Real;
                          dtype=Float64,
                          kwargs...)
# see https://stackoverflow.com/questions/40140699/the-proper-way-to-declare-c-void-pointers-in-julia for how to declare c-void pointers in julia
#   one can also use Array/Vector for cvoid pointer, Array and Ref both work
#   plan_p = Array{finufft_plan_c}(undef,1)
    plan_p = Ref{finufft_plan_c}()

    opts = finufft_default_opts(dtype)
    setkwopts!(opts;kwargs...)

    n_modes = ones(BIGINT,3)
    if type==3
        @assert ndims(n_modes_or_dim) == 0
        dim = n_modes_or_dim
    else
        @assert length(n_modes_or_dim)<=3 && length(n_modes_or_dim)>=1
        dim = length(n_modes_or_dim)
        n_modes[1:dim] .= n_modes_or_dim
    end

    if dtype==Float64
        tol = Float64(eps)
        ret = ccall( (:finufft_makeplan, libfinufft),
                      Cint,
                      (Cint,
                       Cint,
                       Ref{BIGINT},
                       Cint,
                       Cint,
                       Cdouble,
                       Ptr{finufft_plan_c},
                       Ref{nufft_opts}),
                      type,dim,n_modes,iflag,ntrans,tol,plan_p,opts
                      )

    else
        tol = Float32(eps)
        ret = ccall( (:finufftf_makeplan, libfinufft),
                      Cint,
                      (Cint,
                       Cint,
                       Ref{BIGINT},
                       Cint,
                       Cint,
                       Cfloat,
                       Ptr{finufft_plan_c},
                       Ref{nufft_opts}),
                      type,dim,n_modes,iflag,ntrans,tol,plan_p,opts
                      )
    end
    check_ret(ret)

    ms = n_modes[1]
    mt = n_modes[2]
    mu = n_modes[3]
    plan = finufft_plan{dtype}(type,ntrans,dim,ms,mt,mu,0,0,plan_p[])
    return plan
end

"""
    finufft_setpts(plan, xj [, yj[, zj[, s[, t[, u]]]]])

Input nonuniform points for general FINUFFT transform(s).

 Given an already-planned `finufft_plan`, this reads in nonuniform point
 coordinate arrays `xj` (and `yj` if 2D or 3D, and `zj` if 3D), and additionally in the type 3 case, nonuniform
 frequency target coordinate arrays `s` (and `t` if 2D or 3D, and `u` if 3D).
Empty arrays may be passed in the case of
 unused dimensions. For all types, sorting is done to internally store a
 reindexing of points, and for type 3 the spreading and FFTs are planned.
 These nonuniform points may then be used for multiple transforms.

 # Inputs
  -  `plan`   a `finufft_plan` object for one/many general nonuniform FFTs
  -  `xj` Array{Float32} or Array{Float64}, vector of x-coords of all nonuniform points
  -  `yj`     empty (if dim<2), or vector of y-coords of all nonuniform points
-  `zj`     empty (if dim<3), or vector of z-coords of all nonuniform points
 -   `s`      vector of x-coords of all nonuniform frequency targets
  -  `t`      empty (if dim<2), or vector of y-coords of all frequency targets
  -  `u`      empty (if dim<3), or vector of z-coords of all frequency targets
"""
function finufft_setpts(plan::finufft_plan{T},
                        xj::Array{T},
                        yj::Array{T}=T[],
                        zj::Array{T}=T[],
                        s::Array{T}=T[],
                        t::Array{T}=T[],
                        u::Array{T}=T[]) where T <: finufftReal

    (M, N) = valid_setpts(plan.type, plan.dim, xj, yj, zj, s, t, u)

    plan.nj = M
    plan.nk = N

    if T==Float64
        T_real = Array{Float64}
        ret = ccall( (:finufft_setpts, libfinufft),
                     Cint,
                     (finufft_plan_c,
                      BIGINT,
                      Ref{Cdouble},
                      Ref{Cdouble},
                      Ref{Cdouble},
                      BIGINT,
                      Ref{Cdouble},
                      Ref{Cdouble},
                      Ref{Cdouble}),
                     plan.plan_ptr,M,T_real(xj),T_real(yj),T_real(zj),N,T_real(s),T_real(t),T_real(u)
                     )
    else
        T_real = Array{Float32}
        ret = ccall( (:finufftf_setpts, libfinufft),
                     Cint,
                     (finufft_plan_c,
                      BIGINT,
                      Ref{Cfloat},
                      Ref{Cfloat},
                      Ref{Cfloat},
                      BIGINT,
                      Ref{Cfloat},
                      Ref{Cfloat},
                      Ref{Cfloat}),
                     plan.plan_ptr,M,T_real(xj),T_real(yj),T_real(zj),N,T_real(s),T_real(t),T_real(u)
                     )
    end

    check_ret(ret)
    return ret
end

"""
    output::Array{Complex{T}} = finufft_exec(plan::finufft_plan{T},
                      input::Array{Complex{T}}) where T <: finufftReal

Execute single or many-vector FINUFFT transforms in a plan.

  output = finufft_exec(plan, input)

  For `plan` a previously created `finufft_plan` object also containing all
  needed nonuniform point coordinates, do a single (or if `ntrans>1` in the
  plan stage, multiple) NUFFT transform(s), with the strengths or Fourier
  coefficient inputs vector(s) from `input`. The result of the transform(s)
  is returned as a (possibly multidimensional) array.

 # Inputs
    - `plan`     `finufft_plan` object, already planned and containing
    nonuniform points.
    - `input`    strengths (types 1 or 3) or Fourier coefficients (type 2)
              vector, matrix, or array of appropriate size. For type 1 and 3,
              this is either a length-M vector (where M is the length of `xj`),
              or an `(M,ntrans)` matrix when `ntrans>1`. For type 2, in 1D this is size `(ms,)`, in 2D size `(ms,mt)`, or in 3D size `(ms,mt,mu)`, or
              each of these with an extra last dimension `ntrans` if `ntrans>1`.
# Output
     `output`   vector of output strengths at targets (types 2 or 3), or array
              of Fourier coefficients (type 1), or, if `ntrans>1`, a stack of
              such vectors or arrays, of appropriate size.
              Specifically, if `ntrans=1`, for type 1, in 1D
              this is size `(ms,)`, in 2D size
              `(ms,mt)`, or in 3D size `(ms,mt,mu)`; for types 2 and 3
              it is a column vector of length `M` (the length of `xj` in type 2),
              or `nk` (the length of `s` in type 3). If `ntrans>1` it is a stack
              of such objects, ie, it has an extra last dimension `ntrans`.
"""
function finufft_exec(plan::finufft_plan{T},
                      input::Array{Complex{T}}) where T <: finufftReal
    ret = 0
    type = plan.type
    ntrans = plan.ntrans
    dim = plan.dim
    n_modes = Array{BIGINT}(undef,3)
    n_modes[1] = plan.ms
    n_modes[2] = plan.mt
    n_modes[3] = plan.mu
    if type==1
        if dim==1
            output = Array{Complex{T}}(undef,n_modes[1],ntrans)
        elseif dim==2
            output = Array{Complex{T}}(undef,n_modes[1],n_modes[2],ntrans)
        elseif dim==3
            output = Array{Complex{T}}(undef,n_modes[1],n_modes[2],n_modes[3],ntrans)
        else
            ret = ERR_DIM_NOTVALID
        end
    elseif type==2
        nj = plan.nj
        output = Array{Complex{T}}(undef,nj,ntrans)
    elseif type==3
        nk = plan.nk
        output = Array{Complex{T}}(undef,nk,ntrans)
    else
        ret = ERR_TYPE_NOTVALID
    end
    check_ret(ret)
    finufft_exec!(plan,input,output)
    return output
end

"""
    status = finufft_destroy(plan::finufft_plan{T}) where T <: finufftReal

This destroys a FINUFFT plan object: it
deallocates all stored FFTW plans, nonuniform point sorting arrays,
kernel Fourier transforms arrays, and any other allocations, and
nulls the plan pointer.

An integer status code is returned, that is 0 if successful.
If one attempts to destroy an already-destroyed plan, 1 is returned
(see FINUFFT documentation for finufft_destroy).
"""
function finufft_destroy(plan::finufft_plan{T}) where T <: finufftReal
    if plan.plan_ptr!=C_NULL     # this test should not be needed since
        # ccall should handle C_NULL returning 1, but ...that failed one CI test
        if T==Float64
            ret = ccall( (:finufft_destroy, libfinufft),
                         Cint,
                         (finufft_plan_c,),
                         plan.plan_ptr
                         )
        else
            ret = ccall( (:finufftf_destroy, libfinufft),
                         Cint,
                         (finufft_plan_c,),
                         plan.plan_ptr
                         )
        end
        plan.plan_ptr = C_NULL       # signifies destroyed
        return ret                   # pass out status of the ccall
    else
        return 1                     # signifies plan was already destroyed
    end
end

"""
    finufft_exec!(plan::finufft_plan{T},
                      input::Array{Complex{T}},
                      output::Array{Complex{T}}) where T <: finufftReal

Execute single or many-vector FINUFFT transforms in a plan, with output written
to preallocated array. See `finufft_exec` for arguments.
"""
function finufft_exec!(plan::finufft_plan{T},
                      input::Array{Complex{T}},
                      output::Array{Complex{T}}) where T <: finufftReal
    type = plan.type
    ntrans = plan.ntrans
    dim = plan.dim
    n_modes = Array{BIGINT}(undef,3)
    n_modes[1] = plan.ms
    n_modes[2] = plan.mt
    n_modes[3] = plan.mu
    if type==1
        if dim==1
            if ntrans==1
                @assert size(output)==(n_modes[1],) || size(output)==(n_modes[1],ntrans)
            else
                @assert size(output)==(n_modes[1],ntrans)
            end
        elseif dim==2
            if ntrans==1
                @assert size(output)==(n_modes[1],n_modes[2]) || size(output)==(n_modes[1],n_modes[2],ntrans)
            else
                @assert size(output)==(n_modes[1],n_modes[2],ntrans)
            end
        elseif dim==3
            if ntrans==1
                @assert size(output)==(n_modes[1],n_modes[2],n_modes[3]) || size(output)==(n_modes[1],n_modes[2],n_modes[3],ntrans)
            else
                @assert size(output)==(n_modes[1],n_modes[2],n_modes[3],ntrans)
            end
        else
            ret = ERR_DIM_NOTVALID
            check_ret(ret)
        end
        if T==Float64
            ret = ccall( (:finufft_execute, libfinufft),
                         Cint,
                         (finufft_plan_c,
                          Ref{ComplexF64},
                          Ref{ComplexF64}),
                         plan.plan_ptr,input,output
                         )
        else
            ret = ccall( (:finufftf_execute, libfinufft),
                         Cint,
                         (finufft_plan_c,
                          Ref{ComplexF32},
                          Ref{ComplexF32}),
                         plan.plan_ptr,input,output
                         )
        end
    elseif type==2
        nj = plan.nj
        if ntrans==1
            @assert size(output)==(nj,ntrans) || size(output)==(nj,)
        else
            @assert size(output)==(nj,ntrans)
        end
        if T==Float64
            ret = ccall( (:finufft_execute, libfinufft),
                         Cint,
                         (finufft_plan_c,
                          Ref{ComplexF64},
                          Ref{ComplexF64}),
                         plan.plan_ptr,output,input
                         )
        else
            ret = ccall( (:finufftf_execute, libfinufft),
                         Cint,
                         (finufft_plan_c,
                          Ref{ComplexF32},
                          Ref{ComplexF32}),
                         plan.plan_ptr,output,input
                         )
        end
    elseif type==3
        nk = plan.nk
        if ntrans==1
            @assert size(output)==(nk,ntrans) || size(output)==(nk,)
        else
            @assert size(output)==(nk,ntrans)
        end
        if T==Float64
            ret = ccall( (:finufft_execute, libfinufft),
                         Cint,
                         (finufft_plan_c,
                          Ref{ComplexF64},
                          Ref{ComplexF64}),
                         plan.plan_ptr,input,output
                         )
        else
            ret = ccall( (:finufftf_execute, libfinufft),
                         Cint,
                         (finufft_plan_c,
                          Ref{ComplexF32},
                          Ref{ComplexF32}),
                         plan.plan_ptr,input,output
                         )
        end
    else
        ret = ERR_TYPE_NOTVALID
    end
    check_ret(ret)
end