text:
"""
	makecpt(cmd0::String="", arg1=nothing; kwargs...)

Make static color palette tables (CPTs).

Full option list at [`makecpt`]($(GMTdoc)makecpt.html)

- **A** | **alpha** | **transparency** :: [Type => Str]

    Sets a constant level of transparency (0-100) for all color slices.
    ($(GMTdoc)makecpt.html#a)
- $(GMT.opt_C)
- **D** | **bg** | **background** :: [Type => Str | []]			`Arg = [i|o]`

    Select the back- and foreground colors to match the colors for lowest and highest
    z-values in the output CPT.
    ($(GMTdoc)makecpt.html#d)
- **E** | **nlevels** :: [Type => Int | []]		`Arg = [nlevels]`

    Implies reading data table(s) from file or arrays. We use the last data column to
    determine the data range
    ($(GMTdoc)makecpt.html#e)
- **F** | **color_model** :: [Type => Str | []]		`Arg = [R|r|h|c][+c]]`

    Force output CPT to written with r/g/b codes, gray-scale values or color name.
    ($(GMTdoc)makecpt.html#f)
- **G** | **truncate** :: [Type => Str]              `Arg = zlo/zhi`

    Truncate the incoming CPT so that the lowest and highest z-levels are to zlo and zhi.
    ($(GMTdoc)makecpt.html#g)
- **H** | **save** :: [Type => Bool]

    Modern mode only: Write the CPT to disk [Default saves the CPT as the session current CPT].
    Required for scripts used to make animations via movie where we must pass named CPT files.
    ($(GMTdoc)makecpt.html#h)
- **I** | **inverse** | **reverse** :: [Type => Str]	`Arg = [c][z]`

    Reverse the sense of color progression in the master CPT.
    ($(GMTdoc)makecpt.html#i)
- **M** | **overrule_bg** :: [Type => Bool]

    Overrule background, foreground, and NaN colors specified in the master CPT with the values of
    the parameters COLOR_BACKGROUND, COLOR_FOREGROUND, and COLOR_NAN.
    ($(GMTdoc)makecpt.html#m)
- **N** | **no_bg** | **nobg** :: [Type => Bool]

    Do not write out the background, foreground, and NaN-color fields.
    ($(GMTdoc)makecpt.html#n)
- **Q** | **log** :: [Type => Bool | Str]			`Arg = [i|o]`

    Selects a logarithmic interpolation scheme [Default is linear].
    ($(GMTdoc)makecpt.html#q)
- **S** | **auto** :: [Type => Bool | Str]			`Arg = [mode]`

    Determine a suitable range for the -T option from the input table(s) (or stdin).
    ($(GMTdoc)makecpt.html#s)
- **T** | **range** :: [Type => Str]			`Arg = [min/max/inc[+b|l|n]|file|list]`

    Defines the range of the new CPT by giving the lowest and highest z-value and interval.
    ($(GMTdoc)makecpt.html#t)
- **W** | **wrap** | **categorical** :: [Type => Bool | Str | []]      `Arg = [w]`

    Do not interpolate the input color table but pick the output colors starting at the
    beginning of the color table, until colors for all intervals are assigned.
    ($(GMTdoc)makecpt.html#w)
- **Z** | **continuous** :: [Type => Bool]

    Creates a continuous CPT [Default is discontinuous, i.e., constant colors for each interval].
    ($(GMTdoc)makecpt.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_h)
- $(GMT.opt_i)
"""
function makecpt(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("makecpt", cmd0, arg1)	# Monolithic mode

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:V_params])

    # If file name sent in, read it and compute a tight -R if this was not provided
    cmd, arg1, = read_data(d, cmd0, cmd, arg1, " ")
	cmd, arg1, = add_opt_cpt(d, cmd, CPTaliases, 'C', 0, arg1)
	cmd, Tvec = helper_cpt(d, cmd)
	cmd = parse_E_mkcpt(d, [:E :nlevels], cmd, arg1)

	cmd = "makecpt " * cmd
	(dbg_print_cmd(d, cmd) !== nothing) && return cmd
	r = gmt(cmd, arg1, !isempty(Tvec) ? Tvec : nothing)
	current_cpt[1] = (r !== nothing) ? r : GMTcpt()
    return r
end

# ---------------------------------------------------------------------------------------------------
function parse_E_mkcpt(d::Dict, symbs::Array{<:Symbol}, cmd::String, arg1)
	(show_kwargs[1]) && return print_kwarg_opts(symbs, "Number")
	if ((val = find_in_dict(d, symbs)[1]) !== nothing)
		(arg1 === nothing) && error("E option requires that a data table is provided as well")
		cmd *= " -E" * arg2str(val)
	end
	return cmd
end

# ---------------------------------------------------------------------------------------------------
function helper_cpt(d::Dict, cmd::String)
	# Common to both make & grd cpt
	cmd = parse_these_opts(cmd, d, [[:A :alpha :transparency], [:D :bg :background], [:F :color_model], [:G :truncate],
	                                [:I :inverse :reverse], [:M :overrule_bg], [:N :no_bg :nobg], [:Q :log], [:S :auto], [:W :wrap :categorical], [:Z :continuous]])
	cmd, Tvec = parse_opt_range(d, cmd, "T")
	if ((val = find_in_dict(d, [:cptname :cmapname])[1]) !== nothing)
		(IamModern[1]) && (cmd *= " -H")
		cmd *=  " > " * string(val)
	elseif (IamModern[1])  cmd *= " -H"
	end
	return cmd, Tvec
end

# -------------------------------------------------------------------------------------------
function parse_opt_range(d::Dict, cmd::String, opt::String="")::Tuple{String, Vector{Float64}}
	symbs = [:T :range :inc :bin]
	(show_kwargs[1]) && return print_kwarg_opts(symbs, "Tuple | Array | String | Number")	# Just print the options
	Tvec = Vector{Float64}()
	if ((val = find_in_dict(d, symbs)[1]) !== nothing)
		if (isa(val, Tuple))
			n = length(val)
			out = arg2str(val[1:min(n,3)])
			if (n > 3)
				for k = 4:n			# N should be at most = 5 (e.g. +n+b)
					_opt = string(val[k])
					if     (startswith(_opt, "sli") || startswith(_opt, "num"))   out *= "+n"
					elseif (startswith(_opt, "log2"))  out *= "+b"
					elseif (startswith(_opt, "log1"))  out *= "+l"
					else   @warn("Unkown option \"$_opt\" in range option")
					end
				end
			end
		elseif (isa(val, Vector{<:Real}) || isa(val, Matrix{<:Real}))
			out = arg2str(val,',')
			if (length(val) == 1)  out *= ","  end
			if (length(out) > 1023)  Tvec, out = vec(val), ""  end		# Should be temporary till all that use -T catch up with Vec
		else
			out = arg2str(val)		# Everything fits here if given as a string
		end
		if (opt != "")  out = " -" * opt  * out end
		cmd *= out
	end
	return cmd, Tvec
end

# ---------------------------------------------------------------------------------------------------
# Version to use with the -E option
makecpt(arg1; kw...) = makecpt("", arg1; kw...)


text:
export average, decimation, optimalDecimation,
       acqStart, acqStop, acqReset, bufferSize

### Sampling Rate and Decimation ###

decimation(rp::RedPitaya) = query(rp, "ACQ:DEC?", Int)

decimation(rp::RedPitaya, dec::Integer) = send(rp,"ACQ:DEC $dec")

function optimalDecimation(rp::RedPitaya, freq)
  if freq > 100e3#2e6
    return 1
  elseif freq > 2e3#50e3
    return 8
  elseif freq > 1e3
    return 64
  elseif freq > 0.5e3
    return 1024
  else
    return 8192
  end
end

function average(rp::RedPitaya)
  a = query(rp,"ACQ:AVG?")
  return a == "ON" ? true : false
end

function average(rp::RedPitaya, avg::Bool)
  avgStr = avg ? "ON" : "OFF"
  cmd ="ACQ:AVG $avgStr"
  send(rp,cmd)
end

### Acquisition Control ###

"""
Start acquisition of RedPitaya
"""
acqStart(rp::RedPitaya) = send(rp,"ACQ:START")

"""
Stop acquisition of RedPitaya
"""
acqStop(rp::RedPitaya) = send(rp,"ACQ:STOP")

"""
Reset acquisition of RedPitaya
"""
acqReset(rp::RedPitaya) = send(rp,"ACQ:RST")

### Acquisition Functions ###

bufferSize(rp::RedPitaya) = query(rp, "ACQ:BUF:SIZE?", Int)

export receiveAnalogSignal
function receiveAnalogSignal(rp::RedPitaya, chan::Integer, from::Integer=-1, to::Integer=-1;
                             dec::Integer=1, delay=0.2, typ="STA", binary=false, raw=false)

  acqReset(rp)
  decimation(rp,dec)
  acqStart(rp)
  sleep(delay) # fill buffers

  return receiveAnalogSignalLowLevel(rp,chan,from,to,typ,binary,raw)
end

function _awg_trigger_delay(rp::RedPitaya, dec::Integer)
  if dec == 1
    return 25
  elseif dec == 8
    return 3
  else
    return 0
  end
end

export receiveAnalogSignalWithTrigger
function receiveAnalogSignalWithTrigger(rp::RedPitaya, chan::Integer, from::Integer=-1, to::Integer=-1;
                             dec::Integer=1, delay=0.2, typ="OLD", trigger="AWG_NE",
                             triggerLevel=0.2, binary=false, raw=false, triggerDelay=0)

  acqReset(rp)
  decimation(rp,dec)
  send(rp,"ACQ:TRIG:LEV $triggerLevel")
  additionalDelay = occursin("AWG", trigger) ? _awg_trigger_delay(rp, dec) : 0
  send(rp,"ACQ:TRIG:DLY $(8192+additionalDelay+triggerDelay)")
  acqStart(rp)
  sleep(delay) # fill buffers

  send(rp,"ACQ:TRIG $trigger")
  if trigger != "NOW" # Ugly
    send(rp,"SOUR1:TRIG:IMM")
  end

  while true
    trig_rsp = query(rp,"ACQ:TRIG:STAT?")
    println(trig_rsp)
    if trig_rsp[1:2] == "TD"
       break
     end
  end

  if chan > 0
    return receiveAnalogSignalLowLevel(rp,chan,from,to,typ,binary,raw)
  else
    u1 = receiveAnalogSignalLowLevel(rp,1,from,to,typ,binary,raw)
    u2 = receiveAnalogSignalLowLevel(rp,2,from,to,typ,binary,raw)
    return u1, u2
  end

