text:
module ARIMA
#-----------------------------------------------------------
#PURPOSE
#	Colin T. Bowers module for ARIMA modelling
#NOTES
#	This package contains all the functions from the RARIMA package that don't depend on the R language via RCall. Currently, this includes simulation of ARIMA models, and point forecasts (no confidence bounds) of ARIMA models. For estimation and forecasts with confidence bounds, RARIMA is still necessary.
#SUB-OPTIMAL ROUTINES
#	The arima_statdiff and arima_statdiffinv routines create a temporary array for every order of differencing/integration. There is probably a smarter way to do this.
#-----------------------------------------------------------

#Load any specific variables/functions that are needed (use import ModuleName1.FunctionName1, ModuleName2.FunctionName2, etc)
import 	Base.show,
  		Base.copy,
		Base.deepcopy

#Specify the variables/functions to export (use export FunctionName1, FunctionName2, etc)
export	ARIMAInput, #Type that specifies input parameters for estimating an ARIMA model
		AutoARIMAInput, #Type that specifies input parameters for automatically selecting an ARIMA model then estimating it
		ARIMAModel, #Output of estimation (or auto-estimation) of an ARIMA model, ie complete specification of a model
		estimate, #Function for estimating a model given some data
		estimateARIMA, #Keyword wrapper on estimate function
		estimateAutoARIMA, #Keyword wrapper on estimate function
		forecast, #Function for forecasting using a model and data
		simulate, #Function for simulating a model
		simulateARIMA, #Keyword wrapper on simulate function
		coef, #Function for retrieving coefficients of ARIMAModel
		vcov, #Function for retrieving variance-covariance matrix of parameters of ARIMAModel
		coefVar #Function for retrieving coefficients of ARIMAModel along with their variance estimates

#******************************************************************************

#----------------------------------------------------------
#SET CONSTANTS FOR MODULE
#----------------------------------------------------------
const validMethodType = Union(ASCIIString, Symbol)::UnionType #type union nesting valid ways to express an optimization/estimation method

#----------------------------------------------------------
#SET ABSTRACT TYPES FOR MODULE
#----------------------------------------------------------
#None

#----------------------------------------------------------
#TYPE
#	ARIMAInput
#FIELDS
#	p::Int: Number of AR lags
#	d::Int: Order of integration
#	q::Int: Number of MA lags
#	P::Int: Number of seasonal AR lags
#	D::Int: Order of seasonal integration
#	Q::Int: Number of seasonal MA lags
#	seasonalPeriod::Int: Seasonal period of data (only relevant if P, D, or Q are > 0)
#	includeIntercept::Bool: Include intercept term in model. Not relevant for d > 0.
#	transformParams::Bool: AR parameters transformed to ensure they remain in region of stationarity
#	fixedParams::Vector{Float64}: Array of size p+q+P+Q containing fixed values for coefficients. Set value to NaN if you don't want to fix a particular parameter. Set to empty array to not fix any parameters during estimation.
#	initParams::Vector{Float64}: Array of size p+q+P+Q of initial values for estimation. Set value to NaN if you don't have a initial value for a particular parameter (0 will be used). Set to empty array to not supply any initial values.
#	method<:validMethodType: Variable that indicates method of estimating parameters. If ASCIIString, set to "CSS-ML", "CSS", or "ML". See R docs for arima function for more detail. Other types for this variable not coded up yet.
#PURPOSE
#	This type stores all of the input, excluding the data, needed to estimate the parameters of an ARIMA model.
#CONSTRUCTORS
#	ARIMAInput{T<:validMethodType}(p::Int, d::Int, q::Int, P::Int, D::Int, Q::Int, seasonalPeriod::Int, includeIntercept::Bool, transformParams::Bool, fixedParams::Vector{Float64}, initParams::Vector{Float64}, method::T): Constructor for every field
#	ARIMAInput(): White noise constructor
#KEYWORD CONSTRUCTOR
#	ARIMAInput(; p::Int=0, d::Int=0, q::Int=0, P::Int=0, D::Int=0, Q::Int=0, seasonalPeriod::Int=1, includeIntercept::Bool=true, transformParams::Bool=true, fixedParams::Vector{Float64}=Array(Float64, 0), initParams::Vector{Float64}=Array(Float64, 0), method::ASCIIString="CSS-ML")
#METHODS
#	copy(a::ARIMAInput): Standard copy
#	deepcopy(a::ARIMAInput): Standard deepcopy
#	show(io::IO, a::ARIMAInput): Standard show
#	show(a::ARIMAInput): Show defaulting to STDOUT
#NOTES
#	Any new estimation methods can be allowed for by using a custom type for the method field
#----------------------------------------------------------
#--------------------
#Type definition
#--------------------
type ARIMAInput{T<:validMethodType}
	p::Int
	d::Int
	q::Int
	P::Int
	D::Int
	Q::Int
	seasonalPeriod::Int
	includeIntercept::Bool
	transformParams::Bool
	fixedParams::Vector{Float64}
	initParams::Vector{Float64}
	method::T
	function ARIMAInput{T<:validMethodType}(p::Int, d::Int, q::Int, P::Int, D::Int, Q::Int, seasonalPeriod::Int, includeIntercept::Bool, transformParams::Bool, fixedParams::Vector{Float64}, initParams::Vector{Float64}, method::T)
		(p < 0 || d < 0 || q < 0 || P < 0 || D < 0 || Q < 0) && error("p, d, q, P, D, and Q must all be non-negative")
		seasonalPeriod < 0 && error("Seasonal period must be non-negative")
		if includeIntercept == true
			if d > 0 || D > 0
				println("WARNING: includeIntercept altered to false since either d > 0 or D > 0")
				includeIntercept = false
			end
		end
		numCoefTemp = p + q + P + Q
		if length(fixedParams) > 0
			length(fixedParams) != numCoefTemp && error("Fixed parameter vector is different length to number of model coefficients")
		end
		if length(initParams) > 0
			length(initParams) != numCoefTemp && error("Initialization parameter vector is different length to number of model coefficients")
		end
		if typeof(method) == ASCIIString
			(method != "CSS-ML" && method != "ML" && method != "CSS") && error("Invalid estimation method")
		elseif typeof(method) == Symbol
			error("method field of type Symbol has not been coded yet")
		else
			error("method field is an unexpected type")
		end
		new(p, d, q, P, D, Q, seasonalPeriod, includeIntercept, transformParams, fixedParams, initParams, method)
	end
end
#--------------------
#Constructors (Note, one constructor is defined after ARIMAModel type definition)
#--------------------
ARIMAInput{T<:validMethodType}(p::Int, d::Int, q::Int, P::Int, D::Int, Q::Int, seasonalPeriod::Int, includeIntercept::Bool, transformParams::Bool, fixedParams::Vector{Float64}, initParams::Vector{Float64}, method::T) = ARIMAInput{typeof(method)}(p, d, q, P, D, Q, seasonalPeriod, includeIntercept, transformParams, fixedParams, initParams, method)
ARIMAInput() = ARIMAInput(0, 0, 0, 0, 0, 0, 1, true, true, Array(Float64, 0), Array(Float64, 0), "CSS-ML")
#--------------------
#Keyword constructor
#--------------------
ARIMAInput(; p::Int=0, d::Int=0, q::Int=0, P::Int=0, D::Int=0, Q::Int=0, seasonalPeriod::Int=1, includeIntercept::Bool=true, transformParams::Bool=true, fixedParams::Vector{Float64}=Array(Float64, 0), initParams::Vector{Float64}=Array(Float64, 0), method::ASCIIString="CSS-ML") = ARIMAInput(p, d, q, P, D, Q, seasonalPeriod, includeIntercept, transformParams, fixedParams, initParams, method)
#--------------------
#Methods
#--------------------
#copy, deepcopy
copy(a::ARIMAInput) = ARIMAInput(copy(a.p), copy(a.d), copy(a.q), copy(a.P), copy(a.D), copy(a.Q), copy(a.seasonalPeriod), copy(a.includeIntercept), copy(a.transformParams), copy(a.fixedParams), copy(a.initParams), copy(a.method))
deepcopy(a::ARIMAInput) = ARIMAInput(deepcopy(a.p), deepcopy(a.d), deepcopy(a.q), deepcopy(a.P), deepcopy(a.D), deepcopy(a.Q), deepcopy(a.seasonalPeriod), deepcopy(a.includeIntercept), deepcopy(a.transformParams), deepcopy(a.fixedParams), deepcopy(a.initParams), deepcopy(a.method))
#show
function show(io::IO, a::ARIMAInput)
	println(io, "ARIMAInput contents:")
	println(io, "    p = " * string(a.p))
	println(io, "    d = " * string(a.d))
	println(io, "    q = " * string(a.q))
	println(io, "    P = " * string(a.P))
	println(io, "    D = " * string(a.D))
	println(io, "    Q = " * string(a.Q))
	println(io, "    seasonal period = " * string(a.seasonalPeriod))
	println(io, "    include intercept = " * string(a.includeIntercept))
	println(io, "    transform params for stationarity = " * string(a.transformParams))
	if length(a.fixedParams) == 0
		println(io, "    fixed parameters = none")
	else
		println(io, "    fixed parameters = some allocated")
	end
	if length(a.initParams) == 0
		println(io, "    initial parameters = none")
	else
		println(io, "    initial parameters = some allocated")
	end
	println(io, "    optimisation method = " * string(a.method))
end
show(a::ARIMAInput) = show(STDOUT, a)

#----------------------------------------------------------
#TYPE
#	AutoARIMAInput
#FIELDS
# 	d::Int: order of integration (use to fix a specific value). Set to negative value if d not fixed.
# 	D::Int: order of seasonal integration (use to fix a specific value). Set to negative value if D not fixed.
# 	pMax::Int: maximum number of AR lags allowed in search procedure
# 	qMax::Int: maximum number of MA lags allowed in search procedure
# 	PMax::Int: maximum number of seasonal AR lags allowed in search procedure
# 	QMax::Int: maximum number of seasonal MA lags allowed in search procedure
# 	numCoefMax::Int: maximum number of coefficients total allowed in search procedure (not relevant for step-wise search procedure)
# 	dMax::Int: maximum order of integration allowed in search procedure (over-ridden by d)
# 	DMax::Int: maximum seasonal order of integration allowed in search procedure (over-ridden by D)
# 	pStart::Int: Initial value for p in step-wise search procedure
# 	qStart::Int: Initial value for q in step-wise search procedure
# 	PStart::Int: Initial value for P in step-wise search procedure
# 	QStart::Int: Initial value for Q in step-wise search procedure
# 	stationary::Bool: Set to true to restrict search procedure to stationary models
# 	seasonal::Bool: Set to false to restrict search procedure to non-seasonal models
# 	ic::ASCIIString: Information criteria to use in search procedure. Use "aic", "aicc", or "bic" (see R auto-arima docs for more detail)
# 	stepwise::Bool: If true, use step-wise search procedure. If false, check every feasible model (slower).
# 	parallel::Bool: If true and stepwise=false, then do checks in parallel if possible.
# 	numCore::Int: Number of cores to use if doing checks in parallel.
#PURPOSE
#	This type stores all of the input, excluding the data, needed to automatically select an ARIMA model and estimate it
#CONSTRUCTORS
#	AutoARIMAInput(d::Int, D::Int, pMax::Int, qMax::Int, PMax::Int, QMax::Int, numCoefMax::Int, dMax::Int, DMax::Int, pStart::Int, qStart::Int, PStart::Int, QStart::Int, stationary::Bool, seasonal::Bool): Default methodology options
#	AutoARIMAInput(): No seasonality, all method options set to default, all order and step-wise start values set to default.
#KEYWORD CONSTRUCTOR
#	AutoARIMAInput(; d::Int=-1, D::Int=-1, pMax::Int=5, qMax::Int=5, PMax::Int=2, QMax::Int=2, numCoefMax::Int=5, dMax::Int=2, DMax::Int=1, pStart::Int=2, qStart::Int=2, PStart::Int=1, QStart::Int=5, stationary::Bool=false, seasonal::Bool=false, ic::ASCIIString="aic", stepwise::Bool=true, parallel::Bool=false, numCore::Int=2)
#METHODS
#	copy(a::AutoARIMAInput): Standard copy
#	deepcopy(a::AutoARIMAInput): Standard deepcopy
#	show(io::IO, a::AutoARIMAInput): Standard show
#	show(a::AutoARIMAInput): Show defaulting to STDOUT
#NOTES
#----------------------------------------------------------
#--------------------
#Type definition
#--------------------
type AutoARIMAInput
	d::Int
	D::Int
	pMax::Int
	qMax::Int
	PMax::Int
	QMax::Int
	numCoefMax::Int
	dMax::Int
	DMax::Int
	pStart::Int
	qStart::Int
	PStart::Int
	QStart::Int
	stationary::Bool
	seasonal::Bool
	ic::ASCIIString
	stepwise::Bool
	parallel::Bool
	numCore::Int
	function AutoARIMAInput(d::Int, D::Int, pMax::Int, qMax::Int, PMax::Int, QMax::Int, numCoefMax::Int, dMax::Int, DMax::Int, pStart::Int, qStart::Int, PStart::Int, QStart::Int, stationary::Bool, seasonal::Bool, ic::ASCIIString, stepwise::Bool, parallel::Bool, numCore::Int)
		(pMax < 0 || qMax < 0 || PMax < 0 || QMax < 0 || numCoefMax < 0 || dMax < 0 || DMax < 0) && error("maximum value non-negativity constraint violated in inner constructor")
		(pStart < 0 || qStart < 0 || PStart < 0 || QStart < 0) && error("step-wise starting value non-negativity constraint violated in inner constructor")
		(ic != "aicc" && ic != "aic" && ic != "bic") && error("invalid value for ic field")
		if parallel == true && stepwise == true
			println("WARNING: parallel computation incompatible with stepwise search method. parallel altered to false.")
			parallel = false
		end
		numCore < 1 && error("numCore must be strictly positive")
		new(d, D, pMax, qMax, PMax, QMax, numCoefMax, dMax, DMax, pStart, qStart, PStart, QStart, stationary, seasonal, ic, stepwise, parallel, numCore)
	end