end

export receiveAnalogSignalLowLevel
function receiveAnalogSignalLowLevel(rp::RedPitaya, chan::Integer, from::Integer=-1,
                                     to::Integer=-1, typ="OLD", binary::Bool=false,
                                     raw::Bool=false)

  if binary
    send(rp,"ACQ:DATA:FORMAT BIN")
  else
    send(rp,"ACQ:DATA:FORMAT ASCII")
  end

  if raw
    send(rp,"ACQ:DATA:UNITS RAW")
  end

  if from == to == -1
    send(rp,"ACQ:SOUR$(chan):DATA?")
  elseif typ == "STA"
    send(rp,"ACQ:SOUR$(chan):DATA:STA:N? $(from),$(to)")
  else
    send(rp,"ACQ:SOUR$(chan):DATA:$(typ):N? $(to)")
  end

  if binary
    if raw
      uFl = receiveBinaryInt16(rp)
    else
      uFl = receiveBinaryFloat(rp)
    end
  else
    uFl = receiveASCIIArray(rp)
  end

  return uFl
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BisectPointPartition(normal, point)

A method for partitioning spatial data into two half spaces
defined by a `normal` direction and a reference `point`.
"""
struct BisectPointPartition{Dim,T} <: PartitionMethod
  normal::Vec{Dim,T}
  point::Point{Dim,T}

  function BisectPointPartition{Dim,T}(normal, point) where {Dim,T}
    new(normalize(normal), point)
  end
end

BisectPointPartition(normal::Vec{Dim,T}, point::Point{Dim,T}) where {Dim,T} =
  BisectPointPartition{Dim,T}(normal, point)

BisectPointPartition(normal::NTuple{Dim,T}, point::NTuple{Dim,T}) where {Dim,T} =
 BisectPointPartition(Vec(normal), Point(point))

function partition(object, method::BisectPointPartition)
  Dim = embeddim(object)
  T = coordtype(object)

  n = method.normal
  p = method.point

  x = MVector{Dim,T}(undef)

  left  = Vector{Int}()
  right = Vector{Int}()
  for location in 1:nelements(object)
    coordinates!(x, object, location)
    if (Point(x) - p) ⋅ n < zero(T)
      push!(left, location)
    else
      push!(right, location)
    end
  end

  Partition(object, [left,right])
end


text:
abstract type SupervisedLearner end

"""
    measureaccuracy(learner, features, labels[, confusion])

Measure the accuracy of a learner given a set of features and labels.

If the problem is a classification problem, this method returns the percent
of instances which are classified correctly. If it is a regression problem,
it calculates the root mean squared error.

If a confusion matrix is included, and the problem is a classification problem,
the confusion matrix will be updated to show common misclassifications.
"""
function measureaccuracy(learner::SupervisedLearner, features::Matrix, labels::Matrix, confusion::Union{Matrix,Nothing}=nothing)
	if rows(features) != rows(labels)
		error("Expected the features and labels to have the same number of rows")
	elseif columns(labels) != 1
		error("Sorry, this method currently only supports one-dimensional labels")
	elseif rows(features) == 0
		error("Expected at least one row")
	end

	if iscontinuous(labels, 1)
		# The label is continuous, so measure root mean squared error
		sse = 0.
		foreach(features, labels) do feature, target
			prediction = predict(learner, feature)
			delta = target[1] - prediction
			sse += delta^2
		end
		√(sse / rows(features))
	else
		correctcount = 0
		labelvalues = valuecount(labels, 1)
		foreach(features, labels[:,1]) do feature, target
			if target >= labelvalues
				error("The label is out of range")
			end
			prediction = Int(predict(learner, feature))
			if confusion != nothing
				confusion[Int(target) + 1, prediction + 1] += 1
			end
			correctcount += prediction == target
		end
		correctcount / rows(features)
	end
end

"""
    meansquarederror(learner, features, labels)

Calculate the mean squared error of a learner given a set of features and a
corresponding set of labels.
"""
function meansquarederror end

"""
    train(learner, features, labels)

Train a learner given a matrix of features and the corresponding labels
"""
function train end

"""
    predict(learner, instance)

Use a learner to predict the class of a specific instance.
"""
function predict end

max_stars_repo_name: coreywoodfield/mltoolkit-julia

max_issues_repo_name: coreywoodfield/mltoolkit-julia

max_forks_repo_name: coreywoodfield/mltoolkit-julia


text:
# Matlab compability
function ndgrid(x, y)
    nx = length(x)
    ny = length(y)
    vartype = typeof(x[1])
    @assert typeof(y[1])==vartype
    X = Array{vartype}(nx, ny)
    Y = Array{vartype}(nx, ny)
    for i=1:nx
        for j=1:ny
            X[i,j] = x[i]
            Y[i,j] = y[j]
        end
    end
    return X,Y
end


text:
mutable struct TRSinfo
	hard_case::Bool  # Flag indicating if we are in the hard case
	niter::Int # Number of iterations in eigs
	nmul::Int  # Number of multiplication with P in eigs
	λ::Vector  # Lagrange Multiplier(s)

	function TRSinfo(hard_case::Bool, niter::Int, nmul::Int, λ::Vector)
		new(hard_case, niter, nmul, λ)
	end
end

function trs_boundary(P, q::AbstractVector{T}, r::T, C::AbstractMatrix{T}; kwargs...) where {T}
	check_inputs(P, q, r, C)
	return trs_boundary((nev; kw...) -> gen_eigenproblem(P, q, r, C, nev; kw...),
		   (λ, V; kw...) -> pop_solution!(λ, V, P, q, r, C; kw...); kwargs...)
end

function trs_boundary(P, q::AbstractVector{T}, r::T; kwargs...) where {T}
	check_inputs(P, q, r)
	return trs_boundary((nev; kw...) -> eigenproblem(P, q, r, nev; kw...),
		   (λ, V; kw...) -> pop_solution!(λ, V, P, q, r, I; kw...); kwargs...)
end

function trs_boundary(solve_eigenproblem::Function, pop_solution!::Function;
	compute_local=false, tol_hard=2e-7, kwargs...)

	if compute_local
		nev=2  # We will need the two rightmost eigenvalues.
	else
		nev=1  # We will only need the rightmost eigenvalue
	end

	l, V, niter, nmult = solve_eigenproblem(nev; kwargs...)
	X, λ, hard_case = pop_solution!(l, V; tol_hard=tol_hard) # Pop global minimizer(s).
	if !compute_local || hard_case
		return X, TRSinfo(hard_case, niter, nmult, λ)
	else
		X_local, λ_local, _ = pop_solution!(l, V; tol_hard=tol_hard) # Pop local-no-global minimizer.
		return [X X_local], TRSinfo(hard_case, niter, nmult, [λ; λ_local])
	end
end

function check_inputs(P, q::AbstractVector{T}, r::T) where {T}
	@assert(issymmetric(P), "The cost matrix must be symmetric.")
	@assert(eltype(P) == T, "Inconsistent element types.")
	@assert(size(P, 1) == size(P, 2) == length(q), "Inconsistent matrix dimensions.")
end

function check_inputs(P, q::AbstractVector{T}, r::T, C::AbstractMatrix{T}) where {T}
	check_inputs(P, q, r)
	@assert(issymmetric(C), "The norm must be defined by a symmetric positive definite matrix.")
end

function pop_solution!(λ, V, P, q::AbstractVector{T}, r::T, C; tol_hard, direct=false) where {T}
	idx = argmax(real(λ)) 	# Pop rightmost eigenvector
	is_first_pop = all(isfinite.(λ))
	n = length(q)
	if abs(imag(λ[idx])) >= 1e-9 && !is_first_pop # A solution exists only on real eigenvalues
		return zeros(T, n, 0), zeros(T, 0), false
	end
	l = real(λ[idx]);
	complex_v = view(V, :, idx)
	if norm(real(complex_v)) > norm(imag(complex_v))
		v = real(complex_v)
	else
		v = imag(complex_v)  # Sometimes the retuned eigenvector is complex
	end
	v ./= norm(v)
	v1 = view(v, 1:n); v2 = view(v, n+1:2*n)

	if norm(q) < 1e-9
		# Solution the eigenvector of Px = lCx, which is parallel to v2
		x = r*v2/sqrt(v2'*C*v2)
		return [x -x], [l; l], true
	end

	# Extract solution
	norm_v1 = sqrt(dot(v1, C*v1))
	X = zeros(T, n, 0)
	hard_case = false
	if norm_v1 >= tol_hard
		if norm_v1 >= 1e-11
			x = -r*v1/norm_v1
			if sign(q'*v2) < 0
				x .= -x
			end

			residual = (P*x + l*C*x + q)
			# Polish the solution if necessary.
			if norm(residual) > 1e-7
				D = LinearMap{T}((x) -> P*x + l*C*x, n; issymmetric=true)
				lsqr!(x, D, -q, maxiter=20)
				x /= sqrt(x'*C*x)/r
			end
			X = reshape(x, n, 1)
		end
	elseif is_first_pop # hard case
		hard_case = true
		if !direct
			X = extract_solution_hard_case(P, q, r, C, l, v1, v2, tol_hard)
		else
			X = extract_solution_hard_case_direct(P, q, r, C, l, v1, v2)
		end
		# Sometimes extract_solution_hard_case fails in this case only return the Lagrange multipliers
		if !isreal(X)
			X = zeros(T, n, 0)
		end
	end
	λ[idx] = -Inf  # This ensures that the next pop_solution! would not consider the same solution.

	return X, l*ones(size(X, 2)), hard_case
end

function extract_solution_hard_case(P, q::AbstractVector{T}, r::T, C, l::T,
	v1::AbstractVector{T}, v2::AbstractVector{T}, tol::T) where T

	n = length(q)
	D = LinearMap{T}((x) -> P*x + l*C*x + v2*dot(v2, x)/100, n; issymmetric=true)
	y = minres(D, -q, maxiter=100, verbose=false)
	residual = P*y + l*C*y + q
	if sqrt(y'*C*y) > r || norm(residual) > tol
		l += T(1e-10) # Usually perturbing the multiplier and solving again works
		minres!(y, D, -q, tol=norm(residual)/tol, verbose=false)
	end

	if sqrt(y'*C*y) < r
		α = roots(Polynomial([y'*(C*y) - r^2, 2*(C*v2)'*y, v2'*(C*v2)]))
		x1 = y + α[1]*v2
		x2 = y + α[2]*v2
		return [x1 x2]
	else
		x = r*y/sqrt(y'*C*y)
		return reshape(x, n, 1)
	end
end

function extract_solution_hard_case_direct(P, q::AbstractVector{T}, r::T, C, l::T,
	v1::AbstractVector{T}, v2::AbstractVector{T}) where T

	n = length(q)
	λ, V = eigen(Symmetric(P + l*C))
	λ[abs.(λ) .< 1e-9] .= 0
	A = V*diagm(0 => λ)*V' # Essentially P + l*C with a "refined" nullspace
	F = qr(A, Val(true))
	y = -(F\q) # y is the minimum norm solution of P + l*C = q
	if sqrt(y'*C*y) < r
		α = roots(Polynomial([y'*(C*y) - r^2, 2*(C*v2)'*y, v2'*(C*v2)]))
		x1 = y + α[1]*v2
		x2 = y + α[2]*v2
		return [x1 x2]
	else
		x = r*y/sqrt(y'*C*y)
		return reshape(x, n, 1)
	end
end


text:
using Pkg
Pkg.activate(".")
using Revise
using DCMQ
using DICOM
using NIfTI
using DSCMRI

mprage_uri = Dict{String,String}()
perfusion_uri = Dict{String,String}()

function dcmhandler(channel, d, uri)
    study_uid = d[tag"Study Instance UID"]
    if DSCMRI.isdsc(d)
        println("found perfusion series")
        perfusion_uri[study_uid] = uri
        if !(study_uid in keys(mprage_uri))
            return
        end
    elseif DSCMRI.ismprage(d)
        println("found mprage series")
        mprage_uri[study_uid] = uri
        if !(study_uid in keys(perfusion_uri))
            return
        end
    else
        return
    end

    filepath_perf = perfusion_uri[study_uid]
    filepath_mprage = mprage_uri[study_uid]

    println("reading file $filepath_perf")
    nii_perf = NIfTI.niread(filepath_perf)
    println("reading file $filepath_mprage")
    nii_mprage = NIfTI.niread(filepath_mprage)
    te = d[tag"Echo Time"]
    rawvol = DSCMRI.nifti2img(nii_perf)
    mprage = DSCMRI.nifti2img(nii_mprage)
    println("starting perfusion tumor workflow")
    @time perfmaps = DSCMRI.tumor_workflow(rawvol, mprage, te)

    for (key, val) in perfmaps
        vox = DSCMRI.img2nifti(val)
        outdir = ENV["HOME"] * "/.dimseweb/nii/" * d[tag"Study Instance UID"]
        filepath = "$outdir/$(d[tag"Series Instance UID"])_$key.nii"
        NIfTI.niwrite(filepath, vox)
        dcmout = copy(d)
        dcmout[tag"Series Description"] = "$key NOT FOR CLINICAL USE"
        dcmout[tag"Image Type"][1] = "DERIVED"
        dcmout[tag"Image Type"][2] = "SECONDARY"
        dcmout[tag"Image Type"][3] = "RESAMPLED"

        publish(channel, "stored.series.nii", dcmout, uri=filepath)
    end

    delete!(mprage_uri, study_uid)
    delete!(perfusion_uri, study_uid)
    return true
end

consumer_loop("localhost", "tumorperfusion", ["stored.series.nii"], dcmhandler)


text:
"""
```
prepare_forecast_inputs!(m, input_type, cond_type, output_vars;
    df = DataFrame(), verbose = :none)
```