end
#--------------------
#Constructors
#--------------------
AutoARIMAInput() = AutoARIMAInput((-1, -1, 5, 5, 2, 2, 5, 2, 1, 2, 2, 1, 1, false, false, "aic", true, false, 2))
#--------------------
#Keyword constructor
#--------------------
AutoARIMAInput(; d::Int=-1, D::Int=-1, pMax::Int=5, qMax::Int=5, PMax::Int=2, QMax::Int=2, numCoefMax::Int=5, dMax::Int=2, DMax::Int=1, pStart::Int=2, qStart::Int=2, PStart::Int=1, QStart::Int=5, stationary::Bool=false, seasonal::Bool=false, ic::ASCIIString="aic", stepwise::Bool=true, parallel::Bool=false, numCore::Int=2) = AutoARIMAInput(d, D, pMax, qMax, PMax, QMax, numCoefMax, dMax, DMax, pStart, qStart, PStart, QStart, stationary, seasonal, ic, stepwise, parallel, numCore)
#--------------------
#Methods
#--------------------
#copy, deepcopy
copy(a::AutoARIMAInput) = AutoARIMAInput(copy(a.d), copy(a.D), copy(a.pMax), copy(a.qMax), copy(a.PMax), copy(a.QMax), copy(a.numCoefMax), copy(a.dMax), copy(a.DMax), copy(a.pStart), copy(a.qStart), copy(a.PStart), copy(a.QStart), copy(a.stationary), copy(a.seasonal), copy(a.ic), copy(a.stepwise), copy(a.parallel), copy(a.numCore))
deepcopy(a::AutoARIMAInput) = AutoARIMAInput(deepcopy(a.d), deepcopy(a.D), deepcopy(a.pMax), deepcopy(a.qMax), deepcopy(a.PMax), deepcopy(a.QMax), deepcopy(a.numCoefMax), deepcopy(a.dMax), deepcopy(a.DMax), deepcopy(a.pStart), deepcopy(a.qStart), deepcopy(a.PStart), deepcopy(a.QStart), deepcopy(a.stationary), deepcopy(a.seasonal), deepcopy(a.ic), deepcopy(a.stepwise), deepcopy(a.parallel), deepcopy(a.numCore))
#show
function show(io::IO, a::AutoARIMAInput)
	println(io, "AutoARIMAInput contents:")
	a.d < 0 ? println(io, "    order of integration = not set") : println(io, "    order of integration = " * string(a.d))
	a.D < 0 ? println(io, "    seasonal order of integration = not set") : println(io, "    seasonal order of integration = " * string(a.D))
	println(io, "    pMax = " * string(a.pMax))
	println(io, "    qMax = " * string(a.qMax))
	println(io, "    PMax = " * string(a.PMax))
	println(io, "    QMax = " * string(a.QMax))
	println(io, "    numCoefMax = " * string(a.numCoefMax))
	println(io, "    dMax = " * string(a.dMax))
	println(io, "    DMax = " * string(a.DMax))
	println(io, "    pStart = " * string(a.pStart))
	println(io, "    qStart = " * string(a.qStart))
	println(io, "    PStart = " * string(a.PStart))
	println(io, "    QStart = " * string(a.QStart))
	println(io, "    restrict search to stationary models = " * string(a.stationary))
	println(io, "    include seasonal models in search path = " * string(a.seasonal))
	println(io, "    information criterion = " * a.ic)
	println(io, "    use step-wise search procedure = " * string(a.stepwise))
	println(io, "    use parallel computation (only possible when step-wise is false) = " * string(a.parallel))
	println(io, "    number of cores for parallel use = " * string(a.numCore))
end
show(a::AutoARIMAInput) = show(STDOUT, a)

#----------------------------------------------------------
#TYPE
#	ARIMAModel
#FIELDS
# 	arCoef::Vector{Float64}: Vector of AR coefficients
# 	maCoef::Vector{Float64}: Vector of MA coefficients
# 	sarCoef::Vector{Float64}: Vector of seasonal AR coefficients
# 	smaCoef::Vector{Float64}: Vector of seasonal MA coefficients
# 	intercept::Float64: Model intercept. Equal to NaN if no intercept.
# 	sigma2::Float64: Variance of residuals
# 	covCoef::Matrix{Float64}: Covariance matrix of estimated coefficients (and intercept if applicable)
# 	covCoefName::Vector{ASCIIString}: Names corresponding to the rows/columns of covCoef
# 	d::Int: Order of integration
# 	D::Int: Seasonal order of integration
# 	seasonalPeriod::Int: Seasonal period.
# 	logLik::Float64: Log-likelihood at maximum as reported by R in estimation step.
# 	aic::Float64: aic at maximum as reported by R in estimation step.
# 	resid::Vector{Float64}: Vector of residuals from estimation step.
# 	optimCode::Int: Optimisation output code as reported by R. 0=normal convergence obtained, other values imply something wrong (see R docs for more detail)
# 	rModelName::ASCIIString: Variable name assigned to model output from R estimation step in embedded R workspace. This field is used to find the R output from estimation step and feed it into R predict.Arima (if the user calls Julia forecast function)
#PURPOSE
#	This type stores a complete ARIMA specification.
#CONSTRUCTORS
#	ARIMAModel(): Default values
#KEYWORD CONSTRUCTOR
#	ARIMAModel(; arCoef::Vector{Float64}=Array(Float64, 0), maCoef::Vector{Float64}=Array(Float64, 0), sarCoef::Vector{Float64}=Array(Float64, 0), smaCoef::Vector{Float64}=Array(Float64, 0), intercept::Float64=0.0, sigma2::Float64=1.0, covCoef::Matrix{Float64}==Array(Float64, 0, 0), covCoefName::Vector{ASCIIString}=Array(ASCIIString, 0), d::Int=0, D::Int=0, seasonalPeriod::Int=1, logLik::Float64=0.0, aic::Float64=0.0, resid::Vector{Float64}=Array(Float64, 0), optimCode::Int=0, rModelName::ASCIIString="")
#METHODS
#	copy(a::ARIMAModel): Standard copy
#	deepcopy(a::ARIMAModel): Standard deepcopy
#	show(io::IO, a::ARIMAModel): Standard show
#	show(a::ARIMAModel): Show defaulting to STDOUT
#	coef(a::ARIMAModel): Returns all coefficients (and intercept if applicable) in one vector ordered [arCoef, maCoef, sarCoef, smaCoef, intecept (if applicable)]
#	vcov(a::ARIMAModel): Return the variance covariance matrix of coefficients
#	coefVar(a::ARIMAModel): Return a 2x? matrix containing all coefficients in first row (ie output of coef functions) and the corresponding variances of these estimators in the second row.
#NOTES
#	Typically this type is created when the user calls the estimate function on an ARIMAInput or AutoARIMAInput. However, the user may also want to create a specific ARIMAModel to feed into the simulate function.
#	As well as serving as output from estimate step, this type serves as input to simulate function, although obviously not all the fields will be relevant.
#----------------------------------------------------------
type ARIMAModel
	arCoef::Vector{Float64}
	maCoef::Vector{Float64}
	sarCoef::Vector{Float64}
	smaCoef::Vector{Float64}
	intercept::Float64
	sigma2::Float64
	covCoef::Matrix{Float64}
	covCoefName::Vector{ASCIIString}
	d::Int
	D::Int
	seasonalPeriod::Int
	logLik::Float64
	aic::Float64
	resid::Vector{Float64}
	optimCode::Int
	rModelName::ASCIIString
	function ARIMAModel(arCoef::Vector{Float64}, maCoef::Vector{Float64}, sarCoef::Vector{Float64}, smaCoef::Vector{Float64}, intercept::Float64, sigma2::Float64, covCoef::Matrix{Float64}, covCoefName::Vector{ASCIIString}, d::Int, D::Int, seasonalPeriod::Int, logLik::Float64, aic::Float64, resid::Vector{Float64}, optimCode::Int, rModelName::ASCIIString)
		(d < 0 || D < 0) && error("d and D must be non-negative")
		sigma2 < 0 && error("sigma2 must be non-negative")
		size(covCoef, 1) != size(covCoef, 2) && error("Coefficient covariance matrix is non-square")
		numCoefTemp = length(arCoef) + length(maCoef) + length(sarCoef) + length(smaCoef) + convert(Int, !(isnan(intercept)))
		size(covCoef, 1) != numCoefTemp && error("Logic fail. Number of coefficients is inconsistent with size of coefficient covariance matrix")
		length(covCoefName) != numCoefTemp && error("Logic fail. Number of coefficients is inconsistent with number of coefficient names")
		seasonalPeriod < 1 && error("seasonalPeriod must be strictly positive")
		funcOut = new(arCoef, maCoef, sarCoef, smaCoef, intercept, sigma2, covCoef, covCoefName, d, D, seasonalPeriod, logLik, aic, resid, optimCode, rModelName)
		return(funcOut)
	end
end
#--------------------
#Constructors
#--------------------
ARIMAModel() = ARIMAModel(Array(Float64, 0), Array(Float64, 0), Array(Float64, 0), Array(Float64, 0), 0.0, 1.0, Array(Float64, 0, 0), Array(ASCIIString, 0), 0, 0, 1, 0.0, 0.0, Array(Float64, 0), 0, "")
#ARIMAInput constructor for converting an ARIMAModel to ARIMAInput (necessary for code to be here after ARIMAModel type definition)
ARIMAInput(a::ARIMAModel) = ARIMAInput(length(a.arCoef), a.d, length(a.maCoef), length(a.sarCoef), a.D, length(a.smaCoef), a.seasonalPeriod, !(isnan(a.intercept)), true, Array(Float64, 0), Array(Float64, 0), "CSS-ML")
#--------------------
#Keyword constructor
#--------------------
function ARIMAModel(; arCoef::Vector{Float64}=Array(Float64, 0), maCoef::Vector{Float64}=Array(Float64, 0), sarCoef::Vector{Float64}=Array(Float64, 0), smaCoef::Vector{Float64}=Array(Float64, 0), intercept::Float64=0.0, sigma2::Float64=1.0, covCoef::Matrix{Float64}=Array(Float64, 0, 0), covCoefName::Vector{ASCIIString}=Array(ASCIIString, 0), d::Int=0, D::Int=0, seasonalPeriod::Int=1, logLik::Float64=0.0, aic::Float64=0.0, resid::Vector{Float64}=Array(Float64, 0), optimCode::Int=0, rModelName::ASCIIString="")
	numCoefTemp = length(arCoef) + length(maCoef) + length(sarCoef) + length(smaCoef) + convert(Int, !(isnan(intercept)))
	if size(covCoef, 1) == 0 && size(covCoef, 2) == 0
		covCoef = zeros(Float64, numCoefTemp, numCoefTemp)
	end
	if length(covCoefName) == 0
		covCoefName = Array(ASCIIString, numCoefTemp)
	end
	return(ARIMAModel(arCoef, maCoef, sarCoef, smaCoef, intercept, sigma2, covCoef, covCoefName, d, D, seasonalPeriod, logLik, aic, resid, optimCode, rModelName))
end
#--------------------
#Methods
#--------------------
#copy, deepcopy
copy(a::ARIMAModel) = ARIMAModel(copy(a.arCoef), copy(a.maCoef), copy(a.sarCoef), copy(a.smaCoef), copy(a.intercept), copy(a.sigma2), copy(a.covCoef), copy(a.d), copy(a.D), copy(a.seasonalPeriod), copy(a.logLik), copy(a.aic), copy(a.resid), copy(a.optimCode), copy(a.rModelName))
deepcopy(a::ARIMAModel) = ARIMAModel(deepcopy(a.arCoef), deepcopy(a.maCoef), deepcopy(a.sarCoef), deepcopy(a.smaCoef), deepcopy(a.intercept), deepcopy(a.sigma2), deepcopy(a.covCoef), deepcopy(a.d), deepcopy(a.D), deepcopy(a.seasonalPeriod), deepcopy(a.logLik), deepcopy(a.aic), deepcopy(a.resid), deepcopy(a.optimCode), deepcopy(a.rModelName))
#show
function show(io::IO, a::ARIMAModel)
	println("ARIMAModel contents:")
	println("    AR coefficients = " * string(a.arCoef))
	println("    MA coefficients = " * string(a.maCoef))
	println("    SAR coefficients = " * string(a.sarCoef))
	println("    SMA coefficients = " * string(a.smaCoef))
	println("    intercept = " * string(a.intercept))
	println("    residual variance = " * string(a.sigma2))
	println("    order of integration = " * string(a.d))
	println("    seasonal order of integration = " * string(a.D))
	println("    seasonal period = " * string(a.seasonalPeriod))
	println("    log-likelihood at maximum = " * string(a.logLik))
	println("    aic at maximum = " * string(a.aic))
	println("    optimisation output code = " * string(a.optimCode))
	println("    R model name = " * string(a.rModelName))
end
show(a::ARIMAModel) = show(STDOUT, a)
#coef
coef(a::ARIMAModel) = isnan(a.intercept) ? [a.arCoef, a.maCoef, a.sarCoef, a.smaCoef] : [a.arCoef, a.maCoef, a.sarCoef, a.smaCoef, [a.intercept]]
#vcov (to match R function)
vcov(a::ARIMAModel) = a.covCoef
#coefVar (get coefficients and their variances)
coefVar(a::ARIMAModel) = cat(2, coef(a), diag(a.covCoef))'

#----------------------------------------------------------
#FUNCTION
#	estimate
#	estimateARIMA
#	estimateAutoARIMA
#METHODS
#KEYWORD METHODS
#OUTPUT
#	ARIMAModel containing the estimated ARIMA coefficients and other variables of interest.
#PURPOSE
#	The purpose of this function is to estimate the parameters of the input model and return a complete specification of the model
#NOTES
#	If I write any custom optimization routines, make them their own type, and then add a method to estimate#NON-EXPORTED FUNCTIONS
#----------------------------------------------------------
#estimate method that uses arima functionality in R (distinguished via method being an ASCIIString)
function estimate{T<:Number}(x::Vector{T}, a::ARIMAInput, method::ASCIIString)
	error("I have not coded up this method of estimate yet")
end
#estimate method that uses Julia Optim package (distinguished via method being a Symbol)
function estimate{T<:Number}(x::Vector{T}, a::ARIMAInput, method::Symbol)
	error("I have not coded up this method of estimate yet")
end
#POTENTIALLY include another method here that allows optimisation via minimisation of arbitrary function of model residuals
#estimate wrappers (ARIMAInput)
estimate{T<:Number}(x::Vector{T}, a::ARIMAInput) = estimate(x, a, a.method)
estimate{T<:Number}(a::ARIMAInput, x::Vector{T}) = estimate(x, a)
#Keyword versions
estimateARIMA{T<:Number}(x::Vector{T}; p::Int=0, d::Int=0, q::Int=0, P::Int=0, D::Int=0, Q::Int=0, seasonalPeriod::Int=1, includeIntercept::Bool=true, transformParams::Bool=true, fixedParams::Vector{Float64}=Array(Float64, 0), initParams::Vector{Float64}=Array(Float64, 0), method::ASCIIString="CSS-ML") = estimate(x, ARIMAInput(p=p, d=d, q=q, P=P, D=D, Q=Q, seasonalPeriod=seasonalPeriod, includeIntercept=includeIntercept, transformParams=transformParams, fixedParams=fixedParams, initParams=initParams, method=method))
estimateAutoARIMA{T<:Number}(x::Vector{T}; d::Int=-1, D::Int=-1, pMax::Int=5, qMax::Int=5, PMax::Int=2, QMax::Int=2, numCoefMax::Int=5, dMax::Int=2, DMax::Int=1, pStart::Int=2, qStart::Int=2, PStart::Int=1, QStart::Int=5, stationary::Bool=false, seasonal::Bool=false, ic::ASCIIString="aic", stepwise::Bool=true, parallel::Bool=false, numCore::Int=2) = estimate(x, AutoARIMAInput(d=d, D=D, pMax=pMax, qMax=qMax, PMax=PMax, QMax=QMax, numCoefMax=numCoefMax, dMax=dMax, DMax=DMax, pStart=pStart, qStart=qStart, PStart=PStart, QStart=QStart, stationary=stationary, seasonal=seasonal, ic=ic, stepwise=stepwise, parallel=parallel, numCore=numCore))
#Estimate ARIMA model with automatic order detection
function estimate{T<:Number}(x::Vector{T}, a::AutoARIMAInput)
	error("I have not coded up this method of estimate yet")
end
estimate{T<:Number}(a::AutoARIMAInput, x::Vector{T}) = estimate(x, a)

#----------------------------------------------------------
#FUNCTION
#	forecast
#INPUT
#	forecast{T<:Number}(x::Vector{T}, a::ARIMAModel, numStep::Int=1): Forecast numStep ahead using ARIMAModel and data in x. This method is implemented entirely in Julia, but is not able to handle seasonal ARIMA models or bounds computations.
#	forecast{T<:Number}(a::ARIMAModel, x::Vector{T}, numStep::Int=1): Wrapper so that order of inputs doesn't matter.
#OUTPUT
#	Output is tuple (pointFore::Vector{Float64}, boundFore::Matrix{Float64}), where pointFore provides the point forecast up to numStep ahead, and boundFore provides 95% lower and upper bounds for the forecasts (lower bound in first column, upper bound in second column)
#PURPOSE
#	The purpose of this function is to provide forecasts for the input dataset and ARIMA model
#NOTES
#	When d > 0 (order of integration greater than 0), the methods that do not wrap R functions use the non-exported function arima_statdiff, arima_statdiffinv, arima_statdiffinv_one. These functions difference or inverse difference data up to the specified order.
#NON-EXPORTED FUNCTIONS
#	forecast_arma{T<:Number}(x::Vector{T}, resid::Vector{T}, arCoef::Vector{Float64}, maCoef::Vector{Float64}, intercept::Float64=NaN, numStep::Int=1): This function provides an ARMA point forecast implemented entirely in Julia.
#----------------------------------------------------------
#Method implemented entirely in Julia (note, bounds computations not yet implemented)
function forecast{T<:Number}(x::Vector{T}, a::ARIMAModel, numStep::Int=1)
	if length(a.sarCoef) > 0 || length(a.smaCoef) > 0 || a.D > 0
		error("This method does not currently support seasonal ARIMA models. Try using the function method that does not include underlying data (ARIMA model must be available in embedded R workspace).")
	end
	if a.d == 0
		(pointFore, boundFore) = forecast_arma(x, a.resid, a.arCoef, a.maCoef, a.intercept, numStep)
	else
		maxLag = max(length(a.arCoef), length(a.maCoef)) + a.d
		xS = x[end-maxLag:end]
		(xD, init) = arima_statdiff(xS, a.d)
		(pointFore, boundFore) = forecast_arma(xD, a.resid[end-length(xD)+1:end], a.arCoef, a.maCoef, a.intercept, numStep)
		pointFore = arima_statdiffinv([xD, pointFore], init)
		pointFore = pointFore[end-numStep+1:end]
		#bound calcuations will need to be added here if they are ever implemented in forecast_arma
	end
	return(pointFore, boundFore)
end
#Wrapper to make order of input immaterial
forecast{T<:Number}(a::ARIMAModel, x::Vector{T}, numStep::Int=1) = forecast(x, a, numStep)
#forecast_arma
function forecast_arma{T<:Number}(x::Vector{T}, resid::Vector{T}, arCoef::Vector{Float64}, maCoef::Vector{Float64}, intercept::Float64=NaN, numStep::Int=1)
	length(x) != length(resid) && error("Input x vector and residual vector must be the same length")
	!(isnan(intercept)) && (x = x - intercept) #de-mean data if we have an intercept
	maxLag = max(length(arCoef), length(maCoef))
	eFore = [resid[end-maxLag+1:end], zeros(Float64, numStep)]
	xFore = [x[end-maxLag+1:end], zeros(Float64, numStep)]
	for k = maxLag+1:maxLag+numStep
		for j = 1:length(arCoef)
			xFore[k] += arCoef[j] * xFore[k-j]
		end
		for j = 1:length(maCoef)
			xFore[k] += maCoef[j] * eFore[k-j]
		end
	end
	pointFore = xFore[end-numStep+1:end]
	!(isnan(intercept)) && (pointFore = pointFore + intercept) #add mean back in if it exists
	boundFore = NaN * Array(Float64, numStep, 2) #bound calculations could be added here
	return(pointFore, boundFore)
end

#----------------------------------------------------------
#FUNCTION
#	simulate
#	simulateARIMA
#METHODS
#	simulate(numObs::Int, a::ARIMAModel, numBurn::Int=-1, resid::Vector{Float64}=[NaN]): Simulate numObs observations using ARIMAModel a, with burn-in length numBurn and optionally use resid as the source of random numbers in simulation.
#	simulate(a::ARIMAModel, numObs::Int, numBurn::Int=-1, resid::Vector{Float64}=[NaN]): Wrapper so order doesn't matter.
#KEYWORD METHOD
#	simulateARIMA(numObs::Int; arCoef::Vector{Float64}=Array(Float64, 0), maCoef::Vector{Float64}=Array(Float64, 0), d::Int=0, intercept::Float64=NaN, sigma2::Float64=1.0, numBurn::Int=-1, resid::Vector{Float64}=[NaN]): Keyword activated simulate method. First input specifiying number of observations to simulate is non-optional. Default keywords result in white noise model.
#OUTPUT
#	Output is Vector{Float64} of data simulated using the specified ARIMA model.
#PURPOSE
#	The purpose of this function is to simulate the specified ARIMA model
#NOTES
#	This function is not currently capable of simulating seasonal ARIMA models, ie non-seasonal only.
#NON-EXPORTED FUNCTION
#	simulate_arma(numObs::Int, arCoef::Vector{Float64}, maCoef::Vector{Float64}, intercept::Float64=0, sigma2::Float64=1, numBurn::Int=-1, resid::Vector{Float64}=[NaN]): Simulates an ARMA model using specified inputs. simulate calls this function for simulating the ARMA part and uses arima_statdiff and arima_statdiffinv to take care of the order of integration
#----------------------------------------------------------
function simulate(numObs::Int, a::ARIMAModel, numBurn::Int=-1, resid::Vector{Float64}=[NaN])
	if length(a.sarCoef) > 0 || length(a.smaCoef) > 0 || a.D > 0
		error("Simluation of seasonal ARIMA models has not been implemented yet")
	end
	isnan(a.intercept) ? intercept = 0.0 : intercept = a.intercept
	xSim = simulate_arma(numObs, a.arCoef, a.maCoef, intercept, a.sigma2, numBurn, resid)
	if a.d > 0
		xSim = arima_statdiffinv(xSim, zeros(Float64, a.d))[a.d+1:end]
	end
	return(xSim)
end
#Wrapper so order doesn't matter
simulate(a::ARIMAModel, numObs::Int, numBurn::Int=-1, resid::Vector{Float64}=[NaN]) = simulate(numObs, a, numBurn, resid)
#Keyword activated function
simulateARIMA(numObs::Int; arCoef::Vector{Float64}=Array(Float64, 0), maCoef::Vector{Float64}=Array(Float64, 0), d::Int=0, intercept::Float64=NaN, sigma2::Float64=1.0, numBurn::Int=-1, resid::Vector{Float64}=[NaN]) = simulate(numObs, ARIMAModel(arCoef=arCoef, maCoef=maCoef, d=d, intercept=intercept, sigma2=sigma2), numBurn, resid)
#non-exported function that simulates an arma
function simulate_arma(numObs::Int, arCoef::Vector{Float64}, maCoef::Vector{Float64}, intercept::Float64=0, sigma2::Float64=1, numBurn::Int=-1, resid::Vector{Float64}=[NaN])
	p = length(arCoef)
	q = length(maCoef)
	if numBurn < 0
		#IMPLEMENT HERE: A better method for choosing the number of observations to burn. Should be straightforward with Polynomials package, e.g. p+q+ceil(6/log(minroots)), where minroots = min(modulus(roots(p))) where p is polynomial [-1, ar coefficients]
		numBurn = p + q + 200
	elseif numBurn < p + q
		error("Input number of burn-in observations must be of at least length p + q")
	end
	length(resid) == 0 && error("Invalid optional input resid")
	isnan(resid[1]) && (resid = sqrt(sigma2) * randn(numObs + numBurn))
	if length(resid) > numObs + numBurn
		#resid = resid[end-numObs-numBurn+1:end]
		resid = resid[1:numObs+numBurn]
	elseif length(resid) < numObs + numBurn
		error("Number of supplied residuals must be >= numObs + numBurn")
	end
	xSim = intercept * ones(Float64, numObs + numBurn)
	for n = max(p, q)+1:length(xSim)
		for k = 1:length(arCoef) #
			xSim[n] += arCoef[k] * xSim[n-k]
		end
		for k = 1:length(maCoef)
			xSim[n] += maCoef[k] * resid[n-k]
		end
		xSim[n] += resid[n]
	end
	return(xSim[end-numObs+1:end])
end

#----------------------------------------------------------
#FUNCTION
#	arima_statdiff
#	arima_statdiffinv_one
#	arima_statdiffinv
#INPUT
#	arima_statdiff{T<:Number}(x::Vector{T}, order::Int=1): x is a vector of numbers to discretely differentiate. order is the order of differentiation. If order=1, this is equivalent to calling diff from Juia base.
#	arima_statdiffinv_one{T<:Number}(x::Vector{T}, init::T): x is a vector of numbers to discretely integrate ONCE. init is the initial value needed for discrete integration
#	arima_statdiffinv{T<:Number}(x::Vector{T}, init::Vector{T}): x is a vector of numbers to discretely integrate length(init) times. init provides the initial values needed for each successive discrete integration, where the first integration of x uses init[1] and the last integration uses init[end]
#OUTPUT
#	arima_statdiff outputs the tuple (xOut, init), where xOut is the differentiated series and init is the first value of each series prior to the last differentiated series, so init[1] = x[1], init[2] = (delta x)[1], ..., and init[end] = (delta^(order-1) x)[1]
#	arima_statdiffinv_one and arima_statdiffinv both output Vector{Float64} storing the discretely integrated series.
#PURPOSE
#	These three function provide discrete differentiation and discrete integration.
#NOTES
#	These functions are not exported.
#	More efficient versions of arima_statdiff and arima_statdiffinv will almost certainly eventually become a part of Julia Base. The functions provided here can be used in the meantime.
#----------------------------------------------------------
#arima_statdiff
function arima_statdiff{T<:Number}(x::Vector{T}, order::Int=1)
	#WARNING: Temporary array allocated for each order of differencing. Highly inefficient when order is large.
	order < 1 && error("order of statistical differencing must be at least 1")
	length(x) < order + 1 && error("input vector x too short for specified order of differencing")
	xOut = diff(x)
	init = Array(T, order)
	init[1] = x[1]
	for n = 2:order
		init[n] = xOut[1]
		xOut = diff(xOut)
	end
	return(xOut, init)
end
#arima_statdiffinv_one
function arima_statdiffinv_one{T<:Number}(x::Vector{T}, init::T)
	length(x) < 1 && error("input vector x must contain at least one element")
	xOut = Array(T, length(x) + 1)
	xOut[1] = init
	for n = 2:length(xOut)
		xOut[n] = xOut[n-1] + x[n-1]
	end
	return(xOut)
end
#arima_statdiffinv
function arima_statdiffinv{T<:Number}(x::Vector{T}, init::Vector{T})
	#WARNING: Temporary array allocated for each order of integration. Highly inefficient when order is large.
	length(init) < 1 && error("order of statistical differencing must be at least 1")
	xOut = arima_statdiffinv_one(x, init[end])
	for n = 2:length(init)
		xOut = arima_statdiffinv_one(xOut, init[end-n+1])
	end
	return(xOut)
end

end # module


text:
@generated function Base.zero{U<:AbstractUnum}(::Type{U})
  c = unumConstants(U)
  :($(c.zero))
end

Base.one{U<:AbstractUnum}(::Type{U}) = convert(U, 1)

"exponent of the unum"
@generated function Base.exponent{U<:AbstractUnum}(u::U)
  c = unumConstants(U)
  :(reinterpret($(c.INT), u & $(c.emask)) >> $(c.fpos))
end

"significand (fraction) of an unum"
@generated function Base.significand{U<:AbstractUnum}(u::U)
  c = unumConstants(U)
  :(reinterpret($(c.INT), u & $(c.fmask)) >> $(c.ubitpos))
end

"Number of bits in the unum's exponent"
@generated function esize{U<:AbstractUnum}(u::U)
  c = unumConstants(U)
  :(1 + reinterpret($(c.INT), u & $(c.esizemask)) >> $(c.fsizepos))
end

"Number of bits in the unum's significand (fraction)"
@generated function fsize{U<:AbstractUnum}(u::U)
  c = unumConstants(U)
  :(1 + reinterpret($(c.INT), u & $(c.fsizemask)))
end

@generated function isnegative{U<:AbstractUnum}(u::U)
  c = unumConstants(U)
  :(reinterpret($(c.UINT), u & $(c.signbitmask)) != zero($(c.UINT)))
end

ispositive(u::AbstractUnum) = !isnegative(u)