Add required outputs using `add_requisite_output_vars` and load data if
necessary.

### Inputs

- `m::AbstractDSGEModel`: model object
- `input_type::Symbol`: See `?forecast_one`.
- `cond_type::Symbol`: See `?forecast_one`.
- `output_vars::Vector{Symbol}`: vector of desired output variables. See
  `?forecast_one_draw`

### Keyword Arguments

- `df::DataFrame`: historical data. If `cond_type in [:semi, :full]`, then the
   final row of `df` should be the period containing conditional data. If not
   provided, then `df` will be loaded using `load_data` with the appropriate
   `cond_type`
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`

### Outputs

- `output_vars`
- `df`
"""
function prepare_forecast_inputs!(m::AbstractDSGEModel{S},
    input_type::Symbol, cond_type::Symbol, output_vars::Vector{Symbol};
    df::DataFrame = DataFrame(), subset_inds::AbstractRange{Int64} = 1:0,
    check_empty_columns::Bool = true, bdd_fcast::Bool = true,
    verbose::Symbol = :none) where {S<:AbstractFloat}

    @assert cond_type in [:none, :semi, :full] "cond_type must be one of :none, :semi, or :full"

    implied_horizons = subtract_quarters(date_forecast_end(m), date_forecast_start(m)) + 1
    if implied_horizons != get_setting(m, :forecast_horizons)
        @warn "Setting :forecast_horizons does not match the number of periods implied by Settings :date_forecast_end" *
            " and :date_forecast_start. Updating :forecast_horizons to match"
        m <= Setting(:forecast_horizons, implied_horizons, "Number of periods to forecast ahead")
    end

    # Compute everything that will be needed to plot original output_vars
    output_vars = add_requisite_output_vars(output_vars, bdd_fcast = bdd_fcast)
    if input_type == :prior
        output_vars = setdiff(output_vars, [:bddforecastobs])
    end

    # Get products and classes computed
    output_prods   = unique(map(get_product, output_vars))
    output_classes = unique(map(get_class,   output_vars))

    # Throw error if input_type = :subset but no subset_inds provided
    if input_type == :subset && isempty(subset_inds)
        error("Must supply nonempty subset_inds if input_type = :subset")
    end

    # Throw error if trying to compute shock decompositions under an alternative
    # policy rule
    if alternative_policy(m).key != :historical &&
        any(prod -> prod in [:shockdec, :dettrend, :trend], output_prods)

        error("Only histories, forecasts, and IRFs can be computed under an alternative policy")
    end

    # Determine if we are only running IRFs. If so, we won't need to load data
    # below
    irfs_only = all(prod -> prod == :irf, output_prods)

    # Load data if not provided, else check data well-formed
    if !irfs_only
        if isempty(df)
            data_verbose = verbose == :none ? :none : :low
            df = load_data(m; cond_type = cond_type, try_disk = true, check_empty_columns = check_empty_columns, verbose = data_verbose)
        else
            @assert df[1, :date] == date_presample_start(m)
            @assert df[end, :date] == (cond_type == :none ? date_mainsample_end(m) : date_conditional_end(m))
        end
    end

    # If regime switching, check settings are consistent
    if haskey(m.settings, :regime_switching)
        if get_setting(m, :regime_switching)
            @assert get_setting(m, :n_regimes) == length(get_setting(m, :regime_dates)) "The number" *
                " in Setting :n_regimes ($(string(get_setting(m, :n_regimes))))" *
                " does not match the length of Setting :regime_dates ($(string(length(get_setting(m, :regime_dates)))))."
            @assert get_setting(m, :regime_dates)[1] == date_presample_start(m) "The first regime" *
                " date ($(string(get_setting(m, :regime_dates)[1]))) must match the presample start date " *
                "($(string(date_presample_start(m))))."
        end
    end

    return output_vars, df
end

"""
```
load_draws(m, input_type; subset_inds = 1:0, verbose = :low)

load_draws(m, input_type, block_inds; verbose = :low)
```

Load and return parameter draws from Metropolis-Hastings or SMC.

### Inputs

- `m::AbstractDSGEModel`: model object
- `input_type::Symbol`: one of the options for `input_type` described in the
  documentation for `forecast_one`
- `block_inds::AbstractRange{Int64}`: indices of the current block (already indexed by
  `jstep`) to be read in. Only used in second method

### Keyword Arguments

- `subset_inds::AbstractRange{Int64}`: indices specifying the subset of draws to be read
  in. Only used in first method
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`. If `:low` or greater, prints
  location of input file.

### Outputs

- `params`: first method returns a single parameter draw of type
  `Vector{Float64}`. Second method returns a `Vector{Vector{Float64}}` of
  parameter draws for this block.
"""
function load_draws(m::AbstractDSGEModel, input_type::Symbol;
                    subset_inds::AbstractRange{Int64} = 1:0,
                    verbose::Symbol = :low, filestring_addl::Vector{String} =
                    Vector{String}(undef, 0), use_highest_posterior_value::Bool = false,
                    input_file_name::String = "")

    if isempty(input_file_name)
        input_file_name = get_forecast_input_file(m, input_type, filestring_addl = filestring_addl)
    end
    println(verbose, :low, "Loading draws from $input_file_name")

    # Load single draw
    if input_type in [:mean, :mode, :mode_draw_shocks]
        if get_setting(m, :sampling_method) == :MH
            params = convert(Vector{Float64}, h5read(input_file_name, "params"))
        elseif get_setting(m, :sampling_method) == :SMC
            if input_type in [:mode, :mode_draw_shocks]
                # If mode is in forecast overrides, want to use the file specified
                if :mode in collect(keys(forecast_input_file_overrides(m)))
                    params = convert(Vector{Float64}, h5read(input_file_name, "params"))
                    # If not, load it from the cloud
                elseif (occursin("smc_paramsmode", input_file_name) &&
                        !use_highest_posterior_value && input_type == :mode)
                    params = convert(Vector{Float64}, h5read(input_file_name, "params"))
                else
                    # Check that input_file_name is correct. Note that if
                    # it already has smc_cloud, then the following code block does not do anything.
                    if occursin("smc_paramsmode", input_file_name) || occursin(".h5", input_file_name)
                        input_file_name = replace(replace(input_file_name,
                                                          "smc_paramsmode" => "smc_cloud"),
                                                  ".h5" => ".jld2")
                        println(verbose, :low, "Switching estimation file of draws to $input_file_name")
                    end

                    cloud = load(input_file_name, "cloud")
                    params = if typeof(cloud) <: Union{DSGE.Cloud,SMC.Cloud}
                        SMC.get_likeliest_particle_value(SMC.Cloud(cloud))
                    else
                        cloud.particles[argmax(get_logpost(cloud))].value
                    end
                end
            else
                error("SMC mean not implemented yet")
            end
        else
            error("Invalid :sampling method specification. Change in setting :sampling_method")
        end
    # Load full distribution
    elseif input_type == :full

        if get_setting(m, :sampling_method) == :MH
            params = map(Float64, h5read(input_file_name, "mhparams"))
        elseif get_setting(m, :sampling_method) == :SMC
            cloud = load(replace(replace(input_file_name, ".h5" => ".jld2"), "smcsave" => "smc_cloud"), "cloud")
            params_unweighted = get_vals(cloud)
            # Re-sample SMC draws according to their weights
            W = load(replace(replace(input_file_name, "smcsave" => "smc_cloud"), "h5" => "jld2"), "W")
            weights = W[:, end]
            inds = SMC.resample(weights)

            params = Matrix{Float64}(params_unweighted[:, inds]')
        else
            error("Invalid :sampling method specification. Change in setting :sampling_method")
        end

    # Load subset of full distribution
    elseif input_type == :subset

        if isempty(subset_inds)
            error("Must supply nonempty range of subset_inds if input_type == :subset")
        else
            if get_setting(m, :sampling_method) == :MH
                params = map(Float64, h5read(input_file_name, "mhparams", (subset_inds, :)))
            elseif get_setting(m, :sampling_method) == :SMC
                params = map(Float64, h5read(input_file_name, "smcparams", (subset_inds, :)))
            else
                throw("Invalid :sampling method specification. Change in setting :sampling_method")
            end
        end
    elseif input_type == :prior
        params = rand(m.parameters, n_forecast_draws(m, :prior))

    elseif input_type == :init || input_type == :init_draw_shocks
        # Return initial parameters of model object

#=        if m.spec == "het_dsge"
            init_parameters!(m, testing_gamma = false)
        else
            init_parameters!(m)
        end =#
        tmp = map(α -> α.value, m.parameters)
        params = convert(Vector{Float64}, tmp)

    end

    return params
end

function load_draws(m::AbstractDSGEModel, input_type::Symbol, block_inds::AbstractRange{Int64};
                    verbose::Symbol = :low,
                    filestring_addl::Vector{String} = Vector{String}(undef, 0),
                    input_file_name::String = "")

    if isempty(input_file_name)
        input_file_name = get_forecast_input_file(m, input_type, filestring_addl = filestring_addl)
    end
    println(verbose, :low, "Loading draws from $input_file_name")

    if input_type in [:full, :subset]
        if isempty(block_inds)
            error("Must supply nonempty range of block_inds for this load_draws method")
        else
            ndraws = length(block_inds)
            params = Vector{Vector{Float64}}(undef, ndraws)
            if get_setting(m, :sampling_method) == :SMC
                params_unweighted = similar(params)
            end
            for (i, j) in zip(1:ndraws, block_inds)
                if get_setting(m, :sampling_method) == :MH
                    params[i] = vec(map(Float64, h5read(input_file_name, "mhparams", (j, :))))
                elseif get_setting(m, :sampling_method) == :SMC

                    params_unweighted[i] = vec(map(Float64, h5read(input_file_name, "smcparams", (j, :))))
                else
                    throw("Invalid :sampling_method setting specification.")
                end
            end

            # Re-sample draws according to weights if the sampling_method was SMC
            if get_setting(m, :sampling_method) == :SMC
                # Re-sample SMC draws according to their weights
                @load replace(replace(input_file_name, "smcsave" => "smc_cloud"), "h5" => "jld2") W
                weights = W[:, end][block_inds]
                inds = SMC.resample(weights)

                params = params_unweighted[inds]
            end

            return params
        end
    elseif input_type == :init_draw_shocks
        ndraws = length(block_inds)
        params = repeat([map(x -> x.value, m.parameters)], ndraws)
    elseif input_type == :mode_draw_shocks
        ndraws = length(block_inds)
        params = repeat([map(x -> x.value, m.parameters)], ndraws)
    elseif input_type == :prior
        ndraws = length(block_inds)
        params = Vector{Vector{Float64}}(undef, ndraws)
        for (i, j) in zip(1:ndraws, block_inds)
            try_again = true
            param_try = vec(rand(m.parameters, ndraws)[:, i])
            # Need to do this to ensure that draws from prior don't return -Inf likelihood
            while try_again
                param_try = vec(rand(m.parameters, ndraws)[:, i])
                try
                    update!(m, param_try)
                    likelihood(m, rand(length(m.observables), 100))
                    try_again = false
                catch
                    try_again = true
                end
            end
            params[i] = param_try
        end
        return params
    else
        error("This load_draws method can only be called with input_type in [:full, :subset]")
    end
end

function load_draws(m::AbstractDSGEVARModel, input_type::Symbol; subset_inds::AbstractRange{Int64} = 1:0,
                    verbose::Symbol = :low,
                    filestring_addl::Vector{String} = Vector{String}(undef, 0),
                    use_highest_posterior_value::Bool = false,
                    input_file_name::String = "")

    if isempty(input_file_name)
        input_file_name = get_forecast_input_file(m, input_type;
                                                  filestring_addl = filestring_addl)
    end

    return load_draws(get_dsge(m), input_type; subset_inds = subset_inds, verbose = verbose,
                      filestring_addl = filestring_addl, use_highest_posterior_value =
                      use_highest_posterior_value,
                      input_file_name = input_file_name)
end

"""
```
forecast_one(m, input_type, cond_type, output_vars; df = DataFrame(),
    subset_inds = 1:0, forecast_string = "", verbose = :low, ...)
```

Compute and save `output_vars` for input draws given by `input_type` and
conditional data case given by `cond_type`.

### Inputs

- `m::AbstractDSGEModel`: model object

- `input_type::Symbol`: one of:

```
  - `:mode`: forecast using the modal parameters only
  - `:mean`: forecast using the mean parameters only
  - `:init`: forecast using the initial parameter values only
  - `:full`: forecast using all parameters (full distribution)
  - `:subset`: forecast using a well-defined user-specified subset of draws
```

- `cond_type::Symbol`: one of:

```
  - `:none`: no conditional data
  - `:semi`: use \"semiconditional data\" - average of quarter-to-date
    observations for high frequency series
  - `:full`: use \"conditional data\" - semiconditional plus nowcasts for
    desired observables
```

- `output_vars::Vector{Symbol}`: vector of desired output variables. See
  `?forecast_one_draw`.

### Keyword Arguments

- `df::DataFrame`: Historical data. If `cond_type in [:semi, :full]`, then the
   final row of `df` should be the period containing conditional data. If not
   provided, will be loaded using `load_data` with the appropriate `cond_type`
- `subset_inds::AbstractRange{Int64}`: indices specifying the draws we want to use. If a
  more sophisticated selection criterion is desired, the user is responsible for
  determining the indices corresponding to that criterion. If `input_type` is
  not `subset`, `subset_inds` will be ignored