@generated function isinexact{U<:AbstractUnum}(u::U)
  c = unumConstants(U)
  :(reinterpret($(c.UINT), u & $(c.ubitmask)) != zero($(c.UINT)))
end

isexact(u::AbstractUnum) = !isinexact(u)

@generated function typemin{U<:AbstractUnum}(::Type{U})
  c = unumConstants(U)
  :(c.neginf)
end

@generated function typemax{U<:AbstractUnum}(::Type{U})
  c = unumConstants(U)
  :(c.posinf)
end

@generated function maxubits{U<:AbstractUnum}(::Type{U})
  c = unumConstants(U)
  :(2 + $(c.esizesize) + $(c.fsizesize) + 2^$(c.esizesize) + 2^$(c.fsizesize))
end

# for func in [:isnan, :isnull, :isinf, :isfinite, :isposinf, :isneginf,
#               :isnegative, :iszero, :ispositive, :-, :isinteger, :issubnormal]

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

# @generated function +{U<:AbstractUnum}(u1::U, u2::U)

# end

# for op in [:+, :-, :*, :/, :(==), :<, :>, :<=, :>=, :isless, :isequal]

ext: jl


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://twitter.com/physics303/status/1414522243216859143?s=21

# %%
f(A) = sortperm(sort.(collect(eachrow(A))); rev=true)

function g!(A, ix, tmp)
    for i in axes(A, 1) tmp[i] .= @view A[i, :] end
    sort!.(tmp)
    sortperm!(ix, tmp; rev=true)
end

function maketmp(A)
    ix = Vector{Int}(undef, size(A, 1))
    tmp = [Vector{eltype(A)}(undef, size(A, 2)) for _ in 1:size(A, 1)]
    ix, tmp
end

B = [
    2 4 4 4 5 9
    1 2 2 3 5 9
    1 2 3 3 3 9
    1 2 2 5 6 9
    1 3 4 4 4 9
]

ix, tmp = maketmp(B)
@show ix
@show tmp
@show f(B) == g!(B, ix, tmp)
@time f(B)
@time g!(B, ix, tmp)

# %%
using BenchmarkTools
@btime f($B)
@btime g!($B, $ix, $tmp)

# %%

max_stars_repo_path: 0010/sortperm! example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0010/sortperm! example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0010/sortperm! example.jl

max_forks_repo_name: genkuroki/public


text:
"""
    InputData(name::Symbol, option::Symbol, data::Any)

Represents a more structured block of input data.
e.g. `InputData(:k_points, :automatic, (6,6,6,1,1,1))`
would be translated for a QE calculation into
```
K_POINTS(automatic)
6 6 6 1 1 1
```
"""
mutable struct InputData
    name   :: Symbol
    option :: Symbol
    data   :: Any
end

StructTypes.StructType(::Type{InputData}) = StructTypes.Struct()

abstract type Package end
struct NoPackage <: Package end
struct Wannier90 <: Package end
struct QE <: Package end
struct Abinit <: Package end
struct Elk <: Package end
StructTypes.StructType(::Type{<:Package}) = StructTypes.Struct()

"""
    Calculation{P<:Package}(name    ::String;
                            flags   ::AbstractDict = Dict{Symbol, Any}(),
                            data    ::Vector{InputData} = InputData[],
                            exec    ::Exec,
                            run     ::Bool = true,
                            infile  ::String = P == Wannier90 ? name * ".win" : name * ".in",
                            outfile ::String = name * ".out")

The representation of a *DFT* calculation of package `P`,
holding the `flags` that will be written to the `infile`,
the executable `exec` and the output written by the calculation to the `outfile`.
It essentially represents a line in a job script similar to `exec < infile.in > outfile.out`.
`outdata` stores the parsed calculation output after it was read at least once.
The `run` field indicates whether the calculation should be actually performed,
e.g. if `run=false` the corresponding line will be commented out in the job script.

    Calculation{P<:Package}(name::AbstractString, flags::Pair{Symbol, Any}...; kwargs...)

Create a [`Calculation`](@ref) from `name` and `flags`, other `kwargs...` will be passed to the constructor.

    Calculation(template::Calculation, name::AbstractString, flags::Pair{Symbol, Any}...; excs=deepcopy(template.exec), run=true, data=nothing)

Creates a new [`Calculation`](@ref) from the `template`, setting the `flags` of the newly created one to the specified ones.
"""
@with_kw_noshow mutable struct Calculation{P<:Package}
    name::String
    flags::Dict{Symbol,Any} = Dict{Symbol,Any}()
    data::Vector{InputData} = InputData[]
    exec::Exec
    run::Bool = true
    infile::String = P == Wannier90 ? name * ".win" : name * ".in"
    outfile::String = P == Wannier90 ? name * ".wout" : name * ".out"
    function Calculation{P}(name, flags, data, exec, run, infile,
                            outfile) where {P<:Package}
        out = new{P}(name, Dict{Symbol,Any}(), data, exec, run, infile,
                     outfile)
        set_flags!(out, flags...; print = false)
        for (f, v) in flags
            if !hasflag(out, f)
                @warn "Flag $f was not valid for calculation $name."
            end
        end
        return out
    end
end
function Calculation{P}(name, flags, data, exec, run) where {P<:Package}
    return Calculation{P}(name, flags, data, exec, run, Dict{Symbol,Any}(),
                          P == Wannier90 ? name * ".win" : name * ".in",
                          P == Wannier90 ? name * ".wout" : name * ".out")
end

function Calculation{P}(name, flags...; kwargs...) where {P<:Package}
    return Calculation{P}(; name = name, flags = flags, kwargs...)
end

function Calculation(template::Calculation, name, newflags...;
                     excs = deepcopy(template.exec), run  = true, data = nothing)
    newflags = Dict(newflags...)

    calculation       = deepcopy(template)
    calculation.name  = name
    calculation.exec = excs
    calculation.run   = run
    set_flags!(calculation, newflags...; print = false)

    if data !== nothing
        for (name, (option, dat)) in data
            d = getfirst(x -> x.name == name, calculation.data)
            if d !== nothing
                d.data = dat
                d.option = option
            else
                push!(calculation.data, InputData(name, option, dat))
            end
        end
    end
    return calculation
end

# Calculation() = Calculation{NoPackage}(package=NoPackage())
StructTypes.StructType(::Type{<:Calculation}) = StructTypes.Mutable()

# Interface Functions
isbands(c::Calculation)    = false
isnscf(c::Calculation)     = false
isscf(c::Calculation)      = false
isvcrelax(c::Calculation)  = false
isrelax(c::Calculation)    = false
ismagnetic(c::Calculation) = false
issoc(c::Calculation)      = false
outfiles(c::Calculation)   = [c.outfile]
ispw(c::Calculation)       = false
isprojwfc(c::Calculation)  = false

"""
    set_name!(c::Calculation, name::AbstractString)

Sets `calculation.name`, and `calculation.infile` and `calculation.outfile` to conform
with the new `name`.
"""
function set_name!(c::Calculation, name::AbstractString; print = true)
    c.name = name
    c.infile = name * splitext(c.infile)[2]
    c.outfile = name * splitext(c.outfile)[2]
    print && @info "\nname = $name\ninfile = $(c.infile)\noutfile = $(c.outfile)"
    return name
end

"""
    data(calculation::Calculation, n::Symbol)

The former returns `calculation.data`, the later -- the `InputData` with name `n`.
"""
data(calculation::Calculation, n::Symbol) = getfirst(x -> x.name == n, calculation.data)

#
# Directory interface
#
DFC.Utils.searchdir(i::Calculation, glob) = searchdir(i.dir, glob)

Base.eltype(::Calculation{T}) where {T} = T

#
# Flag Dict interace
#

"""
    getindex(c::Calculation, n::Symbol)

Returns the flag with given symbol.

    getindex(job::Job, flag::Symbol)

Searches through the job's calculations for the requested flag.
A `Dict` will be returned with calculation names as the keys and the flags as values.
"""
Base.getindex(c::Calculation, n::Symbol) = haskey(c, n) ? c.flags[n] : throw(KeyError(n))

Base.haskey(c::Calculation, n::Symbol) = haskey(c.flags, n)
Base.get(c::Calculation, args...) = get(c.flags, args...)
Base.pop!(c::Calculation, args...) = pop!(c.flags, args...)

hasflag(c::Calculation, s::Symbol) = haskey(c.flags, s)
hasflag(c::Calculation, s) = false

"""
    setindex!(c::Calculation, value, flag::Symbol)

Sets flags.

    setindex!(job::Job, value, flag::Symbol)

Set `flag` in all the appropriate calculations to the `value`.
"""
Base.setindex!(c::Calculation, dat, key) = set_flags!(c, key => dat)

"""
    set_flags!(c::Calculation, flags::Pair{Symbol, Any}...; print=true)

Sets multiple flags in one go. Flag validity and type are verified.

    set_flags!(job::Job, calculations::Vector{<:Calculation}, flags::Pair{Symbol,<:Any}...; print=true)
    set_flags!(job::Job, flags::Pair{Symbol,<:Any}...; print=true)

Sets the flags in the names to the flags specified.
This only happens if the specified flags are valid for the names.

The values that are supplied will be checked whether they are valid.
"""
function set_flags!(c::Calculation{T}, flags...; print = true) where {T}
    found_keys = Symbol[]
    for (flag, value) in flags
        flag_type = flagtype(c, flag)
        if flag_type !== Nothing
            !(flag in found_keys) && push!(found_keys, flag)
            try
                if isa(value, AbstractVector{<:AbstractVector}) &&
                   flag_type <: AbstractVector
                    value = [convert.(eltype(flag_type), v) for v in value]
                elseif flag == :starting_ns_eigenvalue
                    if length(size(value)) == 3
                        value = convert.(Float32, value)
                    else
                        nat = div(length(value), 7*4)
                        value = convert.(Float32, reshape(value, (7,4,nat)))
                    end
                elseif flag == :Hubbard_occupations
                    if length(size(value)) == 4
                        value = convert.(Float32, value)
                    else
                        nat = div(length(value), 7*7*4)
                        value = convert.(Float32, reshape(value, (7, 7,4,nat)))
                    end
                else
                    value = convert(flag_type, value)
                end
            catch
                print &&
                    (@warn "Filename '$(c.name)':\n  Could not convert '$value' into '$flag_type'.\n    Flag '$flag' not set.\n")
                continue
            end
            old_data = haskey(c.flags, flag) ? c.flags[flag] : ""
            c.flags[flag] = value
            print && (@info "$(c.name): -> $flag:\n      $old_data set to: $value\n")
        else
            print &&
                @warn "Flag $flag was ignored since it could not be found in the allowed flags for calculation $(c.name)"
        end
    end
    return found_keys, c
end

"Runs through all the set flags and checks if they are allowed and set to the correct value"
function convert_flags!(calculation::Calculation)
    for (flag, value) in calculation.flags
        flagtype_ = flagtype(calculation, flag)
        if flagtype_ == Nothing
            @warn "Flag $flag was not found in allowed flags for exec $(calculation.exec). Removing flag."
            rm_flags!(calculation, flag)
            continue
        end
        if !(isa(value, flagtype_) || eltype(value) <: flagtype_)
            try
                if isbitstype(eltype(value))
                    if length(value) > 1
                        calculation[flag] = convert(flagtype_, value)
                    else
                        calculation[flag] = convert(eltype(flagtype_), value)
                    end
                else
                    calculation[flag] = convert.(flagtype_, value)
                end
            catch
                error("Input $(calculation.name): Could not convert :$flag of value $value to the correct type ($flagtype_), please set it to the correct type.")
            end
        end
    end
end

function Base.:(==)(d1::InputData, d2::InputData)
    return all(x -> getfield(d1, x) == getfield(d2, x), fieldnames(InputData))
end

function Base.:(==)(i1::Calculation, i2::Calculation)
    return all(x -> x in (:outdata, :dir, :run) ? true : getfield(i1, x) == getfield(i2, x),
               fieldnames(Calculation))
end

ψ_cutoff_flag(::Calculation) = nothing
ρ_cutoff_flag(::Calculation) = nothing

function set_cutoffs!(c::Calculation, ecutwfc, ecutrho)
    return set_flags!(c, ψ_cutoff_flag(c) => ecutwfc, ρ_cutoff_flag(c) => ecutrho)
end

"""
    set_kpoints!(calculation::Calculation{QE}, k_grid::NTuple{3, Int}; print=true)
    set_kpoints!(calculation::Calculation{QE}, k_grid::NTuple{6, Int}; print=true)
    set_kpoints!(calculation::Calculation{QE}, k_grid::Vector{<:NTuple{4}}; print=true, k_option=:crystal_b)

Convenience function to set the `:k_points` data block of `calculation`.
The three different methods are targeted at `nscf`, `scf` or `vcrelax`,
and `bands` calculations, respectively.
For the `nscf` version an explicit list of `k_points` will be generated.

    set_kpoints!(calculation::Calculation{Wannier90}, k_grid::NTuple{3, Int})

Similar to the `nscf` targeted function in the sense that it will generate
an explicit list of `k_points`, adhering to the same rules as for the `nscf`.
The `mp_grid` flag will also automatically be set.
"""
function set_kpoints!(::Calculation{P}, args...; kwargs...) where {P}
    @error "set_kpoints! not implemented for package $P."
end

#-------- Generating new Calculations ---------- #

function calculation_from_kpoints(template::Calculation, newname, kpoints, newflags...)
    newcalc = Calculation(deepcopy(template); name = newname)
    set_flags!(newcalc, newflags...; print = false)
    set_name!(newcalc, newname)
    set_kpoints!(newcalc, kpoints; print = false)
    return newcalc
end

function gencalc_nscf(::Calculation{P}, args...) where {P}
    @error "gencalc_nscf is not implemented for package $P."
end

function gencalc_scf(::Calculation{P}, args...) where {P}
    @error "gencalc_scf is not implemented for package $P."
end

function gencalc_projwfc(::Calculation{P}, args...) where {P}
    @error "gencalc_projwfc is not implemented for package $P."
end

function gencalc_wan(::Calculation{P}, args...) where {P}
    @error "gencalc_wan is not implemented for package $P."
end

function gencalc_vcrelax(::Calculation{P}, args...) where {P}
    @error "gencalc_vcrelax is not implemented for package $P."
end

function gencalc_bands(::Calculation{P}, args...) where {P}
    @error "gencalc_bands is not implemented for package $P."
end

function isconverged(::Calculation{P}) where {P}
    @error "isconverged is not implemented for package $P."
end