- `forecast_string::String`: short string identifying the subset to be
  appended to the output filenames. If `input_type = :subset` and
  `forecast_string` is empty, an error is thrown.
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`.
- `check_empty_columns::Bool = true`: check empty columns or not when loading data (if `df` is empty)
- `zlb_method::Symbol`: method for enforcing the zero lower bound. Defaults to `:shock`,
    meaning we use a monetary policy shock to enforce the ZLB.

  Other available methods:
  1. `:temporary_altpolicy` -> use a temporary alternative policy to enforce the ZLB.

- `set_regime_vals_altpolicy::Function`: Function that adds new regimes to parameters when
    using temporary alternative policies (if needed). Defaults to identity (which does nothing)
    This function should take as inputs the model object `m` and the total number of regimes
     (after adding the required temporary regimes). It should then
     set up regime-switching parameters for these new additional regimes.
- `temporary_altpolicy_max_iter::Int`: Maximum number of iterations when enforcing
    the ZLB as a temporary alternative policy.
- `pegFFR::Bool = false`: peg the nominal FFR at the value specified by `FFRpeg`
- `FFRpeg::Float64 = -0.25/4`: value of the FFR peg
- `H::Int = 4`: number of horizons for which the FFR is pegged
- `bdd_fcast::Bool = true`: are we computing the bounded forecasts or not?
- `params::AbstractArray{Float64} = Vector{Float64}(undef, 0)`: parameter draws for the forecast.
     If empty, then we load draws from estimation files implied by the settings in `m`.

### Outputs

None. Output is saved to files returned by
`get_forecast_output_files(m, input_type, cond_type, output_vars)`.
"""
function forecast_one(m::AbstractDSGEModel{Float64},
                      input_type::Symbol, cond_type::Symbol, output_vars::Vector{Symbol};
                      df::DataFrame = DataFrame(), subset_inds::AbstractRange{Int64} = 1:0,
                      forecast_string::String = "",
                      use_filtered_shocks_in_shockdec::Bool = false,
                      shock_name::Symbol = :none, shock_var_name::Symbol = :none,
                      shock_var_value::Float64 = 0.0, check_empty_columns = true,
                      zlb_method::Symbol = :shock, set_regime_vals_altpolicy::Function = identity,
                      temporary_altpolicy_max_iter::Int = 10,
                      pegFFR::Bool = false, FFRpeg::Float64 = -0.25/4, H::Int = 4, bdd_fcast::Bool = true,
                      params::AbstractArray{Float64} = Vector{Float64}(undef, 0),
                      show_failed_percent::Bool = false,
                      verbose::Symbol = :low)

    ### Common Setup

    # Add necessary output_vars and load data
    output_vars, df = prepare_forecast_inputs!(m, input_type, cond_type, output_vars;
                                               df = df, verbose = verbose, bdd_fcast = bdd_fcast,
                                               subset_inds = subset_inds,
                                               check_empty_columns = check_empty_columns)

    # Get output file names
    forecast_output = Dict{Symbol, Array{Float64}}()
    forecast_output_files = get_forecast_output_files(m, input_type, cond_type, output_vars;
                                                      forecast_string = forecast_string)
    output_dir = rawpath(m, "forecast")

    # Determine if we are regime_switching
    regime_switching = haskey(get_settings(m), :regime_switching) ?
        get_setting(m, :regime_switching) : false
    n_regimes        = regime_switching && haskey(get_settings(m), :n_regimes) ?
        get_setting(m, :n_regimes) : 1

    # Print
    info_print(verbose, :low, "Forecasting input_type = $input_type, cond_type = $cond_type...")
    println(verbose, :low, "Start time: $(now())")
    println(verbose, :low, "Forecast outputs will be saved in $output_dir")

    ### Single-Draw Forecasts

    if input_type in [:mode, :mean, :init]

        elapsed_time = @elapsed let
            if isempty(params)
                params = load_draws(m, input_type; verbose = verbose)
            end
            forecast_output = forecast_one_draw(m, input_type, cond_type, output_vars,
                                                params, df, verbose = verbose,
                                                shock_name = shock_name,
                                                shock_var_name = shock_var_name,
                                                shock_var_value = shock_var_value,
                                                regime_switching = regime_switching,
                                                n_regimes = n_regimes, zlb_method = zlb_method,
                                                set_regime_vals_altpolicy = set_regime_vals_altpolicy,
                                                temporary_altpolicy_max_iter = temporary_altpolicy_max_iter,
                                                pegFFR = pegFFR, FFRpeg = FFRpeg, H = H)

            write_forecast_outputs(m, input_type, output_vars, forecast_output_files,
                                   forecast_output; df = df, block_number = Nullable{Int64}(),
                                   verbose = verbose)
        end

        write_forecast_outputs(m, input_type, output_vars, forecast_output_files,
                               forecast_output; df = df, block_number = Nullable{Int64}(),
                               verbose = verbose)

        total_forecast_time     = elapsed_time
        total_forecast_time_min = total_forecast_time/60
        println(verbose, :low, "\nTotal time to forecast: $total_forecast_time_min minutes")

    ### Multiple-Draw Forecasts

    elseif input_type in [:full, :subset, :prior, :init_draw_shocks, :mode_draw_shocks]
        if get_setting(m, :sampling_method) == :SMC
            subset_cond = input_type == :subset && get_jstep(m, length(subset_inds)) > 1
            other_cond = input_type in [:full, :prior, :init_draw_shocks, :mode_draw_shocks] &&
                get_jstep(m, n_forecast_draws(m, :full)) > 1
            if subset_cond || other_cond
                answer = ""
                invalid_answer = true
                while invalid_answer
                    println("WARNING: The forecast draws are being thinned by $(get_jstep(m, n_forecast_draws(m,:full))). This is not recommended for forecasts using an SMC estimation. Do you want to continue? (y/n)")
                    answer = readline(stdin)
                    if (answer == "y" || answer == "n")
                        invalid_answer = false
                    else
                        println("Invalid answer! Must be `y` or `n`")
                    end
                end

                if answer == "n"
                    throw(error("Forecast halted by user."))
                end
            end
        end

        if !isempty(params) & isa(params, AbstractMatrix)
            m <= Setting(:regime_switching_ndraws, size(params, 1))
        elseif !isempty(params) & (input_type == :mode_draw_shocks) & (!haskey(get_settings(m), :regime_switching_ndraws))
            error("If using :mode_draw_shocks with a user-defined matrix of parameter draws, " *
                  "then the user must manually add the Setting :regime_switching_ndraws " *
                  "to indicate the number of times shocks will be drawn.")
        end

        # Block info
        block_inds, block_inds_thin = forecast_block_inds(m, input_type; subset_inds = subset_inds)
        nblocks = length(block_inds)
        start_block = forecast_start_block(m)

        # Info needed for printing progress
        total_forecast_time = 0.0
        block_verbose = verbose == :none ? :none : :low

        for block = start_block:nblocks
            println(verbose, :low, )
            info_print(verbose, :low, "Forecasting block $block of $nblocks...")
            begin_time = time_ns()

            # Get to work!
            if isempty(params)
                params_for_map = load_draws(m, input_type, block_inds[block]; verbose = verbose)
            elseif input_type == :mode_draw_shocks
                ndraws = length(block_inds[block])
                @assert isa(params, Vector) "To use mode_draw_shocks with params passed in as a keyword, params must be a Vector."
                params_for_map = Vector{Float64}[params for i in block_inds[block]]
            else
                params_for_map = Vector{Float64}[params[i, :] for i in block_inds[block]]
            end

            mapfcn = use_parallel_workers(m) ? pmap : map
            forecast_outputs = mapfcn(param -> forecast_one_draw(m, input_type, cond_type, output_vars,
                                                                 param, df, verbose = verbose,
                                                                 use_filtered_shocks_in_shockdec =
                                                                 use_filtered_shocks_in_shockdec,
                                                                 shock_name = shock_name,
                                                                 shock_var_name = shock_var_name,
                                                                 shock_var_value = shock_var_value,
                                                                 regime_switching = regime_switching,
                                                                 n_regimes = n_regimes, zlb_method = zlb_method,
                                                                 set_regime_vals_altpolicy = set_regime_vals_altpolicy,
                                                                 temporary_altpolicy_max_iter = temporary_altpolicy_max_iter,
                                                                 pegFFR = pegFFR, FFRpeg = FFRpeg, H = H),
                                      params_for_map)

            # Assemble outputs from this block and write to file
            forecast_outputs = convert(Vector{Dict{Symbol, Array{Float64}}}, forecast_outputs)
            forecast_output = assemble_block_outputs(forecast_outputs; show_failed_percent = show_failed_percent)
            write_forecast_outputs(m, input_type, output_vars, forecast_output_files,
                                   forecast_output; df = df, block_number = Nullable(block),
                                   verbose = block_verbose, block_inds = block_inds_thin[block],
                                   subset_inds = subset_inds)
            GC.gc()

            # Calculate time to complete this block, average block time, and
            # expected time to completion
            block_time = (time_ns() - begin_time)/1e9
            total_forecast_time += block_time
            total_forecast_time_min     = total_forecast_time/60
            blocks_elapsed              = block - start_block + 1
            expected_time_remaining     = (total_forecast_time/blocks_elapsed)*(nblocks - block)
            expected_time_remaining_min = expected_time_remaining/60

            println(verbose, :low, "\nCompleted $block of $nblocks blocks.")
            println(verbose, :low, "Total time elapsed: $total_forecast_time_min minutes")
            println(verbose, :low, "Expected time remaining: $expected_time_remaining_min minutes")
        end # of loop through blocks

    end # of input_type
    combine_raw_forecast_output_and_metadata(m, forecast_output_files, verbose = verbose)

    println(verbose, :low, "\nForecast complete: $(now())")
end

"""
```
forecast_one_draw(m, input_type, cond_type, output_vars; params, df;
    verbose = :low)
```

Compute `output_vars` for a single parameter draw, `params`. Called by
`forecast_one`.

### Inputs

- `m::AbstractDSGEModel{Float64}`: model object
- `input_type::Symbol`: See `?forecast_one`.
- `cond_type::Symbol`: See `?forecast_one`.
- `output_vars::Vector{Symbol}`: vector of desired output variables. See Outputs
  section
- `params::Vector{Float64}`: parameter vector
- `df::DataFrame`: historical data.
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`.

### Output

- `forecast_outputs::Dict{Symbol, Array{Float64}}`: dictionary of forecast
  outputs. Keys are `output_vars`, which is some subset of:

```
  - `:histstates`: `Matrix{Float64}` of smoothed historical states
  - `:histobs`: `Matrix{Float64}` of smoothed historical data
  - `:histpseudo`: `Matrix{Float64}` of smoothed historical
    pseudo-observables
  - `:histshocks`: `Matrix{Float64}` of smoothed historical shocks
  - `:forecaststates`: `Matrix{Float64}` of forecasted states
  - `:forecastobs`: `Matrix{Float64}` of forecasted observables
  - `:forecastpseudo`: `Matrix{Float64}` of forecasted pseudo-observables
  - `:forecastshocks`: `Matrix{Float64}` of forecasted shocks
  - `:bddforecaststates`, `:bddforecastobs`, `:bddforecastpseudo`, and
    `:bddforecastshocks`: `Matrix{Float64}`s of forecasts where we enforce
    the zero lower bound to be `forecast_zlb_value(m)`
  - `:shockdecstates`: `Array{Float64, 3}` of state shock decompositions
  - `:shockdecobs`: `Array{Float64, 3}` of observable shock decompositions
  - `:shockdecpseudo`: `Array{Float64, 3}` of pseudo-observable shock
    decompositions
  - `:dettrendstates`: `Matrix{Float64}` of state deterministic trends
  - `:dettrendobs`: `Matrix{Float64}` of observable deterministic trends
  - `:dettrendpseudo`: `Matrix{Float64}` of pseudo-observable deterministic
    trends
  - `:trendstates`: `Vector{Float64}` of state trends, i.e. the `CCC` vector
  - `:trendobs`: `Vector{Float64}` of observable trends, i.e. the `DD` vector
  - `:trendpseudo`: `Vector{Float64}` of pseudo-observable trends, i.e. the
    `DD_pseudo` vector
  - `:irfstates`: `Array{Float64, 3}` of state impulse responses
  - `:irfobs`: `Array{Float64, 3}` of observable impulse responses
  - `:irfpseudo`: `Array{Float64, 3}` of pseudo-observable impulse responses
```
"""
function forecast_one_draw(m::AbstractDSGEModel{Float64}, input_type::Symbol, cond_type::Symbol,
                           output_vars::Vector{Symbol}, params::Vector{Float64}, df::DataFrame; verbose::Symbol = :low,
                           use_filtered_shocks_in_shockdec::Bool = false,
                           shock_name::Symbol = :none, shock_var_name::Symbol = :none,
                           shock_var_value::Float64 = 0.0, zlb_method::Symbol = :shock,
                           set_regime_vals_altpolicy::Function = identity,
                           temporary_altpolicy_max_iter::Int = 10,
                           pegFFR::Bool = false, FFRpeg::Float64 = -0.25/4, H::Int = 4,
                           param_key::Symbol = :nothing, param_value::Float64 = 0.0,
                           param_key2::Symbol = :nothing, param_value2::Float64 = 0.0,
                           regime_switching::Bool = false, n_regimes::Int = 1)
    ### Setup

    # Re-initialize model indices if forecasting under an alternative policy
    # rule (in case new states or equations were added)
    if alternative_policy(m).key != :historical
        init_model_indices!(m)
    end

    # Are we only running IRFs?
    output_prods = map(get_product, output_vars)
    irfs_only = all(x -> x == :irf, output_prods)

    # Compute state space
    update!(m, params) # Note that params is a Vector{Float64}, not a ParameterVector. This `update!` infers if the forecast is regime-switching if length(params) > length(m.parameters)
    system = compute_system(m)

    # Initialize output dictionary
    forecast_output = Dict{Symbol, Array{Float64}}()

    # Decide whether to draw states/shocks in smoother/forecast
    uncertainty_override = forecast_uncertainty_override(m)
    uncertainty = if Nullables.isnull(uncertainty_override)
        if input_type in [:init, :mode, :mean, :prior]
            false
        elseif input_type in [:full, :subset, :mode_draw_shocks, :init_draw_shocks]
            true
        end
    else
        get(uncertainty_override)
    end

    ### 1. Smoothed Histories

    # Must run smoother for conditional data in addition to explicit cases
    hist_vars = [:histstates, :histpseudo, :histshocks, :histstdshocks]
    shockdec_vars = [:shockdecstates, :shockdecpseudo, :shockdecobs]
    dettrend_vars = [:dettrendstates, :dettrendpseudo, :dettrendobs]
    smooth_vars = vcat(hist_vars, shockdec_vars, dettrend_vars)
    hists_to_compute = intersect(output_vars, smooth_vars)

    run_smoother = !isempty(hists_to_compute) ||
        (cond_type in [:semi, :full] && !irfs_only)
    if run_smoother
        # Call smoother
        histstates, histshocks, histpseudo, initial_states =
            smooth(m, df, system; cond_type = cond_type, draw_states = uncertainty)

        # For conditional data, transplant the obs/state/pseudo vectors from hist to forecast
        if cond_type in [:full, :semi]
            T = n_mainsample_periods(m)

            forecast_output[:histstates] = transplant_history(histstates, T)
            forecast_output[:histshocks] = transplant_history(histshocks, T)
            forecast_output[:histpseudo] = transplant_history(histpseudo, T)
        else
            forecast_output[:histstates] = histstates
            forecast_output[:histshocks] = histshocks
            forecast_output[:histpseudo] = histpseudo
        end

        # Standardize shocks if desired
        if :histstdshocks in output_vars
            if regime_switching
                start_date = max(date_mainsample_start(m), df[1, :date]) # use mainsample b/c shocks only includes mainsample, no presample
                end_date   = prev_quarter(date_forecast_start(m))
                regime_inds = regime_indices(m, start_date, end_date)
                if regime_inds[1][1] < 1
                    regime_inds[1] = 1:regime_inds[1][end]
                end
                forecast_output[:histstdshocks] =
                    standardize_shocks(forecast_output[:histshocks],
                                       Matrix{eltype(system[1, :QQ])}[system[i, :QQ] for i in 1:n_regimes], regime_inds)
            else
                forecast_output[:histstdshocks] = standardize_shocks(forecast_output[:histshocks], system[:QQ])
            end
        end
    end

    ### 2. Forecasts

    unbddforecast_vars = [:forecaststates, :forecastobs, :forecastpseudo, :forecastshocks, :forecaststdshocks]
    bddforecast_vars = [:bddforecaststates, :bddforecastobs, :bddforecastpseudo, :bddforecastshocks, :bddforecaststdshocks]
    forecast_vars = vcat(unbddforecast_vars, bddforecast_vars)
    forecasts_to_compute = intersect(output_vars, forecast_vars)

    if !isempty(forecasts_to_compute)
        # Get initial forecast state vector s_T
        s_T = if run_smoother # ONLY THIS BRANCH WORKS FOR REGIME SWITCHING
            # The last smoothed state is either s_{T|T} (if !uncertainty) or
            # drawn from N(s_{T|T}, P_{T|T}) (if uncertainty)
            histstates[:, end]
        else # THIS BRANCH HAS NOT BEEN EXTENDED YET FOR REGIME SWITCHING
            kal = Kalman(Vector{Float64}(undef,0), Matrix{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0, 0), Matrix{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0, 0), Vector{Float64}(undef, 0), Array{Float64}(undef, 0, 0, 0), Vector{Float64}(undef, 0), Array{Float64}(undef, 0, 0, 0))
            try
                kal = filter(m, df, system; cond_type = cond_type)
            catch err
                if isa(err, DomainError)
                    return Dict{Symbol, Array{Float64}}()
                else
                    rethrow(err)
                end
            end
            if uncertainty
                # If we want to draw s_T but haven't run the smoother, draw from
                # N(s_{T|T}, P_{T|T}) directly
                U, singular_values, _ = svd(kal[:P_T])
                dist = DegenerateMvNormal(kal[:s_T], U*Matrix(Diagonal(sqrt.(singular_values))))
                rand(dist)
            else
                # If we don't want to draw s_T, simply use the mean s_{T|T}
                kal[:s_T]
            end
        end

        if param_key!=:nothing
            m[param_key] = param_value
            steadystate!(m)
            system = compute_system(m)
        end
        if param_key2!=:nothing
            m[param_key2] = param_value2
            steadystate!(m)
            system = compute_system(m)
        end

        # Re-solve model with alternative policy rule, if applicable
        if alternative_policy(m).solve != solve
            system = compute_system(m; apply_altpolicy = true)
        end

        # Adjust the initial state vector for pgap and ygap
        if haskey(m.settings, :pgap_type) && haskey(get_settings(m), :pgap_value)
            if get_setting(m, :pgap_type) == :ngdp
                _, s_T = ngdp_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
            elseif get_setting(m, :pgap_type) == :ait
                _, s_T = ait_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
            end
            histstates[:, end] = s_T
            ZZ_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :ZZ_pseudo] : system[:ZZ_pseudo] # currently assumes these are
            DD_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :DD_pseudo] : system[:DD_pseudo] # the same across regimes
            histpseudo[:, end] = ZZ_pseudo * s_T + DD_pseudo
        end

        if haskey(m.settings, :ygap_type) && haskey(get_settings(m), :ygap_value) &&
            haskey(m.settings, :pgap_type) && haskey(get_settings(m), :pgap_value)
            if get_setting(m, :ygap_type) == :smooth_ait_gdp && get_setting(m, :pgap_type) == :smooth_ait_gdp
                _, s_T = smooth_ait_gdp_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
            elseif get_setting(m, :ygap_type) == :smooth_ait_gdp_alt && get_setting(m, :pgap_type) == :smooth_ait_gdp_alt
                _, s_T = smooth_ait_gdp_alt_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
            end
            histstates[:, end] = s_T
            ZZ_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :ZZ_pseudo] : system[:ZZ_pseudo] # currently assumes these are
            DD_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :DD_pseudo] : system[:DD_pseudo] # the same across regimes
            histpseudo[:, end] = ZZ_pseudo * s_T + DD_pseudo
        end

        if haskey(m.settings, :Rref_type)
            if get_setting(m, :ygap_type) == :rw && get_setting(m, :pgap_type) == :rw &&
                get_setting(m, :Rref_type) == :rw
                _, s_T = rw_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
            end
            histstates[:, end] = s_T
            ZZ_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :ZZ_pseudo] : system[:ZZ_pseudo] # currently assumes these are
            DD_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :DD_pseudo] : system[:DD_pseudo] # the same across regimes
            histpseudo[:, end] = ZZ_pseudo * s_T + DD_pseudo
        end

        # 2A. Unbounded forecasts
        if !isempty(intersect(output_vars, unbddforecast_vars))
            if pegFFR
                nshocks = size(system[:RRR], 2)
                nstates = size(system[:TTT], 1)
                shocks = m.exogenous_shocks
                PsiR1 = 0
                PsiR2 = zeros(nstates)
                PsiR2[m.endogenous_states[:R_t]] = 1.
                Rht = system[:RRR][:, vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")])]
                bb = zeros(H+1, 1)
                MH = zeros(H+1, H+1)
                for hh = 1:H+1
                    bb[hh, 1] = (FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_T)[1]
                    MH[hh, :] = PsiR2'*(system[:TTT])^(hh-1)*Rht
                end
                monshocks = MH\bb
                etpeg = zeros(nshocks, forecast_horizons(m))
                etpeg[vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")]), 1] = monshocks
                forecaststates, forecastobs, forecastpseudo, forecastshocks =
                    forecast(system, s_T, etpeg; cond_type = cond_type, enforce_zlb = false, draw_shocks = uncertainty)
                println("The forecasted interest rate path is $(forecastobs[m.observables[:obs_nominalrate], :])")
            else
                forecaststates, forecastobs, forecastpseudo, forecastshocks =
                    forecast(m, system, s_T;
                             cond_type = cond_type, enforce_zlb = false, draw_shocks = uncertainty)
            end

            # For conditional data, transplant the obs/state/pseudo vectors from hist to forecast
            if cond_type in [:full, :semi]
                forecast_output[:forecaststates] = transplant_forecast(histstates, forecaststates, T)
                forecast_output[:forecastshocks] = transplant_forecast(histshocks, forecastshocks, T)
                forecast_output[:forecastpseudo] = transplant_forecast(histpseudo, forecastpseudo, T)
                # NOTE: ZZ REGIME SWITCHING NOT SUPPORTED, SO JUST TAKE THE FIRST ZZ IN THE SYSTEM
                forecast_output[:forecastobs]    =
                    transplant_forecast_observables(histstates, forecastobs,
                                                    isa(system, RegimeSwitchingSystem) ? system[1] : system, T)
            else
                forecast_output[:forecaststates] = forecaststates
                forecast_output[:forecastshocks] = forecastshocks
                forecast_output[:forecastpseudo] = forecastpseudo
                forecast_output[:forecastobs]    = forecastobs
            end

            # Standardize shocks if desired
            if :forecaststdshocks in output_vars
                forecast_output[:forecaststdshocks] =
                    standardize_shocks(forecast_output[:forecastshocks],
                                       regime_switching ? system[n_regimes, :QQ] : system[:QQ])
            end
        end

        # 2B. Bounded forecasts

        if !isempty(intersect(output_vars, bddforecast_vars))
            if pegFFR
                nshocks = size(system[:RRR], 2)
                nstates = size(system[:TTT], 1)
                shocks = m.exogenous_shocks
                PsiR1 = 0
                PsiR2 = zeros(nstates)
                PsiR2[m.endogenous_states[:R_t]] = 1.
                Rht = system[:RRR][:, vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")])]
                bb = zeros(H+1, 1)
                MH = zeros(H+1, H+1)
                @show s_T[m.endogenous_states[:R_t]] + 400*log(m[:Rstarn])
                for hh = 1:H+1
                    bb[hh, 1] = (FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_T)[1]
                    MH[hh, :] = PsiR2'*(system[:TTT])^(hh-1)*Rht
                end
                monshocks = MH\bb
                etpeg = zeros(nshocks, forecast_horizons(m))
                etpeg[vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")]), 1] = monshocks
                forecaststates, forecastobs, forecastpseudo, forecastshocks = forecast(system, s_T, etpeg)
                @show forecaststates[m.endogenous_states[:R_t], :]
                @show forecastobs[m.observables[:obs_nominalrate], :]
            elseif zlb_method == :temporary_altpolicy
                altpolicy = get_setting(m, :alternative_policy).key
                @assert altpolicy in [:historical, :ait, :ngdp, :smooth_ait_gdp, :smooth_ait_gdp_alt] "altpolicy must be permanent " *
                    "and among [:historical, :ait, :ngdp, :smooth_ait_gdp, :smooth_ait_gdp_alt] for this method of enforcing the ZLB."

                # Run the unbounded forecast if they haven't already been computed
                if isempty(intersect(output_vars, unbddforecast_vars))
                    forecaststates, forecastobs, forecastpseudo, forecastshocks =
                        forecast(m, system, s_T;
                                 cond_type = cond_type, enforce_zlb = false, draw_shocks = uncertainty)
                end

                # Now run the ZLB enforcing forecast
                forecaststates, forecastobs, forecastpseudo =
                    forecast(m, altpolicy, s_T, forecaststates, forecastobs, forecastpseudo, forecastshocks; cond_type = cond_type,
                             temporary_altpolicy_max_iter = temporary_altpolicy_max_iter,
                             set_zlb_regime_vals = set_regime_vals_altpolicy)
            else
                forecaststates, forecastobs, forecastpseudo, forecastshocks =
                    forecast(m, system, s_T;
                             cond_type = cond_type, enforce_zlb = true, draw_shocks = uncertainty)
            end

            # For conditional data, transplant the obs/state/pseudo vectors from hist to forecast
            if cond_type in [:full, :semi]
                forecast_output[:bddforecaststates] = transplant_forecast(histstates, forecaststates, T)
                forecast_output[:bddforecastshocks] = transplant_forecast(histshocks, forecastshocks, T)
                forecast_output[:bddforecastpseudo] = transplant_forecast(histpseudo, forecastpseudo, T)
                # NOTE: ZZ REGIME SWITCHING NOT SUPPORTED, SO JUST TAKE THE FIRST ZZ IN TEH SYSTEM
                forecast_output[:bddforecastobs]    =
                    transplant_forecast_observables(histstates, forecastobs,
                                                    isa(system, RegimeSwitchingSystem) ? system[1] : system, T)
            else
                forecast_output[:bddforecaststates] = forecaststates
                forecast_output[:bddforecastshocks] = forecastshocks
                forecast_output[:bddforecastpseudo] = forecastpseudo
                forecast_output[:bddforecastobs]    = forecastobs
            end

            # Standardize shocks if desired
            if :bddforecaststdshocks in output_vars
                forecast_output[:bddforecaststdshocks] =
                    standardize_shocks(forecast_output[:bddforecastshocks],
                                       regime_switching ? system[n_regimes, :QQ] : system[:QQ])
            end
        end
    end

    ### 3. Shock Decompositions

    shockdecs_to_compute = intersect(output_vars, shockdec_vars)

    if !isempty(shockdecs_to_compute)
        histshocks_shockdec = if use_filtered_shocks_in_shockdec
            filter_shocks(m, df, system, cond_type = cond_type)
        else
            histshocks
        end

        start_date = max(date_mainsample_start(m), df[1, :date]) # smooth doesn't return presample
        end_date   = if cond_type in [:semi, :full] # end date of histshocks includes conditional periods
            max(date_conditional_end(m), prev_quarter(date_forecast_start(m)))
        else
            prev_quarter(date_forecast_start(m)) # this is the end date of history period
        end
        shockdecstates, shockdecobs, shockdecpseudo = isa(system, RegimeSwitchingSystem) ?
            shock_decompositions(m, system, histshocks_shockdec, start_date, end_date) :
            shock_decompositions(m, system, histshocks_shockdec)

        forecast_output[:shockdecstates] = shockdecstates
        forecast_output[:shockdecobs]    = shockdecobs
        forecast_output[:shockdecpseudo] = shockdecpseudo
    end

    ### 4. Trend

    trend_vars = [:trendstates, :trendobs, :trendpseudo]
    trends_to_compute = intersect(output_vars, trend_vars)

    if !isempty(trends_to_compute)
        trendstates, trendobs, trendpseudo = trends(system)

        forecast_output[:trendstates] = trendstates
        forecast_output[:trendobs]    = trendobs
        forecast_output[:trendpseudo] = trendpseudo
    end

    ### 5. Deterministic Trend

    dettrends_to_compute = intersect(output_vars, dettrend_vars)

    if !isempty(dettrends_to_compute)
        if isa(system, RegimeSwitchingSystem)
            start_date = max(date_mainsample_start(m), df[1, :date]) # smooth doesn't return presample
            end_date   = if cond_type in [:semi, :full] # end date of histshocks includes conditional periods
                max(date_conditional_end(m), prev_quarter(date_forecast_start(m)))
            else
                prev_quarter(date_forecast_start(m)) # this is the end date of history period
            end
            dettrendstates, dettrendobs, dettrendpseudo = deterministic_trends(m, system, initial_states, start_date, end_date)
        else
            dettrendstates, dettrendobs, dettrendpseudo = deterministic_trends(m, system, initial_states)
        end

        forecast_output[:dettrendstates] = dettrendstates
        forecast_output[:dettrendobs]    = dettrendobs
        forecast_output[:dettrendpseudo] = dettrendpseudo
    end

    ### 6. Impulse Responses

    irf_vars = [:irfstates, :irfobs, :irfpseudo]
    irfs_to_compute = intersect(output_vars, irf_vars)

    if !isempty(irfs_to_compute)
        if shock_name!=:none
            irfstates, irfobs, irfpseudo = impulse_responses(m, system, impulse_response_horizons(m),
                                                             shock_name,
                                                             shock_var_name,
                                                             shock_var_value)
        else
            irfstates, irfobs, irfpseudo = impulse_responses(m, system)
        end
        forecast_output[:irfstates] = irfstates
        forecast_output[:irfobs] = irfobs
        forecast_output[:irfpseudo] = irfpseudo
    end

    ### Return only desired output_vars

    for key in keys(forecast_output)
        if !(key in output_vars)
            delete!(forecast_output, key)
        end
    end
    return forecast_output
end


text:
using Unitful
using UnitfulEquivalences
using PhysicalConstants.CODATA2018
using Interpolations

import Unitful:
    nm, μm, mm, cm, m, km,
    mg, g, kg,
    ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
    eV,
    μJ, mJ, J,
	μW, mW, W,
	A, N, mol, mmol, V, L, M

	import PhysicalConstants.CODATA2018: N_A

"""
	struct Fov

Represent a field of view

# Fields
- `d::Unitful.Length`  : diameter of Fov
- `z::Unitful.Length`  : thickness
- `a::Unitful.Area`    : area (computed)
- `v::Unitful.Volume`  : volume (computed)

"""
struct Fov
    d::Unitful.Length
    z::Unitful.Length
	a::Unitful.Area
    v::Unitful.Volume

	function Fov(d,z)
		a = π * (d/2.)^2
		v = a * z
		new(d,z,a,v)
	end
end

"""
	struct Fluorophore

Represent a fluorescent molecule

# Fields
- `expeak::Units (nm)`           : peak excitation
- `empeak::Units (nm)`           : peak emission
- `ϵ::Units(``cm^{-1} M^{-1}``)` : molar extinction coefficient
- `Q::Float64`                   : Quantum efficiency
- `σ::Units(``cm^2``)`           : attenuation cross section
"""
struct Fluorophore
    expeak::typeof(1.0nm)
    empeak::typeof(1.0nm)
    ϵ::typeof(1.0/(cm*M))
    Q::Float64
	σ::typeof(1.0cm^2)
	function Fluorophore(ex, en, ϵ, Q)
		σ = log(10) * uconvert(cm^2/mol, ϵ) / N_A
		new(ex, en, ϵ, Q, σ)
	end
end

"""
	struct Monolayer