"""
    kgrid(na, nb, nc, calculation)

Returns an array of k-grid points that are equally spaced, calculation can be either `:wan` or `:nscf`, the returned grids are appropriate as calculations for wannier90 or an nscf calculation respectively.
"""
kgrid(na, nb, nc, ::Calculation{T}) where {T} = kgrid(na, nb, nc, T)

function sanitize_flags!(cs::Vector{<:Calculation}, str::Structure, name, outdir)
    if any(x -> eltype(x) == Wannier90, cs) && any(isnscf, cs)
        nscfcalc = getfirst(isnscf, cs)
        wancalc = getfirst(x -> eltype(x) == Wannier90, cs)
        if eltype(nscfcalc) == Elk
            wancalc.flags[:num_bands]         = length(nscfcalc[:wann_bands])
            nscfcalc.flags[:wann_projections] = Structures.projections_string.(unique(filter(x -> !isempty(x.projections), str.atoms)))
            nscfcalc.flags[:elk2wan_tasks]    = ["602", "604"]
            nscfcalc.flags[:wann_seedname]    = Symbol(name)
            if get(wancalc, :wannier_plot, false)
                push!(nscfcalc[:elk2wan_tasks], "605")
            end
        end
    end

    u_ats = unique(str.atoms)
    isnc = Structures.isnoncolin(str)
    flags_to_set = []
    if any(x -> x.dftu.U != 0 ||
                    x.dftu.J0 != 0.0 ||
                    sum(x.dftu.J) != 0 ||
                    sum(x.dftu.α) != 0, u_ats)
        Jmap = map(x -> copy(x.dftu.J), u_ats)
        Jdim = maximum(length.(Jmap))
        Jarr = zeros(Jdim, length(u_ats))
        for (i, J) in enumerate(Jmap)
            diff = Jdim - length(J)
            if diff > 0
                for d in 1:diff
                    push!(J, zero(eltype(J)))
                end
            end
            Jarr[:, i] .= J
        end
        append!(flags_to_set,
                [:Hubbard_U     => map(x -> x.dftu.U, u_ats),
                 :Hubbard_alpha => map(x -> x.dftu.α, u_ats),
                 :Hubbard_beta  => map(x -> x.dftu.β, u_ats), :Hubbard_J     => Jarr,
                 :Hubbard_J0    => map(x -> x.dftu.J0, u_ats)])
    end
    if !isempty(flags_to_set) || any(x -> hasflag(x, :Hubbard_parameters), cs)
        push!(flags_to_set, :lda_plus_u => true)
        if isnc
            push!(flags_to_set, :lda_plus_u_kind => 1)
        end
    end
    if !isempty(flags_to_set)
        for c in cs
            set_flags!(c, flags_to_set...; print = false)
        end
    else
        for c in cs
            for f in
                (:lda_plus_u, :lda_plus_u_kind, :Hubbard_U, :Hubbard_alpha, :Hubbard_beta,
                 :Hubbard_J, :Hubbard_J0, :U_projection_type)
                pop!(c, f, nothing)
            end
        end
    end

    flags_to_set = []
    mags = map(x -> x.magnetization, u_ats)
    starts = Float64[]
    θs = Float64[]
    ϕs = Float64[]
    ismagcalc = isnc ? true : Structures.ismagnetic(str)
    if (ismagcalc && isnc) || any(x -> get(x, :noncolin, false), cs)
        for m in mags
            tm = normalize(m)
            if norm(m) == 0
                push!.((starts, θs, ϕs), 0.0)
            else
                θ = acos(tm[3]) * 180 / π
                ϕ = atan(tm[2], tm[1]) * 180 / π
                start = norm(m)
                push!(θs, θ)
                push!(ϕs, ϕ)
                push!(starts, start)
            end
        end
        push!(flags_to_set, :noncolin => true)
    elseif ismagcalc
        for m in mags
            push!.((θs, ϕs), 0.0)
            if norm(m) == 0
                push!(starts, 0)
            else
                push!(starts, sign(sum(m)) * norm(m))
            end
        end
    end
    append!(flags_to_set, [:starting_magnetization => starts, :angle1 => θs, :angle2 => ϕs])
    ismagcalc && !isnc && push!(flags_to_set, :nspin => 2)
    for c in cs
        set_flags!(c, :prefix => "$name", :outdir => "$outdir"; print = false)
        if ispw(c)
            set_flags!(c, flags_to_set...; print = false)
            if isnc
                pop!(c, :nspin, nothing)
            end
            if isvcrelax(c)
                #this is to make sure &ions and &cell are there in the calculation
                !hasflag(c, :ion_dynamics) &&
                    set_flags!(c, :ion_dynamics => "bfgs"; print = false)
                !hasflag(c, :cell_dynamics) &&
                    set_flags!(c, :cell_dynamics => "bfgs"; print = false)
            end
            #TODO add all the required flags
            @assert hasflag(c, :calculation) "Please set the flag for calculation with name: $(name(c))"

            set_flags!(c, :pseudo_dir => "."; print = false)
        end
        convert_flags!(c)
    end

    return
end

rm_tmp_flags!(::Calculation) = nothing
function rm_tmp_flags!(c::Calculation{QE})
    pop!(c, :prefix, nothing)
    pop!(c, :outdir, nothing)
    return pop!(c, :nspin, nothing)
end


text:
using Test
using Footings
using StaticModules

@testset "TestSetField" begin
    @info "Running Field TestSet ..."

    p = Parameter(name=:p, type=Int, description="This is a parameter.")
    t = Temporary(name=:t, type=Int, description="This is a Temporary.")
    r = Result(name=:r, type=Int, description="This is a Return.")
    @test isa(p, Parameter)
    @test isa(t, Temporary)
    @test isa(r, Result)
end

@testset "TestSetModelBase" begin
    @info "Running Model TestSet ..."

    @staticmodule ns begin
        @doc "Add x and y" ->
        function func_x_y(x, y)
            return x + y
        end

        @doc "Add t and z" ->
        function func_t_z(t, z)
            return t + z
        end
    end

    model1 = Model(name="Test")
    add_parameter!(model1, name=:x, type=Int, description="This is x.")
    add_parameter!(model1, name=:y, type=Int, description="This is y.")
    add_parameter!(model1, name=:z, type=Int, description="This is z.")
    add_temporary!(model1, name=:t, type=Int, description="This is t.")
    add_result!(model1, name=:r, type=Int, description="This is r.")

    add_step!(model1, ns, assign=:t, func=:func_x_y, args=(:x, :y))
    add_step!(model1, ns, assign=:z, func=:func_t_z, args=(:t, :z))

    run_test = build_runner(model1)
    @test run_test(x=1, y=2, z=3) == 6

    model2 = @model "TestMacro" ns begin
        Parameter(name=:x, type=Int, description="This is x.")
        Parameter(name=:y, type=Int, description="This is y.")
        Parameter(name=:z, type=Int, description="This is z.")
        Temporary(name=:t, type=Int, description="This is t.")
        Result(name=:r, type=Int, description="This is r.")

        Steps(:t=>func_x_y(x, y), :r=>func_t_z(t, z))
    end

    run_testmacro = build_runner(model2)
    @test run_testmacro(x=1, y=2, z=3) == 6
end


text:
abstract type AbstractTigerPOMDP{S, O} <: POMDP{S, Symbol, O} end

# Infinite Horizon

@with_kw struct COTigerPOMDP <: AbstractTigerPOMDP{Symbol, Float64}
    r_wait::Float64             = -1.0
    r_listen::Float64           = -2.0
    r_findtiger::Float64        = -100.0
    r_escapetiger::Float64      = 10.0
    p_correct::Float64          = 0.85
    discount::Float64           = 0.95
end

@with_kw struct DOTigerPOMDP <: AbstractTigerPOMDP{Symbol, Symbol}
    r_wait::Float64             = -1.0
    r_listen::Float64           = -2.0
    r_findtiger::Float64        = -100.0
    r_escapetiger::Float64      = 10.0
    p_correct::Float64          = 0.85
    discount::Float64           = 0.95
end

# Finite Horizon

abstract type AbstractTimedTigerPOMDP{O} <: AbstractTigerPOMDP{Tuple{Symbol,Int}, O} end

@with_kw struct TimedCOTigerPOMDP <: AbstractTimedTigerPOMDP{Float64}
    r_wait::Float64             = -1.0
    r_listen::Float64           = -2.0
    r_findtiger::Float64        = -100.0
    r_escapetiger::Float64      = 10.0
    p_correct::Float64          = 0.85
    horizon::Int                = 5
    discount::Float64           = 0.95
end

@with_kw struct TimedDOTigerPOMDP <: AbstractTimedTigerPOMDP{Symbol}
    r_wait::Float64             = -1.0
    r_listen::Float64           = -2.0
    r_findtiger::Float64        = -100.0
    r_escapetiger::Float64      = 10.0
    p_correct::Float64          = 0.85
    horizon::Int                = 5
    discount::Float64           = 0.95
end

const ContinuousTiger = Union{COTigerPOMDP, TimedCOTigerPOMDP}
const DiscreteTiger = Union{DOTigerPOMDP, TimedDOTigerPOMDP}

const stateinds = (left=1, right=2, done=3)
const actioninds = (left=1, right=2, wait=3, listen=4)
const obsinds = (left=1, right=2)
tigerpos(s::Symbol) = s
tigerpos(st::Tuple{Symbol,Int}) = st[1]
stepindex(st::Tuple{Symbol,Int}) = st[2]

horizon(m::AbstractTimedTigerPOMDP) = m.horizon

POMDPs.states(::AbstractTigerPOMDP) = keys(stateinds)
POMDPs.actions(::AbstractTigerPOMDP) = keys(actioninds)
POMDPs.states(m::AbstractTimedTigerPOMDP) = product(keys(stateinds), 0:horizon(m))
POMDPs.observations(::DiscreteTiger) = keys(obsinds)

POMDPs.stateindex(::AbstractTigerPOMDP, s) = stateinds[s]
function POMDPs.stateindex(m::AbstractTimedTigerPOMDP, s)
    LinearIndices((length(stateinds), horizon(m)+1))[stateinds[tigerpos(s)], stepindex(s)+1]
end
POMDPs.actionindex(::AbstractTigerPOMDP, a) = actioninds[a]
POMDPs.obsindex(::DiscreteTiger, o) = obsinds[o]

POMDPs.initialstate_distribution(m::AbstractTigerPOMDP) = POMDPModelTools.Uniform((:left, :right))
function POMDPs.initialstate_distribution(m::AbstractTimedTigerPOMDP)
    POMDPModelTools.Uniform((t, 0) for t in (:left,:right))
end

function POMDPs.transition(m::AbstractTigerPOMDP, s, a)
    if a in (:left, :right)
        return Deterministic(:done)
    else
        return Deterministic(s)
    end
end

function POMDPs.transition(m::AbstractTimedTigerPOMDP, s, a)
    if a in (:left, :right)
        return Deterministic((:done, min(horizon(m), stepindex(s)+1)))
    else
        return Deterministic((tigerpos(s), min(horizon(m), stepindex(s)+1)))
    end
end

POMDPs.isterminal(m::AbstractTigerPOMDP, s) = s == :done
POMDPs.isterminal(m::AbstractTimedTigerPOMDP, s) = stepindex(s) >= horizon(m) || tigerpos(s) == :done

function POMDPs.observation(m::DiscreteTiger, a, sp)
    if a == :listen
        probs = (m.p_correct, 1.0-m.p_correct)
    else # a is a door open or wait
        probs = (0.5, 0.5)
    end
    if tigerpos(sp) == :left
        obs = (:left, :right)
    else
        obs = (:right, :left)
    end
    return SparseCat(obs, probs)
end

struct MultiUniformDistribution{N}
    probs::NTuple{N, Float64}
    dists::NTuple{N, Distributions.Uniform{Float64}}
end

function POMDPs.rand(rng::AbstractRNG, d::MultiUniformDistribution)
    i = rand(rng, Categorical(SVector(d.probs)))
    return rand(rng, d.dists[i])
end

function POMDPs.pdf(d::MultiUniformDistribution, o)
    for i in 1:length(d.probs)
        u = d.dists[i]
        if minimum(u) <= o <= maximum(u)
            return pdf(u, o)*d.probs[i]
        end
    end
    return 0.0
end

# left = Uniform(0, 0.5), right = Uniform(0.5, 1)
function POMDPs.observation(m::ContinuousTiger, a, sp)
    if a == :listen
        probs = (m.p_correct, 1.0-m.p_correct)
    else # a is a door open or wait
        probs = (0.5, 0.5)
    end
    if tigerpos(sp) == :left
        dists = (Distributions.Uniform(0.0, 0.5), Distributions.Uniform(0.5, 1.0))
    else
        dists = (Distributions.Uniform(0.5, 1.0), Distributions.Uniform(0.0, 0.5))
    end
    return MultiUniformDistribution(probs, dists)
end

function POMDPs.reward(m::AbstractTigerPOMDP, s, a)
    if isterminal(m, s)
        return 0.0
    end
    if a == :wait
        return m.r_wait
    elseif a == :listen
        return m.r_listen
    elseif tigerpos(s) == a # a is open
        return m.r_findtiger
    else
        return m.r_escapetiger
    end
end

POMDPs.discount(m::AbstractTigerPOMDP) = m.discount


text:
@testset "Test RQC module" begin
    @testset "Test 2D indexing" begin
        @testset "2D -> 1D Rectangular (5x7)" begin
            x_dim = 5
            y_dim = 7
            Av=Vector(1:x_dim*y_dim)
            A = reshape(Av, (x_dim,y_dim))
            for j in 1:x_dim
                for i in 1:y_dim
                    @test Av[QXZoo.RQC.lin_idx(i, j, x_dim)] == A[j,i]
                end
            end
        end
        @testset "2D -> 1D Square even (4x4)" begin
            x_dim = 4
            y_dim = 4
            Av=Vector(1:x_dim*y_dim)
            A = reshape(Av, (x_dim,y_dim))
            for j in 1:x_dim
                for i in 1:y_dim
                    @test Av[QXZoo.RQC.lin_idx(i, j, x_dim)] == A[j,i]
                end
            end
        end
        @testset "2D -> 1D Square odd (5x5)" begin
            x_dim = 5
            y_dim = 5
            Av=Vector(1:x_dim*y_dim)
            A = reshape(Av, (x_dim,y_dim))
            for j in 1:x_dim
                for i in 1:y_dim
                    @test Av[QXZoo.RQC.lin_idx(i, j, x_dim)] == A[j,i]
                end
            end
        end

        @testset "1D -> 2D Rectangular (5x7)" begin
            x_dim = 5
            y_dim = 7
            Av=Vector(1:x_dim*y_dim)
            A = reshape(Av, (x_dim,y_dim))
            for i in 1:x_dim*y_dim
                @test Av[i] == A[QXZoo.RQC.grid_idx(i, x_dim)...]
            end
        end
    end

    @testset "RQC struct" begin
        xdim=4
        ydim=5
        rqc_ds = QXZoo.RQC.RQC_DS(xdim, ydim)

        @test rqc_ds.circ.num_qubits == xdim*ydim

        @test rqc_ds.single_qubit_gates[1] == GateSymbol(:z, 2, false)
        @test rqc_ds.single_qubit_gates[2] == GateSymbol(:x, 1, false)
        @test rqc_ds.single_qubit_gates[3] == GateSymbol(:y, 1, false)

        for i in 1:xdim*ydim
            @test rqc_ds.next_gate[i] == -1
        end
    end

    @testset "RQC random gate" begin
        # Need a better test-cases here; difficult due to stochastic nature of RQC
        xdim=4
        ydim=5
        depth = 6

        rqc_ds = QXZoo.RQC.RQC_DS(xdim, ydim)
        rqc = QXZoo.RQC.create_RQC(xdim, ydim, depth)

        @test rqc.gate_set == Set((GateSymbol(:z, 2, false), GateSymbol(:x, 1, false), GateSymbol(:y, 1, false), GateSymbol(:c_z, 0, false), GateSymbol(:h, 0, false)))

    end
end


text:
"""
    SeisBinHeaders(in,out; <keyword arguments>)

Sequentially bin seismic headers using the available grid information.

Keyword arguments should be consistent with SeisGeometry keyword arguments.

# Arguments
* `in`: filename of input, irregularly sampled data
* `out`: filename of output, regularly sampled data

# Keyword arguments
* `style="sxsygxgy"`: bin style. Options: "mxmyhxhy","mxmyhaz","sxsyhxhy","gxgyhxhy","sxsyhaz","gxgyhaz"
* `ang=90`: inline direction measured in degrees CC from East
* `gamma=1`: vp/vs ratio for PS Asymptotic Conversion Point gathers (use gamma=1 for PP data)
* `osx=0`,`osy=0`,`ogx=0`,`ogy=0` : origin for source and receiver coordinate system
* `omx=0`,`omy=0`,`ohx=0`,`ohy=0`: origin for midpoint and offset coordinate system
* `oaz=0`,`oh=0` : origin for azimuth and offset coordinate system
* `dsx=1`,`dsy=1`,`dgx=1`,`dgy=1`: source and receiver step-size
* `dmx=1`,`dmy=1`,`dhx=1`,`dhy=1`: midpoint and offset step-size
* `dh=1`,`daz=1`: offset and azimuth step-size
* `min_isx=0`,`max_isx=0`,`min_isy=0`,`max_isy=0`: grid extreme values for sources
* `min_igx=0`,`max_igx=0`,`min_igy=0`,`max_igy=0`: grid extreme values for receivers
* `min_imx=0`,`max_imx=0`,`min_imy=0`,`max_imy=0`: grid extreme values for midpoints
* `min_ihx=0`,`max_ihx=0`,`min_ihy=0`,`max_ihy=0`: grid extreme values for offsets
* `min_ih=0`,`max_ih=0`,`min_iaz=0`,`max_iaz=0`: grid extreme values for azimuth and offset
* `ntrace=10000`: maximum number of traces processed at a time

# Output
In file `out`, binned headers are created.

*Credits: Aaron Stanton,2017*

"""
function SeisBinHeaders(in,out;style="sxsygxgy",ang=90,gamma=1,osx=0,osy=0,ogx=0,ogy=0,omx=0,omy=0,ohx=0,ohy=0,oh=0,oaz=0,dsx=1,dsy=1,dgx=1,dgy=1,dmx=1,dmy=1,dhx=1,dhy=1,dh=1,daz=1,min_isx=0,max_isx=0,min_isy=0,max_isy=0,min_igx=0,max_igx=0,min_igy=0,max_igy=0,min_imx=0,max_imx=0,min_imy=0,max_imy=0,min_ihx=0,max_ihx=0,min_ihy=0,max_ihy=0,min_ih=0,max_ih=0,min_iaz=0,max_iaz=0,ntrace=10000)

r2d = 180/pi;
d2r = pi/180;
gammainv = 1/gamma;
if (ang > 90)
	ang2=-d2r*(ang-90)
else
	ang2=d2r*(90-ang)
end

naz=convert(Int32,360/daz)

if (style=="sxsygxgy")
	nsx = max_isx - min_isx + 1
	nsy = max_isy - min_isy + 1
	ngx = max_igx - min_igx + 1
	ngy = max_igy - min_igy + 1
	nx1=nsx;nx2=nsy;nx3=ngx;nx4=ngy;
	ox1=osx;ox2=osy;ox3=ogx;ox4=ogy;
	dx1=dsx;dx2=dsy;dx3=dgx;dx4=dgy;
	label2="sx";label3="sy";label4="gx";label5="gy";
	unit2="m";unit3="m";unit4="m";unit5="m";
elseif (style=="mxmyhxhy")
	nmx = max_imx - min_imx + 1
	nmy = max_imy - min_imy + 1
	nhx = max_ihx - min_ihx + 1
	nhy = max_ihy - min_ihy + 1
	nx1=nmx;nx2=nmy;nx3=nhx;nx4=nhy;
	ox1=omx;ox2=omy;ox3=ohx;ox4=ohy;
	dx1=dmx;dx2=dmy;dx3=dhx;dx4=dhy;
	label2="mx";label3="my";label4="hx";label5="hy";
	unit2="m";unit3="m";unit4="m";unit5="m";
elseif (style=="mxmyhaz")
	nmx = max_imx - min_imx + 1
	nmy = max_imy - min_imy + 1
	nh  = max_ih  - min_ih  + 1
	naz = max_iaz - min_iaz + 1
	nx1=nmx;nx2=nmy;nx3=nh;nx4=naz;
	ox1=omx;ox2=omy;ox3=oh;ox4=oaz;
	dx1=dmx;dx2=dmy;dx3=dh;dx4=daz;
	label2="mx";label3="my";label4="h";label5="az";
	unit2="m";unit3="m";unit4="m";unit5="Degrees";
elseif (style=="sxsyhxhy")
	nsx = max_isx - min_isx + 1
	nsy = max_isy - min_isy + 1
	nhx = max_ihx - min_ihx + 1
	nhy = max_ihy - min_ihy + 1
	nx1=nsx;nx2=nsy;nx3=nhx;nx4=nhy;
	ox1=osx;ox2=osy;ox3=ohx;ox4=ohy;
	dx1=dsx;dx2=dsy;dx3=dhx;dx4=dhy;
	label2="sx";label3="sy";label4="hx";label5="hy";
	unit2="m";unit3="m";unit4="m";unit5="m";
elseif (style=="gxgyhxhy")
	ngx = max_igx - min_igx + 1
	ngy = max_igy - min_igy + 1
	nhx = max_ihx - min_ihx + 1
	nhy = max_ihy - min_ihy + 1
	nx1=ngx;nx2=ngy;nx3=nhx;nx4=nhy;
	ox1=ogx;ox2=ogy;ox3=ohx;ox4=ohy;
	dx1=dgx;dx2=dgy;dx3=dhx;dx4=dhy;
	label2="gx";label3="gy";label4="hx";label5="hy";
	unit2="m";unit3="m";unit4="m";unit5="m";

elseif (style=="sxsyhaz")
	nsx = max_isx - min_isx + 1
	nsy = max_isy - min_isy + 1
	nh  = max_ih  - min_ih  + 1
	naz = max_iaz - min_iaz + 1
	nx1=nsx;nx2=nsy;nx3=nh;nx4=naz;
	ox1=osx;ox2=osy;ox3=oh;ox4=oaz;
	dx1=dsx;dx2=dsy;dx3=dh;dx4=daz;
	label2="sx";label3="sy";label4="hx";label5="az";
	unit2="m";unit3="m";unit4="m";unit5="Degrees";

elseif (style=="gxgyhaz")
 	ngx = max_igx - min_igx + 1
	ngy = max_igy - min_igy + 1
	nh  = max_ih  - min_ih  + 1
	naz = max_iaz - min_iaz + 1
	nx1=ngx;nx2=ngy;nx3=nh;nx4=naz;
	ox1=ogx;ox2=ogy;ox3=oh;ox4=oaz;
	dx1=dgx;dx2=dgy;dx3=dh;dx4=daz;
	label2="gx";label3="gy";label4="h";label5="az";
	unit2="m";unit3="m";unit4="m";unit5="Degrees";
else
	error("style not recognized.")
end

nx_out = nx1*nx2*nx3*nx4
stream_in = open(ParseHeaderName(in))

nx_in = round(Int,filesize(stream_in)/(4*length(fieldnames(Header))))

seek(stream_in, header_count["n1"])
nt = read(stream_in,Int32)
println("nt= ",nt)
println("The final binned cube will have an approximate size of ", round(nx1*nx2*nx3*nx4*nt*4*1e-9,digits=3)," Gb")
seek(stream_in, header_count["o1"])
o1 = read(stream_in,Float32)
seek(stream_in, header_count["d1"])
dt = read(stream_in,Float32)
close(stream_in)

extent = Extent(convert(Int32,nt),convert(Int32,nx1),convert(Int32,nx2),convert(Int32,nx3),convert(Int32,nx4),
	   convert(Float32,o1),convert(Float32,ox1),convert(Float32,ox2),convert(Float32,ox3),convert(Float32,ox4),
	   convert(Float32,dt),convert(Float32,dx1),convert(Float32,dx2),convert(Float32,dx3),convert(Float32,dx4),
	   "Time",label5,label4,label3,label2,
	   "s",unit5,unit4,unit3,unit2,
	   "")

DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
filename_d = join([DATAPATH out "@data@"])
filename_h = join([DATAPATH out "@headers@"])

#writes extent file
WriteTextHeader(out,extent,"native_float", 4,filename_d,filename_h)

h = Array{Header}(undef,1)
h[1] = InitSeisHeader()
stream_out = open(filename_h,"a+")

j = 1

if (style=="sxsygxgy")

   for ix4 = 1 : nx4
       for ix3 = 1 : nx3
       	  for ix2 = 1 : nx2
		for ix1 = 1 : nx1
	       	  	h[1].tracenum = convert(typeof(h[1].tracenum),j)
			h[1].o1 = convert(typeof(h[1].o1),o1)
			h[1].n1 = convert(typeof(h[1].n1),nt)
			h[1].d1 = convert(typeof(h[1].d1),dt)
			h[1].isx = convert(typeof(h[1].isx),ix1 - 1 + min_isx)
			h[1].isy = convert(typeof(h[1].isy),ix2 - 1 + min_isy)
			h[1].igx = convert(typeof(h[1].igx),ix3 - 1 + min_igx)
			h[1].igy = convert(typeof(h[1].igy),ix4 - 1 + min_igy)
			sx_rot = convert(Float32,(ix1 - 1 + min_isx)*dsx + osx);
			sy_rot = convert(Float32,(ix2 - 1 + min_isy)*dsy + osy);
			gx_rot = convert(Float32,(ix3 - 1 + min_igx)*dgx + ogx);
			gy_rot = convert(Float32,(ix4 - 1 + min_igy)*dgy + ogy);
			h[1].sx =  (sx_rot-osx)*cos(ang2) + (sy_rot-osy)*sin(ang2) + osx;
			h[1].sy = -(sx_rot-osx)*sin(ang2) + (sy_rot-osy)*cos(ang2) + osy;
			h[1].gx =  (gx_rot-ogx)*cos(ang2) + (gy_rot-ogy)*sin(ang2) + ogx;
			h[1].gy = -(gx_rot-ogx)*sin(ang2) + (gy_rot-ogy)*cos(ang2) + ogy;
			h[1].hx = h[1].gx - h[1].sx
			h[1].hy = h[1].gy - h[1].sy
			h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
			h[1].az = r2d*atan((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
			if (h[1].az < 0)
	      			h[1].az += 360.0
			end
			h[1].mx = h[1].sx + h[1].hx/(1 + gammainv);
			h[1].my = h[1].sy + h[1].hy/(1 + gammainv);
			mx_rot = (h[1].mx-omx)*cos(ang2) - (h[1].my-omy)*sin(ang2) + omx;
			my_rot = (h[1].mx-omx)*sin(ang2) + (h[1].my-omy)*cos(ang2) + omy;
			hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
			hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
			h[1].imx = convert(Int32,round((mx_rot-omx)/dmx))
			h[1].imy = convert(Int32,round((my_rot-omy)/dmy))
			h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
			h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
			h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
			h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))<naz ? convert(Int32,round((h[1].az-oaz)/daz)) : 0
			h[1].selev = convert(typeof(h[1].selev),0)
			h[1].gelev = convert(typeof(h[1].gelev),0)
			h[1].trid = convert(typeof(h[1].trid),0)
			PutHeader(stream_out,h[1],j)
			j += 1
		end
     	  end
      end
   end

   h_in=SeisReadHeaders(in)
   hout = SeisReadHeaders(out)

   for k = 1 : nx_in

	itrace = (h_in[k].igy - min_igy)*nx3*nx2*nx1 + (h_in[k].igx - min_igx)*nx2*nx1 + (h_in[k].isy - min_isy)*nx1 + h_in[k].isx - min_isx + 1

	if (itrace > 0 && itrace <= nx_out)
	    	if (h_in[k].isx >= min_isx && h_in[k].isx <= max_isx && h_in[k].isy >= min_isy && h_in[k].isy <= max_isy && h_in[k].igx >= min_igx && h_in[k].igx <= max_igx && h_in[k].igy >= min_igy && h_in[k].igy <= max_igy)

			hout=GrabHeader(stream_out,itrace)

		     	hout.sx=h_in[k].sx
			hout.sy=h_in[k].sy
			hout.gx=h_in[k].gx
			hout.gy=h_in[k].gy
			hout.mx=h_in[k].mx
			hout.my=h_in[k].my
			hout.hx=h_in[k].hx
			hout.hy=h_in[k].hy
			hout.h=h_in[k].h
			hout.az=h_in[k].az
			hout.ang=h_in[k].ang

			PutHeader(stream_out,hout,itrace)
		end
	end
    end

 elseif (style=="mxmyhxhy")

    for ix4 = 1 : nx4
       for ix3 = 1 : nx3
       	  for ix2 = 1 : nx2
		for ix1 = 1 : nx1

			h[1].tracenum = convert(typeof(h[1].tracenum),j)
			h[1].o1 = convert(typeof(h[1].o1),o1)
			h[1].n1 = convert(typeof(h[1].n1),nt)
			h[1].d1 = convert(typeof(h[1].d1),dt)
			h[1].imx = convert(typeof(h[1].imx),ix1 - 1 + min_imx)
			h[1].imy = convert(typeof(h[1].imy),ix2 - 1 + min_imy)
			h[1].ihx = convert(typeof(h[1].ihx),ix3 - 1 + min_ihx)
			h[1].ihy = convert(typeof(h[1].ihy),ix4 - 1 + min_ihy)
			mx_rot = convert(Float32,(ix1 - 1 + min_imx)*dmx + omx);
			my_rot = convert(Float32,(ix2 - 1 + min_imy)*dmy + omy);
			hx_rot = convert(Float32,(ix3 - 1 + min_ihx)*dhx + ohx);
			hy_rot = convert(Float32,(ix4 - 1 + min_ihy)*dhy + ohy);
			h[1].mx =  (mx_rot-omx)*cos(ang2) + (my_rot-omy)*sin(ang2) + omx;
			h[1].my = -(mx_rot-omx)*sin(ang2) + (my_rot-omy)*cos(ang2) + omy;
			h[1].hx =  (hx_rot-ohx)*cos(ang2) + (hy_rot-ohy)*sin(ang2) + ohx;
			h[1].hy = -(hx_rot-ohx)*sin(ang2) + (hy_rot-ohy)*cos(ang2) + ohy;
			h[1].sx = h[1].mx - h[1].hx/(1 + gammainv);
			h[1].sy = h[1].my - h[1].hy/(1 + gammainv);
			h[1].gx = h[1].mx + h[1].hx*(1-(1/(1 + gammainv)));
			h[1].gy = h[1].my + h[1].hy*(1-(1/(1 + gammainv)));
			sx_rot = (h[1].sx-osx)*cos(ang2) - (h[1].sy-osy)*sin(ang2) + osx;
			sy_rot = (h[1].sx-osx)*sin(ang2) + (h[1].sy-osy)*cos(ang2) + osy;
			gx_rot = (h[1].gx-ogx)*cos(ang2) - (h[1].gy-ogy)*sin(ang2) + ogx;
			gy_rot = (h[1].gx-ogx)*sin(ang2) + (h[1].gy-ogy)*cos(ang2) + ogy;
			h[1].isx = convert(Int32,round((sx_rot-osx)/dsx))
			h[1].isy = convert(Int32,round((sy_rot-osy)/dsy))
			h[1].igx = convert(Int32,round((gx_rot-ogx)/dgx))
			h[1].igy = convert(Int32,round((gy_rot-ogy)/dgy))
			h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
			h[1].az = r2d*atan((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
			if (h[1].az < 0)
				h[1].az += 360.0
			end
			h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
			h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))<naz ? convert(Int32,round((h[1].az-oaz)/daz)) : 0
			h[1].selev = convert(typeof(h[1].selev),0)
			h[1].gelev = convert(typeof(h[1].gelev),0)
			h[1].trid = convert(typeof(h[1].trid),0)
			PutHeader(stream_out,h[1],j)

			j += 1
		 end
     	 end
     end
  end

  h_in=SeisReadHeaders(in)
  hout = SeisReadHeaders(out)

  for k = 1 : nx_in

	itrace = (h_in[k].ihy - min_ihy)*nx3*nx2*nx1 + (h_in[k].ihx - min_ihx)*nx2*nx1 + (h_in[k].imy - min_imy)*nx1 + h_in[k].imx - min_imx + 1
	if (itrace > 0 && itrace <= nx_out)
	 if (h_in[k].imx >= min_imx && h_in[k].imx <= max_imx && h_in[k].imy >= min_imy && h_in[k].imy <= max_imy && h_in[k].ihx >= min_ihx && h_in[k].ihx <= max_ihx && h_in[k].ihy >= min_ihy && h_in[k].ihy <= max_ihy)

		 hout=GrabHeader(stream_out,itrace)
		 hout.sx=h_in[k].sx
		 hout.sy=h_in[k].sy
		 hout.gx=h_in[k].gx
		 hout.gy=h_in[k].gy
		 hout.mx=h_in[k].mx
		 hout.my=h_in[k].my
		 hout.hx=h_in[k].hx
		 hout.hy=h_in[k].hy
		 hout.h=h_in[k].h
		 hout.az=h_in[k].az
		 hout.ang=h_in[k].ang
		 PutHeader(stream_out,hout,itrace)

	 end
       end
 end

elseif (style=="mxmyhaz")
	for ix4 = 1 : nx4
         for ix3 = 1 : nx3
	  for ix2 = 1 : nx2
	   for ix1 = 1 : nx1

       	 	h[1].tracenum = convert(typeof(h[1].tracenum),j)
 		h[1].o1 = convert(typeof(h[1].o1),o1)
		h[1].n1 = convert(typeof(h[1].n1),nt)
		h[1].d1 = convert(typeof(h[1].d1),dt)
		h[1].imx = convert(typeof(h[1].imx),ix1 - 1 + min_imx)
		h[1].imy = convert(typeof(h[1].imy),ix2 - 1 + min_imy)
		h[1].ih = convert(typeof(h[1].ih),  ix3 - 1 + min_ih)
		h[1].iaz = convert(typeof(h[1].iaz),ix4 - 1 + min_iaz)
		mx_rot = convert(Float32,(ix1 - 1 + min_imx)*dmx + omx);
		my_rot = convert(Float32,(ix2 - 1 + min_imy)*dmy + omy);
		h[1].mx =  (mx_rot-omx)*cos(ang2) + (my_rot-omy)*sin(ang2) + omx;
		h[1].my = -(mx_rot-omx)*sin(ang2) + (my_rot-omy)*cos(ang2) + omy;
		h[1].h = convert(Float32,(ix3 - 1 + min_ih)*dh + oh);
		h[1].az = r2d*atan((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
		if (h[1].az <= 90)
			h[1].hx = h[1].h*cos(d2r*h[1].az);
			h[1].hy = h[1].h*sin(d2r*h[1].az);
		elseif (h[1].az > 90 && h[1].az <= 180)
			h[1].hx =-h[1].h*cos(pi-(d2r*h[1].az));
			h[1].hy = h[1].h*sin(pi-(d2r*h[1].az));
		elseif (h[1].az > 180 && h[1].az <= 270)
			h[1].hx =-h[1].h*cos((d2r*h[1].az)-pi);
			h[1].hy =-h[1].h*sin((d2r*h[1].az)-pi);
		else
			h[1].hx = h[1].h*cos(2*pi-(d2r*h[1].az));
			h[1].hy =-h[1].h*sin(2*pi-(d2r*h[1].az));
		end
		h[1].sx = h[1].mx - h[1].hx/(1 + gammainv);
		h[1].sy = h[1].my - h[1].hy/(1 + gammainv);
		h[1].gx = h[1].mx + h[1].hx*(1-(1/(1 + gammainv)));
		h[1].gy = h[1].my + h[1].hy*(1-(1/(1 + gammainv)));
		sx_rot = (h[1].sx-osx)*cos(ang2) - (h[1].sy-osy)*sin(ang2) + osx;
		sy_rot = (h[1].sx-osx)*sin(ang2) + (h[1].sy-osy)*cos(ang2) + osy;
		gx_rot = (h[1].gx-ogx)*cos(ang2) - (h[1].gy-ogy)*sin(ang2) + ogx;
		gy_rot = (h[1].gx-ogx)*sin(ang2) + (h[1].gy-ogy)*cos(ang2) + ogy;
		h[1].isx = convert(Int32,round((sx_rot-osx)/dsx))
		h[1].isy = convert(Int32,round((sy_rot-osy)/dsy))
		h[1].igx = convert(Int32,round((gx_rot-ogx)/dgx))
		h[1].igy = convert(Int32,round((gy_rot-ogy)/dgy))
		hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
		hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
		h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
		h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
		h[1].selev = convert(typeof(h[1].selev),0)
		h[1].gelev = convert(typeof(h[1].gelev),0)
		h[1].trid = convert(typeof(h[1].trid),0)
		PutHeader(stream_out,h[1],j)

		j += 1
	end
     end
    end
   end

 h_in=SeisReadHeaders(in)
 hout = SeisReadHeaders(out)

 for k = 1 : nx_in
	itrace = (h_in[k].iaz - min_iaz)*nx3*nx2*nx1 + (h_in[k].ih - min_ih)*nx2*nx1 + (h_in[k].imy - min_imy)*nx1 + h_in[k].imx - min_imx + 1
	 if (itrace > 0 && itrace <= nx_out)
		 if (h_in[k].imx >= min_imx && h_in[k].imx <= max_imx && h_in[k].imy >= min_imy && h_in[k].imy <= max_imy && h_in[k].ih >= min_ih && h_in[k].ih <= max_ih && h_in[k].iaz >= min_iaz && h_in[k].iaz <= max_iaz)

			 hout[itrace].sx=h_in[k].sx
			 hout[itrace].sy=h_in[k].sy
			 hout[itrace].gx=h_in[k].gx
			 hout[itrace].gy=h_in[k].gy
			 hout[itrace].mx=h_in[k].mx
			 hout[itrace].my=h_in[k].my
			 hout[itrace].hx=h_in[k].hx
			 hout[itrace].hy=h_in[k].hy
			 hout[itrace].h=h_in[k].h
			 hout[itrace].az=h_in[k].az
			 hout[itrace].ang=h_in[k].ang
			 PutHeader(stream_out,hout[itrace],itrace)

		 end
	 end
 end

elseif (style=="sxsyhxhy")
 for ix4 = 1 : nx4
  for ix3 = 1 : nx3
   for ix2 = 1 : nx2
    for ix1 = 1 : nx1

	h[1].tracenum = convert(typeof(h[1].tracenum),j)
	h[1].o1 = convert(typeof(h[1].o1),o1)
	h[1].n1 = convert(typeof(h[1].n1),nt)
	h[1].d1 = convert(typeof(h[1].d1),dt)
	h[1].isx = convert(typeof(h[1].isx),ix1 - 1 + min_isx)
	h[1].isy = convert(typeof(h[1].isy),ix2 - 1 + min_isy)
	h[1].ihx = convert(typeof(h[1].ihx),ix3 - 1 + min_ihx)
	h[1].ihy = convert(typeof(h[1].ihy),ix4 - 1 + min_ihy)
	sx_rot = convert(Float32,(ix1 - 1 + min_isx)*dsx + osx);
	sy_rot = convert(Float32,(ix2 - 1 + min_isy)*dsy + osy);
	hx_rot = convert(Float32,(ix3 - 1 + min_ihx)*dhx + ohx);
	hy_rot = convert(Float32,(ix4 - 1 + min_ihy)*dhy + ohy);
	h[1].sx =  (sx_rot-osx)*cos(ang2) + (sy_rot-osy)*sin(ang2) + osx;
	h[1].sy = -(sx_rot-osx)*sin(ang2) + (sy_rot-osy)*cos(ang2) + osy;
	h[1].hx =  (hx_rot-ohx)*cos(ang2) + (hy_rot-ohy)*sin(ang2) + ohx;
	h[1].hy = -(hx_rot-ohx)*sin(ang2) + (hy_rot-ohy)*cos(ang2) + ohy;
	h[1].gx = h[1].sx + h[1].hx;
	h[1].gy = h[1].sy + h[1].hy;
	h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
	h[1].az = r2d*atan((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
	if (h[1].az < 0)
		h[1].az += 360.0
	end
	h[1].mx = h[1].sx + h[1].hx/(1 + gammainv);
	h[1].my = h[1].sy + h[1].hy/(1 + gammainv);
	mx_rot = (h[1].mx-omx)*cos(ang2) - (h[1].my-omy)*sin(ang2) + omx;
	my_rot = (h[1].mx-omx)*sin(ang2) + (h[1].my-omy)*cos(ang2) + omy;
	hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
	hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
	h[1].imx = convert(Int32,round((mx_rot-omx)/dmx))
	h[1].imy = convert(Int32,round((my_rot-omy)/dmy))
	h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
	h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
	h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
	h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))<naz ? convert(Int32,round((h[1].az-oaz)/daz)) : 0
	h[1].selev = convert(typeof(h[1].selev),0)
	h[1].gelev = convert(typeof(h[1].gelev),0)
	h[1].trid = convert(typeof(h[1].trid),0)
	PutHeader(stream_out,h[1],j)

	j += 1
    end
   end
  end
 end

 h_in=SeisReadHeaders(in)
 hout = SeisReadHeaders(out)

 for k = 1 : nx_in
  itrace = (h_in[k].ihy - min_ihy)*nx3*nx2*nx1 + (h_in[k].ihx - min_ihx)*nx2*nx1 + (h_in[k].isy - min_isy)*nx1 + h_in[k].isx - min_isx + 1
	if (itrace > 0 && itrace <= nx_out)

	 if (h_in[k].isx >= min_isx && h_in[k].isx <= max_isx && h_in[k].isy >= min_isy && h_in[k].isy <= max_isy && h_in[k].ihx >= min_ihx && h_in[k].ihx <= max_ihx && h_in[k].ihy >= min_ihy && h_in[k].ihy <= max_ihy)

	 	hout[itrace].sx=h_in[k].sx
		hout[itrace].sy=h_in[k].sy
		hout[itrace].gx=h_in[k].gx
		hout[itrace].gy=h_in[k].gy
		hout[itrace].mx=h_in[k].mx
		hout[itrace].my=h_in[k].my
		hout[itrace].hx=h_in[k].hx
		hout[itrace].hy=h_in[k].hy
		hout[itrace].h=h_in[k].h
		hout[itrace].az=h_in[k].az
		hout[itrace].ang=h_in[k].ang
                PutHeader(stream_out,hout[itrace],itrace)

         end
      end
  end

elseif (style=="gxgyhxhy")
 for ix4 = 1 : nx4
  for ix3 = 1 : nx3
   for ix2 = 1 : nx2
    for ix1 = 1 : nx1

	h[1].tracenum = convert(typeof(h[1].tracenum),j)
	h[1].o1 = convert(typeof(h[1].o1),o1)
	h[1].n1 = convert(typeof(h[1].n1),nt)
	h[1].d1 = convert(typeof(h[1].d1),dt)
	h[1].igx = convert(typeof(h[1].igx),ix1 - 1 + min_igx)
	h[1].igy = convert(typeof(h[1].igy),ix2 - 1 + min_igy)
	h[1].ihx = convert(typeof(h[1].ihx),ix3 - 1 + min_ihx)
	h[1].ihy = convert(typeof(h[1].ihy),ix4 - 1 + min_ihy)
	gx_rot = convert(Float32,(ix1 - 1 + min_igx)*dgx + ogx);
	gy_rot = convert(Float32,(ix2 - 1 + min_igy)*dgy + ogy);
	hx_rot = convert(Float32,(ix3 - 1 + min_ihx)*dhx + ohx);
	hy_rot = convert(Float32,(ix4 - 1 + min_ihy)*dhy + ohy);
	h[1].gx =  (gx_rot-ogx)*cos(ang2) + (gy_rot-ogy)*sin(ang2) + ogx;
	h[1].gy = -(gx_rot-ogx)*sin(ang2) + (gy_rot-ogy)*cos(ang2) + ogy;
	h[1].hx =  (hx_rot-ohx)*cos(ang2) + (hy_rot-ohy)*sin(ang2) + ohx;
	h[1].hy = -(hx_rot-ohx)*sin(ang2) + (hy_rot-ohy)*cos(ang2) + ohy;
	h[1].sx = h[1].gx - h[1].hx;
	h[1].sy = h[1].gy - h[1].hy;
	h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
	h[1].az = r2d*atan((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
	if (h[1].az < 0)
		h[1].az += 360.0
	end
	h[1].mx = h[1].sx + h[1].hx/(1 + gammainv);
	h[1].my = h[1].sy + h[1].hy/(1 + gammainv);
	mx_rot = (h[1].mx-omx)*cos(ang2) - (h[1].my-omy)*sin(ang2) + omx;
	my_rot = (h[1].mx-omx)*sin(ang2) + (h[1].my-omy)*cos(ang2) + omy;
	hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
	hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
	h[1].imx = convert(Int32,round((mx_rot-omx)/dmx))
	h[1].imy = convert(Int32,round((my_rot-omy)/dmy))
	h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
	h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
	h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
	h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))<naz ? convert(Int32,round((h[1].az-oaz)/daz)) : 0
	h[1].selev = convert(typeof(h[1].selev),0)
	h[1].gelev = convert(typeof(h[1].gelev),0)
	h[1].trid = convert(typeof(h[1].trid),0)
	PutHeader(stream_out,h[1],j)

	j += 1
     end
    end
   end
  end

 h_in=SeisReadHeaders(in)
 hout = SeisReadHeaders(out)

 for k = 1 : nx_in
	itrace = (h_in[k].ihy - min_ihy)*nx3*nx2*nx1 + (h_in[k].ihx - min_ihx)*nx2*nx1 + (h_in[k].igy - min_igy)*nx1 + h_in[k].igx - min_igx + 1
	if (itrace > 0 && itrace <= nx_out)

	 if (h_in[k].igx >= min_igx && h_in[k].igx <= max_igx && h_in[k].igy >= min_igy && h_in[k].igy <= max_igy && h_in[k].ihx >= min_ihx && h_in[k].ihx <= max_ihx && h_in[k].ihy >= min_ihy && h_in[k].ihy <= max_ihy)

		 hout[itrace].sx=h_in[k].sx
		 hout[itrace].sy=h_in[k].sy
		 hout[itrace].gx=h_in[k].gx
		 hout[itrace].gy=h_in[k].gy
		 hout[itrace].mx=h_in[k].mx
		 hout[itrace].my=h_in[k].my
		 hout[itrace].hx=h_in[k].hx
		 hout[itrace].hy=h_in[k].hy
		 hout[itrace].h=h_in[k].h
		 hout[itrace].az=h_in[k].az
		 hout[itrace].ang=h_in[k].ang
		 PutHeader(stream_out,hout[itrace],itrace)

	 end
       end
     end

elseif (style=="sxsyhaz")
	error("sxsyhaz not developed yet.")

elseif (style=="gxgyhaz")
	error("gxgyhaz not developed yet.")

end

close(stream_out)

end


text:
"""
    _calculate_interval_inner_counts(order::Dict{Int,String},
                                          intervals::OrderedDict{String,<:Dates.TimePeriod},
                                          step_resolution::Dates.TimePeriod)

Calculates how many times a problem is executed for every interval of the previous problem
"""
function _calculate_interval_inner_counts(
    intervals::OrderedDict{Symbol, Tuple{Dates.TimePeriod, <:FeedForwardChronology}},
)
    order = collect(keys(intervals))
    reverse_order = length(intervals):-1:1
    interval_run_counts = Dict{Int, Int}(1 => 1)
    for k in reverse_order[1:(end - 1)]
        problem_name = order[k]
        previous_problem_name = order[k - 1]
        problem_interval = intervals[problem_name][1]
        previous_problem_interval = intervals[previous_problem_name][1]
        if Dates.Millisecond(previous_problem_interval % problem_interval) !=
           Dates.Millisecond(0)
            throw(
                IS.ConflictingInputsError(
                    "The interval configuration provided results in a fractional number of executions of problem $problem_name",
                ),
            )
        end
        interval_run_counts[k] = previous_problem_interval / problem_interval
        @debug "problem $k is executed $(interval_run_counts[k]) time within each interval of problem $(k-1)"
    end
    return interval_run_counts
end

""" Function calculates the total number of problem executions in the simulation and allocates the appropiate vector"""
function _allocate_execution_order(interval_run_counts::Dict{Int, Int})
    total_size_of_vector = 0
    for (k, counts) in interval_run_counts
        mult = 1
        for i in 1:k
            mult *= interval_run_counts[i]
        end
        total_size_of_vector += mult
    end
    return -1 * ones(Int, total_size_of_vector)
end

function _fill_execution_order(
    execution_order::Vector{Int},
    interval_run_counts::Dict{Int, Int},
)
    function _fill_problem(index::Int, problem::Int)
        if problem < last_problem
            next_problem = problem + 1
            for i in 1:interval_run_counts[next_problem]
                index = _fill_problem(index, next_problem)
            end
        end
        execution_order[index] = problem
        index -= 1
    end

    index = length(execution_order)
    problems = sort!(collect(keys(interval_run_counts)))
    last_problem = problems[end]
    _fill_problem(index, problems[1])
    return
end

function _get_execution_order_vector(
    intervals::OrderedDict{Symbol, Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}},
)
    length(intervals) == 1 && return [1]
    interval_run_counts = _calculate_interval_inner_counts(intervals)
    execution_order_vector = _allocate_execution_order(interval_run_counts)
    _fill_execution_order(execution_order_vector, interval_run_counts)
    @assert isempty(findall(x -> x == -1, execution_order_vector))
    return execution_order_vector
end

function _check_feedforward(
    feedforward::Dict{Tuple{String, Symbol, Symbol}, <:AbstractAffectFeedForward},
    feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology},
)
    isempty(feedforward) && return
    for problem_key in keys(feedforward)
        @debug problem_key
        if !mapreduce(x -> x.second == problem_key[1], |, keys(feedforward_chronologies))
            throw(
                ArgumentError(
                    "No valid Chronology has been defined for the feedforward added to $(problem_key[1])",
                ),
            )
        end
    end
    return
end

function _check_chronology_consistency(
    problems::SimulationProblems,
    feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology},
    ini_cond_chronology::InitialConditionChronology,
)
    if isempty(feedforward_chronologies)
        @warn(
            "No Feedforward Chronologies have been defined. This configuration assummes that there is no information passing between problems"
        )
    end
    if length(problems) == 1
        if isa(ini_cond_chronology, InterProblemChronology)
            @warn(
                "Single problem detected, the default Initial Condition Chronology is IntraProblemChronology(), other values will be ignored."
            )
        end
    end
    # TODO: Add more consistency checks
    return
end

function _check_cache_definition(cache::Dict{<:Tuple, <:AbstractCache})
    for (problem_names, c) in cache
        if typeof(c) == TimeStatusChange && length(problem_names) > 1
            error(
                "TimeStatusChange cache currently only supports single problem. Please consider changing your cache definitions",
            )
        end
    end
    return
end

function _get_num_executions_by_problem(problems, execution_order)
    problem_names = get_problem_names(problems)
    executions_by_problem = Dict(x => 0 for x in problem_names)
    for problem_number in execution_order
        executions_by_problem[problem_names[problem_number]] += 1
    end
    return executions_by_problem
end

@doc raw"""
    SimulationSequence(horizons::Dict{String, Int}
                        step_resolution::Dates.TimePeriod
                        intervals::Dict{String, <:Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}}
                        order::Dict{Int, String}
                        feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology}
                        feedforward::Dict{Tuple{String, Symbol, Symbol}, <:AbstractAffectFeedForward}
                        ini_cond_chronology::Dict{String, <:FeedForwardChronology}
                        cache::Dict{String, AbstractCache}
                        )
"""
mutable struct SimulationSequence
    horizons::OrderedDict{Symbol, Int}
    # JDNOTE: This field might be able to go away.
    step_resolution::Dates.TimePeriod
    # The string here is the name of the problem
    intervals::OrderedDict{Symbol, Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}}
    feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology}
    feedforward::Dict{Tuple{String, Symbol, Symbol}, <:AbstractAffectFeedForward}
    ini_cond_chronology::InitialConditionChronology
    cache::Dict{Tuple, AbstractCache}
    execution_order::Vector{Int}
    executions_by_problem::Dict{Symbol, Int}
    current_execution_index::Int64
    uuid::Base.UUID

    function SimulationSequence(;
        problems::SimulationProblems,
        # JDNOTE: We could remove interval here later
        intervals::Dict{String, <:Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}},
        feedforward_chronologies = Dict{Pair{String, String}, FeedForwardChronology}(),
        feedforward = Dict{Tuple{String, Symbol, Symbol}, AbstractAffectFeedForward}(),
        ini_cond_chronology = InterProblemChronology(),
        cache = Dict{Tuple, AbstractCache}(),
    )
        horizons = determine_horizons!(problems)
        _intervals =
            OrderedDict{Symbol, Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}}()
        for k in get_problem_names(problems)
            # JDNOTE: Temporary conversion while we re-define how to do this
            k_ = string(k)
            if !(k_ in keys(intervals))
                throw(IS.ConflictingInputsError("Interval not defined for problem $(k_)"))
            end
            _intervals[k] = (IS.time_period_conversion(intervals[k_][1]), intervals[k_][2])
        end
        step_resolution = determine_step_resolution(_intervals)
        _check_feedforward(feedforward, feedforward_chronologies)
        _check_chronology_consistency(
            problems,
            feedforward_chronologies,
            ini_cond_chronology,
        )
        _check_cache_definition(cache)
        if length(problems) == 1
            ini_cond_chronology = IntraProblemChronology()
        end
        execution_order = _get_execution_order_vector(_intervals)
        executions_by_problem = _get_num_executions_by_problem(problems, execution_order)
        sequence_uuid = IS.make_uuid()
        initialize_simulation_internals!(problems, sequence_uuid)
        new(
            horizons,
            step_resolution,
            _intervals,
            feedforward_chronologies,
            feedforward,
            ini_cond_chronology,
            cache,
            execution_order,
            executions_by_problem,
            0,
            sequence_uuid,
        )
    end
end

get_step_resolution(sequence::SimulationSequence) = sequence.step_resolution

function get_problem_interval_chronology(sequence::SimulationSequence, problem)
    return sequence.intervals[Symbol(problem)][2]
end

function get_interval(sequence::SimulationSequence, problem::Symbol)
    return sequence.intervals[problem][1]
end

function get_interval(sequence::SimulationSequence, problem)
    return sequence.intervals[Symbol(problem)][1]
end

get_execution_order(sequence::SimulationSequence) = sequence.execution_order


text:
"""
    file_paths(spec::SpatQSimSpec)

Returns paths for population state HDF5 array, population total CSV, and catch
CSV. Will create directories as needed. For example, for the 5th realization
using the second preference intensity value, this function will return a vector
with the strings:

    repl_05/prefintensity_02_popstate.h5
    repl_05/prefintensity_02_popstate.csv
    repl_05/prefintensity_02_catch.csv
    repl_05/prefintensity_02_popstate.feather
    repl_05/prefintensity_02_catch.feather
"""
function file_paths(spec::SpatQSimSpec; base_dir = ".")
    rlz = realization(spec)

    base_dir = normpath(base_dir,
                        simstudy_dir(spec),
                         "repl_" * string(rlz, pad = 2))
    if !isdir(base_dir)
        mkpath(base_dir)
    end

    file_base = joinpath(base_dir,
                         simstudy_prefix(spec) *
                         string(sim_value_idx(spec), pad = 2))
    file_base .* ["_popstate.h5",
                  "_popstate.csv",
                  "_catch.csv",
                  "_popstate.feather",
                  "_catch.feather",
                  "_prep.h5"]
end

"""
Return the directory used to store simulations for the simulation study.
"""
simstudy_dir(::T) where T<:SpatQSimSpec = simstudy_dir(T)
simstudy_dir(::Type{<:QDevScalingSpec}) = "qdevscaling"
simstudy_dir(::Type{<:SharedQSpec}) = "sharedq"
simstudy_dir(::Type{<:PrefIntensitySpec}) = "prefintensity"
simstudy_dir(::Type{<:DensityDependentQSpec}) = "densdepq"
simstudy_dir(::Type{<:HabQSpec}) = "habq"

"""
Return the file prefix for the simulation study.
"""
simstudy_prefix(::T) where T<:SpatQSimSpec = simstudy_prefix(T)
simstudy_prefix(::Type{<:QDevScalingSpec}) = "qdevscale_"
simstudy_prefix(::Type{<:SharedQSpec}) = "sharedq_"
simstudy_prefix(::Type{<:PrefIntensitySpec}) = "prefintensity_"
simstudy_prefix(::Type{<:DensityDependentQSpec}) = "densdepq_"
simstudy_prefix(::Type{<:HabQSpec}) = "habq_"

function prep_path(simtype::Type{<:HabQSpec}, repl::Integer; base_dir = ".")
    normpath(base_dir,
             simstudy_dir(simtype),
             "repl_" * string(repl, pad = 2),
             simstudy_prefix(simtype) * "prep.h5")
end

function prep_path(spec::HabQSpec; base_dir = ".")
    prep_path(typeof(spec), realization(spec); base_dir = base_dir)
end