Represents a monolayer

# Fields
- `p::Unitful.Length`  : pitch between molecules
- `σ::Units (``molecules/μm^2``)`  : number of molecules per area
"""
struct Monolayer
	p::Unitful.Length
    σ::typeof(1.0*μm^-2)

	function Monolayer(p)
		σ = uconvert(μm^-2, 1.0/p^2)
		new(p, σ)

	end
end

abstract type Laser end

"""
	struct CLaser <: Laser

Simple representation of a continous laser

# Fields
- `λ::typeof(1.0nm)`  : Laser wavelength
- `P::typeof(1.0mW)`  : Power

"""
struct CLaser <: Laser
	λ::Unitful.Length
	P::Unitful.Power
end

"""
	struct PulsedLaser

Simple representation of a pulsed laser

# Fields
- `λ::Unitful.Length`  : Laser wavelength
- `Pk::Unitful.Power`  : peak Power
- `P::Unitful.Power`  : average power
- `f::Unitful.Frequency`  : laser frequency
- `w::Unitful.Time` : pulse width

"""
struct PulsedLaser  <: Laser
	λ::Unitful.Length
	Pk::Unitful.Power
	f::Unitful.Frequency
	P::Unitful.Power
	w::Unitful.Time

	function PulsedLaser(λ::Unitful.Length,
						 Pk::Unitful.Power,
						 f::Unitful.Frequency,
		                 w::Unitful.Time)

		P = uconvert(mW, Pk * w * f)
		new(λ, Pk, f, P,  w)
	end

	function PulsedLaser(λ::Unitful.Length,
						 Pk::Unitful.Power,
						 f::Unitful.Frequency,
						 P::Unitful.Power,
						 w::Unitful.Time)

			new(λ, Pk, f, P,  w)
	end
end

# """
# GaussianLaser

# Representation of a Gaussian laser defined by a laser and a waist

# # Fields
# - `laser::Laser`       : A Laser
# - `w0::Unitful.Length` : Waist of laser

# """
# struct GaussianLaser
# 	laser::Laser
# 	w0::Unitful.Length
# end

"""
	struct TLens

Represents a thin lens

# Fields
- `M::Real`         : Magnification
- `f::Unitful.Length`   : focal distance of the lens
- `d::Unitful.Length`   : diameter of the lens

"""
struct TLens
    M::Real
	f::Unitful.Length
	d::Unitful.Length

end

"""
	propagate(tl::TLens, y1::Real, z1::Real)

Propagate ray from (z1, y1) to (z2, y2) through thin lens tl

# Fields
- `tl::TLens`         : A thin lens
- `z1::Unitful.Length`   : point z1
- `y1::Unitful.Length`   : point y1

"""
function propagate(tl::TLens, y1::Unitful.Length, z1::Unitful.Length)
	y2 = tl.M * y1
	iz2 = (1.0/tl.f) - (1.0/z1)
	y2, 1.0/iz2
end

"""
	struct Objective

Simple representation of a microscope objective

# Fields
- `name::String`       : identifies the objective
- `NA::Float64`        : Numerical aperture
- `M::Float64`         : Magnification
- `f::typeof(1.0mm)`   : focal length
- `d::typeof(1.0mm)`   : entrance pupil diameter

"""
struct Objective
	name::String
    NA::Real
    M::Real
	f::Unitful.Length
	d::Unitful.Length

	function Objective(name::String, NA::Float64,  M::Float64)
		new(name, NA, M, -1.0mm, -1.0mm)
	end

	function Objective(name::String, f::Unitful.Length, d::Unitful.Length, M::Float64)
		ff = f / d   # https://www.eckop.com/resources/optics/numerical-aperture-and-f-number/
		NA = 1.0/(2.0*ff)
		new(name, NA, M, f,d)
	end
end

#CCD is defined in terms of a function which returns the CCD response
#(e.g, function ccd returns a response function, which gives efficiency
# as a function of wavelength)
"""
	ccd(lmin::Float64=350.0, lmax::Float64=1000.0)

Return the efficiency of a generic CCD as a function of wavelength.

# Fields

- `lmin::Float64=350.0` : Minimum wavelength for which efficiency is defined
- `lmax::Float64=350.0` : Maximum wavelength for which efficiency is defined

"""
function ccd(lmin::Float64=350.0, lmax::Float64=1000.0)
	function eff(l::Float64)
		if l < lmin || l > lmax
			return 0.
		else
			wl = 350.:50.:1000.
			ϵ = [0.3, 0.4,0.65,0.78,0.82,0.82,0.8,0.72,0.62,0.5,0.37,
			  0.24,0.12,0.07]
			e = CubicSplineInterpolation(wl, ϵ)
			return e(l)
		end
	end
	return eff
end

#FUNCTIONS

"""
	photon_energy(λ::Unitful.Length)

Given wavelength of photon return its energy.
# Fields

- `λ::Unitful.Length`  : Photon wavelength

"""
function photon_energy(λ::Unitful.Length)
	uconvert(eV, λ, Spectral())
end

"""
	delivered_energy(laser::Laser, t::Unitful.Time)

Delivered energy of a laser in a given time.
# Fields

- `laser::Laser`     : Laser
- `t::Unitful.Time`  : Time in which target is illuminated

"""
function delivered_energy(laser::Laser, t::Unitful.Time)
	laser.P * t
end

"""
	n_photons(laser::Laser)

Rate of photons (number of photons per unit time) produced by a laser
# Fields

- `laser::Laser`     : Laser
- `t::Unitful.Time`  : Time in which target is illuminated

"""
function n_photons(laser::Laser)
	uconvert(Hz, laser.P / photon_energy(laser.λ))
end

"""
	n_photons(λ::Unitful.Length, p::Unitful.Power)

Rate of photons (number of photons per unit time) corresponding to a wavelength
λ and a power P

# Fields

- `λ::Unitful.Length` : photon wavelength
- `p::Unitful.Power`  : Power

"""
function n_photons(λ::Unitful.Length, p::Unitful.Power)
	uconvert(Hz, p / photon_energy(λ))
end

"""
	n_photons_int(laser::Laser, t::Unitful.Time)

Integrated number of photons in a given time emitted by a laser

# Fields

- `laser::Laser`    : Laser
- `t::Unitful.Time` : time of measurement

"""
function n_photons_int(laser::Laser, t::Unitful.Time)
	uconvert(eV,delivered_energy(laser, t)) / photon_energy(laser.λ)
end

"""
	photon_density(λ::Unitful.Length, p::Unitful.Power, a::Unitful.Area)

number of photons per unit time per unit area

# Fields

- `λ::Unitful.Length` : photon wavelength
- `p::Unitful.Power`  : Power
- `a::Unitful.Area`   : Area

"""
function photon_density(λ::Unitful.Length, p::Unitful.Power, a::Unitful.Area)
	return n_photons(λ, p)/ a
end

"""
	photon_density(l::Laser, fov::Fov)

number of photons per unit time per unit area, in a Fov illuminated by a laser

# Fields

- `laser::Laser` : Laser
- `fov::Fov`     : Field of view

"""
function photon_density(laser::Laser, fov::Fov)
	return n_photons(laser) / fov.a
end

"""
	diffraction_limit(l::Laser, obj:: Objective)

Return the diameter of diffractive spot for a laser l
focused with an objective obj

# Fields

- `l::Laser`         : A laser
- `obj:: Objective`  : An objective
"""
function diffraction_limit(l::Laser, obj:: Objective)
    return 1.83 * l.λ/(2 * obj.NA)
end

"""
	geometrical_acceptance(d::Float64, D::Float64)

Compute the fraction of photons that make it through an iris
of diameter D located at a distance d from the emission point.

# Fields

- `d::Float64`   : distance between emission point and iris
- `D::Float64`   : Diameter of iris

"""
function geometrical_acceptance(d::Float64, D::Float64)
	return 0.5(1. - d/sqrt(d^2 + (D/2.)^2))
end

"""
	transmission(objective::Objective)

Compute the transmission of an objective (depends only of NA).

# Fields

- `objective::Objective` : Objective

"""
function transmission(objective::Objective)
	A = objective.NA
	A <=1 ? (1 - sqrt(1 - A^2)) /2 : 0.5
end

"""
	transmission(A::Float64)

Compute the transmission as a function of NA.

# Fields

- `A::Float64` : Numerical acceptance (NA)

"""
function transmission(A::Float64)
	A <=1 ? (1 - sqrt(1 - A^2)) /2 : 0.5
end

# fluorophores

"""
	fluorescence(f::Fluorophore, I::Quantity)

Number of  photons emitted per unit time when fluorosphore F
is illuminated with a laser of photon density I

# Fields

- `f::Fluorophore`    : Define fluorophore
- `I::(``nγ/cm^2``)`  : Photon density
"""
fluorescence(f::Fluorophore, I::typeof(1.0Hz/cm^2)) =  f.σ * f.Q * I

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab