text:
# Elementwise addition of two Datasets, assuming same size
function (+)(a::MultiChannelData, b::MultiChannelData)
    if a.channellayout.nobs_per_channel != b.channellayout.nobs_per_channel
        throw(ErrorException("Data observations per channel not equal"))
    elseif a.nsamples != b.nsamples
        throw(ErrorException("Number of samples not equal"))
    else
        vals = [[c1 .+ c2 for (c1, c2) in zip(samp1, samp2)] for (samp1, samp2) in zip(a.values, b.values)]
        out = MultiChannelData(vals, a.channellayout, a.nsamples)
        return out
    end
end

function eigenvalue_ratio_test(evs :: Array{Float64, 1}; threshold=5)
	# assume evs all real, sorted from largest to smallest, strictly positive
	# Could use random matrix theory to pick thresh.
	evratios = [min(ev1/ev2, typemax(Int64)) for (ev1, ev2) in zip(evs[1:(end-1)], evs[2:end])]
	rmax = maximum(evratios)
	if rmax < threshold
		facnum=0
	else
		facnum = argmax(evratios)
	end
	return facnum
end

function LTOrthog(B, blocksize)
    Blq = lq(B[1:blocksize, 1:blocksize])
    # Extract diagonal elems, sort into descending abs val, put back in diag mat.
    Lsigns = diagm(sign.(sort(diag(Blq.L), lt=(x, y) -> abs(y) < abs(x) )))
    Q = transpose(Blq.Q) * Lsigns  # eqn 11
    return Q
end

function extract_horizontal_blocks(mat, blocksizes; get_view=false)
    rbounds = vcat([0], cumsum(blocksizes))
    if get_view
        blocks = [view(mat, start:stop, :) for (start, stop) in zip(rbounds .+ 1, rbounds[2:end])]
    else
        blocks = [mat[start:stop, :] for (start, stop) in zip(rbounds .+ 1, rbounds[2:end])]
    end
    return blocks
end

function extract_diagonal_blocks(mat, rowblocksizes, colblocksizes; get_view=false)
    rbounds = vcat([0], cumsum(rowblocksizes))
    cbounds = vcat([0], cumsum(colblocksizes))
    if get_view
        blocks =  [view( mat, rstart:rstop, cstart:cstop)
            for (rstart, rstop, cstart, cstop) in zip(rbounds .+ 1, rbounds[2:end], cbounds .+ 1, cbounds[2:end])]
    else
        blocks = [mat[rstart:rstop, cstart:cstop]
            for (rstart, rstop, cstart, cstop) in zip(rbounds .+ 1, rbounds[2:end], cbounds .+ 1, cbounds[2:end])]
    end
    return blocks
end

# function extract_blocks(mat, Ls, pjs; diag=false)
#     rbounds = vcat([0], cumsum(Ls))
#     cbounds = vcat([0], cumsum(pjs))
#     if diag
#         # Pull out diagonal blocks
#         blocks = [mat[rstart:rstop, cstart:cstop]
#                     for (rstart, rstop, cstart, cstop) in zip(rbounds .+ 1, rbounds[2:end], cbounds .+ 1, cbounds[2:end])]
#     else
#         # Split rows of matrix into pieces.
#         blocks = [mat[start:stop, :] for (start, stop) in zip(rbounds .+ 1, rbounds[2:end])]
#     end
#     return blocks
# end

function flip_signs(mat)
    diag_signs = sign.(diag(mat))
    flip= mat * diagm(diag_signs)
    return flip
end

function blockdiag(blks; sparse=false)
    if sparse
        throw("unimplemented")
    else
        ncols = [size(blk)[2] for blk in blks]
        nrows = [size(blk)[1] for blk in blks]
        blkmat = zeros(sum(nrows), sum(ncols))
        colbounds = vcat([0], cumsum(ncols))
        rowbounds = vcat([0], cumsum(nrows))
        for (i, blk) in enumerate(blks)
            cstart, cstop = colbounds[i]+ 1, colbounds[i+1]
            rstart, rstop = rowbounds[i]+ 1, rowbounds[i+1]
            blkmat[rstart:rstop, cstart:cstop] = blk
        end
    end
    return blkmat
end

function stack(data :: MultiChannelData)
    X = hcat([vcat(samp...) for samp in data.values ]...)
    @assert size(X) == (sum(data.channellayout.nobs_per_channel), data.nsamples)
    return X
end

function unstack(X :: Matrix, channellayout::MultiChannelLayout; obs_ids = [], sample_ids = [])
    nsamples = size(X)[2]
    values = Array{Array{Array{Float64, 1}, 1}, 1}(undef, nsamples)
    for samp in 1:nsamples
        v = extract_horizontal_blocks(X[:, samp], channellayout.nobs_per_channel)
        values[samp] = [dropdims(vb; dims=2) for vb in v]
    end
    data = MultiChannelData(values, channellayout, nsamples; obs_ids, sample_ids)
    return data
end

#
# function stack_and_view(data)
#     N = length(data)
#     x1, xview1 = stack_and_view_single(data[1])
#     M = length(xview1)
#     L = sum([length(xc) for xc in x1])
#     X = Matrix(zeros((L, N))) # M x N, flat matrix
#     Xobs = repeat([xview1], N) # Obs-first, N x M x L_m
#     Xchan = [[Xobs[n][m] for n in 1:N] for m in 1:M ] # Chan first, M x N x L_m
#     X[:, 1] = x1
#     for n in 2:length(data)
#         #TODO: This doesn't accomplish the desired view, as view is to data
#         # underlying x, which is copied into matrix X. Need view into X itself,
#         # which requires rewriting this.
#         x, xobs = stack_and_view_single(data[n])
#         X[:, n] = x
#         Xobs[n] = xobs
#         for m in 1:M
#             Xchan[m][n] = xobs[m]
#         end
#     end
#     return X, Xobs, Xchan
# end
#
#
# # Probably need to separaete viewing and stacking code.
# function stack_and_view_single(xs)
#     # Combine ragged 3-dim array into matrix + view with struct.
#     lengs = [length(x) for x in xs]
#     bounds = vcat([0], cumsum(lengs))
#     x = vcat(xs...)
#     xview = [view(x, start:stop) for (start, stop) in zip(bounds .+ 1, bounds[2:end])]
#     return x, xview
# end

function data_fromlong(df; should_interpolate=true, max_missingpcent = 0.3 )
	#df must have 4 columns: "sample_id", "channel", "observation_id", "value"
	samp_ids = sort(unique(df.sample_id))
	channels = unique(df.channel)
	chanouts = Array{Array{Float64, 2}, 1}(undef, length(channels))
	nobs_per_channel = zeros(length(channels))
	k = 1
	local included_obs
	for chandf in groupby(df, :channel)
		obs_ids  = sort(unique(chandf.observation_id))
		included_rows = []
		included_obs = []
		chanout = zeros(length(obs_ids), length(samp_ids))
		for (j ,obs) in enumerate(obs_ids)
			obsdf = filter(:observation_id => id -> id == obs, chandf)
			sort!(obsdf, :sample_id)
			ts = obsdf.sample_id
			vals = obsdf.value
			outvals = zeros(length(samp_ids))
			if length(ts) != length(samp_ids)
				if should_interpolate
					if length(ts) <= (1-max_missingpcent)*length(samp_ids)
						# too many missing, should not be included in model.
						print("$obs has too many missing values, removed.\n")
						# Some misingness pattern problem though. What if different
						# numbers of obs for same obs_id are missing in different channels?
						continue
					else
						int = LinearInterpolation(ts, vals, extrapolation_bc=Line())
						interpolated_values = int(samp_ids)
						# outvals is full-length vector, where the observed values
						# are present, and missing values filled in
						for (i, t) in enumerate(samp_ids)
							outvals[i] = t in ts ? vals[findfirst(ts .== t)] : interpolated_values[i]
						end
					end
				else
					throw(ErrorException("Missing observations and not interpolating!"))
				end
			else
				outvals = vals
			end
			push!(included_obs, obs)
			push!(included_rows, j)
			chanout[j, :] = outvals
		end
		nobs_per_channel[k] = length(included_obs)
		chanouts[k] = chanout[included_rows, :]
		k += 1
	end
	stacked_X = cat(chanouts...; dims=1) :: Array{Float64, 2}
	convert(Array{Float64, 2}, stacked_X)
	channellayout = MultiChannelLayout(nobs_per_channel)
	data = unstack(stacked_X, channellayout; obs_ids = included_obs, sample_ids = samp_ids)
	return data
end

function to_long(params::MCFMParams; obs_ids = Int[])
	nchannels = params.channellayout.nchannels
	h_longs = Array{DataFrame, 1}(undef, nchannels)
	g_longs = Array{DataFrame, 1}(undef, nchannels)
	for (c, Hchan) in enumerate(params.H_by_channel)
		df = DataFrame(Hchan)
		rename!(df, ["common" * string(i) for i in 1:ncol(df)])
		df.obs = length(obs_ids) == nrow(df) ? obs_ids : axes(df, 1)
		df_long = DataFrames.stack(df, Not(:obs))
		df_long[!, :channel] .= c
		rename!(df_long, :variable => "factor")
		h_longs[c] = df_long
	end
	for (c, Gchan) in enumerate(params.G_by_channel)
		df = DataFrame(Gchan)
		rename!(df, ["channel$(c)_unique" * string(i) for i in 1:ncol(df)])
		df.obs = length(obs_ids) == nrow(df) ? obs_ids : axes(df, 1)
		df_long = DataFrames.stack(df, Not(:obs))
		df_long[!, :channel] .= c
		rename!(df_long, :variable => "factor")
		g_longs[c] = df_long
	end

	Hlong = vcat(h_longs...)
	Glong = vcat(g_longs...)
	out = vcat(Hlong, Glong)
	return out
end

function to_long(factors :: MultiChannelFactors; sample_ids = Int[])
	factor_names = ["common" * string(i) for i in  1:factors.factorlayout.nchannelshared]
	for (c, nc) in enumerate(factors.factorlayout.nspecific_per_channel)
		factor_names = vcat(factor_names, ["channel$(c)_unique" * string(i) for i in 1:nc])
	end
	df = DataFrame(transpose(factors.factors))
	rename!(df, factor_names)
	df.sample_id = length(sample_ids) == nrow(df) ? sample_ids : axes(df, 1)
	df_long = DataFrames.stack(df, Not(:sample_id))
	rename!(df_long, :variable => "factor")
	return df_long
end

function to_long(data :: MultiChannelData; obs_ids = Int[], sample_ids = Int[])
	X = MultiChannelFactorAnalysis.stack(data)
	nchannels = data.channellayout.nchannels
	chanblocks = extract_horizontal_blocks(X, data.channellayout.nobs_per_channel)
	Xlongs = Array{DataFrame, 1}(undef, nchannels)
	for (c, blk) in enumerate(chanblocks)
		df = DataFrame(transpose(blk))
		df.sample_id = length(sample_ids) == nrow(df) ? sample_ids : axes(df, 1)
		df_long = DataFrames.stack(df, Not(:sample_id))
		df_long[!, :channel] .= c
		rename!(df_long, :variable => "obs_str")
		# Drop the "x" in front of obs_id (stored as Categorical by default.)
		fixed_obs = [parse(Int, levels(x)[1][2:end]) for x in df_long[:obs_str]]
		df_long[!, :obs] = fixed_obs
		if length(obs_ids) == ncol(df) -1 # added sample_id
			# index into the given ids (assuming in order).
			df_long[!, :obs] = obs_ids[df_long[!, :obs]]
		end
		select!(df_long, Not(:obs_str))
		Xlongs[c] = df_long
	end
	out = vcat(Xlongs...)
	return out
end

function parse_df(df)
    #df must have 3 columns: "sample", "channel", "observation"
    # must be sorted appropriately
    nsamples = length(unique(df.sample))
    nchannel = length(unique(df.channel))

    values = Array{Array{Array{Float64, 1}, 1}, 1}(undef, nsamples)
    n = 1
    for sampdf in groupby(df, :sample)
        sample_dat = Array{Array{Float64, 1}, 1}(undef, nchannel)
        i = 1
        for chandf in groupby(sampdf, :channel) # Potential sorting bug here
            obs = copy(chandf.observation)
            sample_dat[i] = obs
            i += 1
        end
        values[n] = sample_dat
        n +=1
    end
    channellayout = MultiChannelLayout([size(c)[1] for c in values[1]])
    data = MultiChannelData(values, channellayout, nsamples)
    return data
end

"""
	varimax(A; gamma = 1.0, minit = 20, maxit = 1000, reltol = 1e-12)
VARIMAX perform varimax (or quartimax, equamax, parsimax) rotation to the column vectors of the input matrix.
# Input Arguments
- `A::Matrix{Float64}`: input matrix, whose column vectors are to be rotated. d, m = size(A).
- `gamma::Float64`: default is 1. gamma = 0, 1, m/2, and d(m - 1)/(d + m - 2), corresponding to quartimax, varimax, equamax, and parsimax.
- `minit::Int`: default is 20. Minimum number of iterations, in case of the stopping criteria fails initially.
- `maxit::Int`: default is 1000. Maximum number of iterations.
- `reltol::Float64`: default is 1e-12. Relative tolerance for stopping criteria.
# Output Argument
- `B::Matrix{Float64}`: output matrix, whose columns are already been rotated.
Implemented by Haotian Li, Aug. 20, 2019
"""
function varimax(A; gamma = 1.0, minit = 20, maxit = 1000, reltol = 1e-12)
	# Get the sizes of input matrix
	d, m = size(A)

	# If there is only one vector, then do nothing.
	if m == 1
		return A
	end

	# Warm up step: start with a good initial orthogonal matrix T by SVD and QR
	T = Matrix{Float64}(I, m, m)
	B = A * T
	L,_,M = svd(A' * (d*B.^3 - gamma*B * Diagonal(sum(B.^2, dims = 1)[:])))
	T = L * M'
	if norm(T-Matrix{Float64}(I, m, m)) < reltol
		T,_ = qr(randn(m,m)).Q
		B = A * T
	end

	# Iteration step: get better T to maximize the objective (as described in Factor Analysis book)
	D = 0
	for k in 1:maxit
		Dold = D
		L,s,M = svd(A' * (d*B.^3 - gamma*B * Diagonal(sum(B.^2, dims = 1)[:])))
		T = L * M'
		D = sum(s)
		B = A * T
		if (abs(D - Dold)/D < reltol) && k >= minit
			break
		end
	end

	# Adjust the sign of each rotated vector such that the maximum absolute value is positive.
	for i in 1:m
		if abs(maximum(B[:,i])) < abs(minimum(B[:,i]))
			B[:,i] .= - B[:,i]
		end
	end

	return B
end


text:
module CellValuesAppend

using Base: @propagate_inbounds

using Gridap
using Gridap.Helpers
using Gridap.CachedArrays

export append
export local_append

import Base: iterate
import Base: length
import Base: eltype
import Base: size
import Base: getindex
import Base: IndexStyle

function append(cvs::Vararg{IndexCellValue{T,1}}) where T
  IndexCellValueByGlobalAppend(cvs...)
end

function local_append(offset::Int,a::IndexCellVector,b::IndexCellVector)
  IndexCellVectorByLocalAppend(offset,a,b)
end

function local_append(offset::NTuple{M,Int},a::Vararg{<:IndexCellVector}) where M
  IndexCellVectorByLocalAppend(offset,a...)
end

struct IndexCellValueByGlobalAppend{
  T,V<:IndexCellValue{T,1},W<:IndexCellValue{T,1}} <: IndexCellValue{T,1}
  cvs1::V
  cvs2::W
  offset::Int
end

function IndexCellValueByGlobalAppend(cvs1::IndexCellValue{T,1}, cvs2::IndexCellValue{T,1}) where {T}
  offset = length(cvs1)
  IndexCellValueByGlobalAppend(cvs1, cvs2, offset)
end

function IndexCellValueByGlobalAppend(cvs::Vararg{IndexCellValue{T,1},N}) where {T,N}
  aux = cvs[1]
  for i in 2:N
    aux = IndexCellValueByGlobalAppend(aux, cvs[i])
  end
  return aux
end

@propagate_inbounds function getindex(self::IndexCellValueByGlobalAppend,cell::Int)
  if cell <= self.offset
    return self.cvs1[cell]
  else
    return self.cvs2[cell-self.offset]
  end
end

size(self::IndexCellValueByGlobalAppend) = (length(self.cvs1) + length(self.cvs2),)

struct IndexCellVectorByLocalAppend{
  T,V<:IndexCellVector{T},W<:IndexCellVector{T}} <: IndexCellVector{T,CachedArray{T,1,Array{T,1}},1}
  offset::Int
  cvs1::V
  cvs2::W
  cache::CachedArray{T,1,Array{T,1}}
end

function IndexCellVectorByLocalAppend(
  offset::Int,cvs1::IndexCellVector{T}, cvs2::IndexCellVector{T}) where T
  @assert length(cvs1) == length(cvs2)
  cache = CachedArray(T,1)
  IndexCellVectorByLocalAppend(offset, cvs1, cvs2, cache)
end

function IndexCellVectorByLocalAppend(
  offset_vec::NTuple{M,Int}, cvs::Vararg{<:IndexCellVector{T},N}) where {T,N,M}
  @assert M == N-1
  aux = cvs[1]
  for i in 2:N
    aux = IndexCellVectorByLocalAppend(offset_vec[i-1], aux, cvs[i])
  end
  return aux
end

@propagate_inbounds function getindex(self::IndexCellVectorByLocalAppend,cell::Integer)
  v1 = self.cvs1[cell]
  v2 = self.cvs2[cell]
  l1 = length(v1)
  l2 = length(v2)
  l = l1 + l2
  s = (l,)
  setsize!(self.cache,s)
  for i in 1:l1
    @inbounds self.cache[i] = v1[i]
  end
  for i in 1:l2
    @inbounds self.cache[i+l1] = v2[i] + self.offset
  end
  self.cache
end

length(self::IndexCellVectorByLocalAppend) = length(self.cvs1)

size(self::IndexCellVectorByLocalAppend) = (length(self),)

end # module

ext: jl


text:
using EllipsisNotation
using Base.Test

A = Array{Int}(2,4,2)

A[..,1] = [2 1 4 5
           2 2 3 6]

A[..,2] = [3 2 6 5
          3 2 6 6]

@test A[:,:,1] == [2 1 4 5
                   2 2 3 6]

@test A[:,:,2] == [3 2 6 5
                   3 2 6 6]

A[1,..] = reshape([3 4
                   5 6
                   4 5
                   6 7],1,4,2)

B = [3 4
     5 6
     4 5
     6 7]

@test B == reshape(A[1,..],4,2)


text:
using Pkg;
Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate()

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter, Zygote

makedocs(
  sitename="Zygote",
  pages = [
        "Home" => "index.md",
        "Custom Adjoints" => "adjoints.md",
        "Utilities" => "utils.md",
        "Complex Differentiation" => "complex.md",
        "Profiling" => "profiling.md",
        "Internals" => "internals.md",
        "Glossary" => "glossary.md"],
  format = Documenter.HTML(prettyurls = haskey(ENV, "CI"), analytics = "UA-36890222-9"))

deploydocs(
    repo = "github.com/FluxML/Zygote.jl.git",
)


text:
# declare the concrete subtypes of AbstractParamType and AbstractSolutionData

include("conversion_interface.jl")
include("flux_types.jl")
include("abstract_shock_sensor.jl")
include("abstract_diffusion.jl")
include("abstract_diffusion_penalty.jl")
include("shock_mesh_types.jl")
include("abstract_shock_capturing.jl")
include("shock_sensor_types.jl")
include("shock_capturing_types.jl")

@doc """
### EulerEquationMod.ParamType

  This type holds the values of any constants or parameters needed during the
  computation.  These parameters can be specified in the opts dictionary or
  have default values set here.  If there is no reasonable default, values
  are initialized to -1

  There are also a bunch of arrays that are used as temporaries by low
  level functions (to avoid having to allocate arrays themselves, which is
  a performance trap).  In general, this Type is used as a container to pass
  around values.

  gamma and R are the independent themodynamic variables

  Whether this type should be immutable or not is an open question

  This type is paramaterized on the dimension of the equation for purposes
  of multiple dispatch

  **Static Parameters**:

   * Tdim : dimensionality of the equation, integer, (used for dispatch)
   * var_type : type of variables used used in the weak form, symbol, (used for
             dispatch), currently supported values: :conservative, :entropy
   * Tsol : datatype of solution variables q
   * Tres : datatype of residual
   * Tmsh : datatype of mesh related quantities (mapping jacobian etc.)

  **Fields (with default values)**:

   * cv  : specific heat constant
   * R : specific gas constant (J/(Kg*K))
   * gamma : ratio of specific heats
   * gamma_1 : gamma - 1

  **Fields (without default values)**:

   * Ma  : free stream Mach number
   * Re  : free stream Reynolds number
   * aoa : angle of attack (radians)

"""->
mutable struct ParamType{Tdim, var_type, Tsol, Tres, Tmsh} <: AbstractParamType{Tdim}
  f::BufferedIO
  t::Float64  # current time value
  order::Int  # accuracy of elements (p=1,2,3...)

  numDofPerNode::Int
  numNodesPerElement::Int
  numNodesPerFace::Int
  #TODO: consider making these vectors views of a matrix, to guarantee
  #      spatial locality

  # numDofPerNode x stencilsize arrays for entropy variables

  # temporary storage for flux functions
  eulerfluxdata::CalcEulerFluxData{Tsol}
  bcdata::BCData{Tsol, Tres}
  roefluxdata::RoeFluxData{Tsol, Tres, Tmsh}
  calcsatdata::CalcSatData{Tres}
  lffluxdata::LFFluxData{Tres}
  irfluxdata::IRFluxData{Tsol}
  hllfluxdata::HLLFluxData{Tres}
  tworwavespeeddata::TwoRWaveSpeedData{Tsol, Tres}
  apply_entropy_kernel_diagE_data::ApplyEntropyKernel_diagEData{Tsol, Tres}
  get_lambda_max_simple_data::GetLambdaMaxSimpleData{Tsol}
  get_lambda_max_data::GetLambdaMaxData{Tsol}
  calc_vorticity_data::CalcVorticityData{Tsol, Tres, Tmsh}
  contract_res_entropy_vars_data::ContractResEntropyVarsData{Tsol}
  get_tau_data::GetTauData{Tsol, Tres}

  # temporary storage for face element integral functions
  calc_ec_face_integral_data::CalcECFaceIntegralData{Tres, Tmsh}
  calc_entropy_penalty_integral_data::CalcEntropyPenaltyIntegralData{Tsol, Tres}
  interpolate_element_staggered_data::InterpolateElementStaggeredData{Tsol}

  # temporary storage for higher level functions
  calc_volume_integrals_data::CalcVolumeIntegralsData{Tres, Tmsh}
  face_element_integral_data::FaceElementIntegralData{Tsol, Tres}
  calc_face_integrals_data::CalcFaceIntegralsData{Tsol, Tres}

  # stabilization functions
  lps_data::LPSData{Tsol, Tres}

  # entropy kernels
  entropy_lf_kernel::LFKernel{Tsol, Tres, Tmsh}
  entropy_lw2_kernel::LW2Kernel{Tsol, Tres, Tmsh}
  entropy_identity_kernel::IdentityKernel{Tsol, Tres, Tmsh}

  get_ira0data::GetIRA0Data{Tsol}

  shockmesh::ShockedElements{Tres}

  h::Float64 # temporary: mesh size metric
  cv::Float64  # specific heat constant
  R::Float64  # specific gas constant used in ideal gas law (J/(Kg * K))
  R_ND::Float64  # specific gas constant, nondimensionalized
  gamma::Float64 # ratio of specific heats
  gamma_1::Float64 # = gamma - 1

  Ma::Float64  # free stream Mach number

  # these quantities are dimensional (ie. used for non-dimensionalization)
  aoa::Tsol  # angle of attack (radians)
  sideslip_angle::Tsol
  rho_free::Float64  # free stream density
  p_free::Float64  # free stream pressure
  T_free::Float64 # free stream temperature
  E_free::Float64 # free stream energy (4th conservative variable)
  a_free::Float64 # free stream speed of sound (computed from p_free and rho_free)

  edgestab_gamma::Float64  # edge stabilization parameter
  # debugging options
  writeflux::Bool  # write Euler flux
  writeboundary::Bool  # write boundary data
  writeq::Bool # write solution variables
  use_edgestab::Bool  # use edge stabilization
  use_filter::Bool  # use filtering
  use_res_filter::Bool # use residual filtering

  filter_mat::Array{Float64, 2}  # matrix that performs filtering operation
                                 # includes transformations to/from modal representation

  use_dissipation::Bool  # use artificial dissipation
  dissipation_const::Float64  # constant used for dissipation filter matrix

  tau_type::Int  # type of tau to use for GLS stabilization

  use_Minv::Int  # apply Minv to explicit jacobian calculation, 0 = do not
                 # apply, 1 = do apply
  vortex_x0::Float64  # vortex center x coordinate at t=0
  vortex_strength::Float64  # strength of the vortex

  #TODO: get rid of these, move to NewtonData
  krylov_itr::Int  # Krylov iteration number for iterative solve
  krylov_type::Int # 1 = explicit jacobian, 2 = jac-vec prod

  Rprime::Array{Float64, 2}  # numfaceNodes x numNodesPerElement interpolation matrix
                             # this should live in sbpface instead
  x_design::Array{Tsol, 1}  # design variables

  #=
  # timings
  t_volume::Float64  # time for volume integrals
  t_face::Float64 # time for surface integrals (interior)
  t_source::Float64  # time spent doing source term
  t_sharedface::Float64  # time for shared face integrals
  t_bndry::Float64  # time spent doing boundary integrals
  t_send::Float64  # time spent sending data
  t_wait::Float64  # time spent in MPI_Wait
  t_allreduce::Float64 # time spent in allreduce
  t_barrier::Float64  # time spent in MPI_Barrier
  t_jacobian::Float64 # time spend computing Jacobian
  t_solve::Float64 # linear solve time
  =#
  time::Timings

  function ParamType{Tdim, var_type, Tsol, Tres, Tmsh}(mesh, sbp, opts, order::Integer; open_file::Bool=true) where {Tdim, var_type, Tsol, Tres, Tmsh}
  # create values, apply defaults

    # all the spatial computations happen on the *flux* grid when using
    # the staggered grid algorithm, so make the temporary vectors the
    # right size
    if opts["use_staggered_grid"]
      numNodesPerElement = mesh.mesh2.numNodesPerElement
      stencilsize = size(mesh.mesh2.sbpface.perm, 1)
      numNodesPerFace = mesh.mesh2.numNodesPerFace
    else
      numNodesPerElement = mesh.numNodesPerElement
      stencilsize = size(mesh.sbpface.perm, 1)
      numNodesPerFace = mesh.numNodesPerFace
    end

    t = 0.0
    myrank = mesh.myrank
    if DB_LEVEL >= 1 && open_file
      f = BufferedIO("log_$myrank.dat", "w")
    else
      f = BufferedIO(DevNull)
    end

    nd = 2*mesh.numDofPerNode
    eulerfluxdata = CalcEulerFluxData{Tsol}(mesh.numDofPerNode)
    bcdata = BCData{Tsol, Tres}(mesh.numDofPerNode)
    roefluxdata = RoeFluxData{Tsol, Tres, Tmsh}(mesh.numDofPerNode, mesh.dim)
    calcsatdata = CalcSatData{Tres}(mesh.numDofPerNode)
    lffluxdata = LFFluxData{Tres}(mesh.numDofPerNode, mesh.numDofPerNode)
    irfluxdata = IRFluxData{Tsol}(mesh.numDofPerNode)
    hllfluxdata = HLLFluxData{Tres}(mesh.numDofPerNode)
    tworwavespeeddata = TwoRWaveSpeedData{Tsol, Tres}(mesh.numDofPerNode)

    apply_entropy_kernel_diagE_data = ApplyEntropyKernel_diagEData{Tsol, Tres}(mesh.numDofPerNode, 2*mesh.numDofPerNode)
    get_lambda_max_simple_data = GetLambdaMaxSimpleData{Tsol}(mesh.numDofPerNode)
    get_lambda_max_data = GetLambdaMaxData{Tsol}(mesh.numDofPerNode)

    calc_vorticity_data = CalcVorticityData{Tsol, Tres, Tmsh}(mesh.numNodesPerElement, mesh.dim)
    contract_res_entropy_vars_data = ContractResEntropyVarsData{Tsol}(mesh.numDofPerNode)
    get_tau_data = GetTauData{Tsol, Tres}(mesh.numDofPerNode, mesh.dim)

    calc_ec_face_integral_data = CalcECFaceIntegralData{Tres, Tmsh}(mesh.numDofPerNode, mesh.dim)
    calc_entropy_penalty_integral_data = CalcEntropyPenaltyIntegralData{Tsol, Tres}(mesh.numDofPerNode, numNodesPerFace, stencilsize, numNodesPerElement, 2*mesh.numDofPerNode)
    interpolate_element_staggered_data = InterpolateElementStaggeredData{Tsol}(mesh.numDofPerNode, mesh.numNodesPerElement, numNodesPerElement)

    calc_volume_integrals_data = CalcVolumeIntegralsData{Tres, Tmsh}(
                                mesh.numDofPerNode, mesh.dim,
                                mesh.numNodesPerElement, numNodesPerElement, sbp)
    face_element_integral_data = FaceElementIntegralData{Tsol, Tres}(
                                 mesh.numDofPerNode, mesh.numNodesPerElement,
                                 numNodesPerElement)
    calc_face_integrals_data = CalcFaceIntegralsData{Tsol, Tres}(
                               mesh.numDofPerNode, mesh.numNodesPerFace,
                               mesh.numNodesPerElement)
    lps_data = LPSData{Tsol, Tres}(mesh, sbp, opts)
    entropy_lf_kernel = LFKernel{Tsol, Tres, Tmsh}(mesh.numDofPerNode, nd)
    entropy_lw2_kernel = LW2Kernel{Tsol, Tres, Tmsh}(mesh.numDofPerNode, mesh.dim)
    entropy_identity_kernel = IdentityKernel{Tsol, Tres, Tmsh}()
    get_ira0data = GetIRA0Data{Tsol}(mesh.numDofPerNode)

    shockmesh = ShockedElements{Tres}(mesh)

    h = maximum(mesh.jac)

    gamma = opts["gamma"]
    gamma_1 = gamma - 1
    R = opts["R"]
    cv = R/gamma_1

    Ma = opts["Ma"]
    aoa = opts[ "aoa"]*pi/180
    sideslip_angle = opts["sideslip_angle"]
    rho_free = 1.0
    p_free = opts["p_free"]
    T_free = opts["T_free"]
    E_free = p_free/gamma_1 + 0.5*rho_free*Ma*Ma
    a_free = sqrt(p_free/rho_free)  # free stream speed of sound
                                    # should be sqrt(p_free*gamma/rho_free) ?
    R_ND = R*a_free*a_free/T_free

    edgestab_gamma = opts["edgestab_gamma"]

    # debugging options
    writeflux = opts[ "writeflux"]
    writeboundary = opts[ "writeboundary"]
    writeq = opts["writeq"]
    use_edgestab = opts["use_edgestab"]
    if use_edgestab println("edge stabilization enabled") end

    use_filter = opts["use_filter"]
    if use_filter println("solution variables filter enabled") end

    use_res_filter = opts["use_res_filter"]
    if use_res_filter println("residual filter enabled") end

    if use_filter || use_res_filter || opts["use_filter_prec"]
      filter_fname = opts["filter_name"]
      filter_mat = calcFilter(sbp, filter_fname, opts)
    else
      filter_mat = zeros(Float64, 0,0)
    end

    use_dissipation = opts["use_dissipation"]
    if use_dissipation println("artificial dissipation enabled") end

    dissipation_const = opts["dissipation_const"]

    tau_type = opts["tau_type"]

    use_Minv = opts["use_Minv"] ? 1 : 0

    vortex_x0 = opts["vortex_x0"]
    vortex_strength = opts["vortex_strength"]

    krylov_itr = 0
    krylov_type = 1 # 1 = explicit jacobian, 2 = jac-vec prod

    sbpface = mesh.sbpface

    if typeof(sbpface) <: SparseFace
      Rprime = zeros(0, 0)
    else
      Rprime = zeros(size(sbpface.interp, 2), numNodesPerElement)
      # expand into right size (used in SBP Gamma case)
      for i=1:size(sbpface.interp, 1)
        for j=1:size(sbpface.interp, 2)
          Rprime[j, i] = sbpface.interp[i, j]
        end
      end
    end  # end if

    x_design = zeros(Tsol, 0)  # this can be resized later

    time = Timings()

    return new(f, t, order, mesh.numDofPerNode, mesh.numNodesPerElement,
               mesh.numNodesPerFace,
               # flux functions
               eulerfluxdata, bcdata,
               roefluxdata, calcsatdata, lffluxdata, irfluxdata, hllfluxdata,
               tworwavespeeddata,
               apply_entropy_kernel_diagE_data, get_lambda_max_simple_data,
               get_lambda_max_data,
               calc_vorticity_data, contract_res_entropy_vars_data,
               get_tau_data,
               # face level functions
               calc_ec_face_integral_data, calc_entropy_penalty_integral_data,
               interpolate_element_staggered_data,
               # entire mesh functions
               calc_volume_integrals_data,
               face_element_integral_data, calc_face_integrals_data,
               lps_data,
               entropy_lf_kernel, entropy_lw2_kernel, entropy_identity_kernel,
               get_ira0data,
               shockmesh,
               h, cv, R, R_ND, gamma, gamma_1, Ma, aoa, sideslip_angle,
               rho_free, p_free, T_free, E_free, a_free,
               edgestab_gamma, writeflux, writeboundary,
               writeq, use_edgestab, use_filter, use_res_filter, filter_mat,
               use_dissipation, dissipation_const, tau_type, use_Minv, vortex_x0,
               vortex_strength,
               krylov_itr, krylov_type,
               Rprime,
               x_design, time)

    end   # end of ParamType function

end  # end type declaration

# now that EulerData is declared, include other files that use it
@doc """
### EulerEquationMod.EulerData_

  This type is an implementation of the abstract EulerData.  It is
  parameterized by the residual datatype Tres and the mesh datatype Tmsh
  because it stores some arrays of those types.  Tres is the 'maximum' type of
  Tsol and Tmsh, where Tsol is the type of the solution variables.
  It is also paramterized by `var_type`, which should be a symbol describing
  the set of variables stored in eqn.q.  Currently supported values are
  `:conservative` and `:entropy`, which indicate the conservative variables and
  the entropy variables described in:

  'A New Finite Element Formulation for
  Computational Fluid Dynamics: Part I' by Hughes et al.`

  *Note*: this constructor does not fully populate all fields.  The
          [`init`])@ref) function must be called to finish initialization.

  **Static Parameters**:

   * Tsol : datatype of variables solution variables, ie. the
           q vector and array
   * Tres : datatype of residual. ie. eltype(res_vec)
   * Tdim : dimensionality of equation, integer, (2 or 3, currently only 2 is
           supported).
   * Tmsh : datatype of mesh related quantities
   * var_type : symbol describing variables used in weak form, (:conservative
               or :entropy)

  **Fields**

  This type has many fields, not all of them are documented here.  A few
  of the most important ones are:

   * comm: MPI communicator
   * commsize: size of MPI communicator
   * myrank: MPI rank of this process

  When computing the jacobian explicitly (options key `calc_jac_explicit`),
  Tsol and Tres are typically `Float64`, however node-level operations
  sometime use complex numbers or dual numbers.  Also, some operations on
  entropy variables require doing parts of the computation with
  conservative variables.  To support these use-cases, the fields

   * params: ParamType object with `Tsol`, `Tres`, `Tmsh`, and `var_type` matching the equation object
   * params_conservative: ParamType object with `Tsol, Tres`, and `Tmsh` matching the `EulerData_` object, but `var_type = :conservative`
   * params_entropy: similar to `param_conservative`, but `var_type = :entropy`
   * params_complex: ParamType object with `Tmsh` and `var_type` matching the `EulerData_` object, but `Tsol = Tres = Complex128`

 exist.
"""->
mutable struct EulerData_{Tsol, Tres, Tdim, Tmsh, var_type} <: EulerData{Tsol, Tres, Tdim, var_type}
# hold any constants needed for euler equation, as well as solution and data
#   needed to calculate it
# Formats of all arrays are documented in SBP.
# Only the constants are initialized here, the arrays are not.

  # this is the ParamType object that uses the same variables as
  # the EulerData_ object
  params::ParamType{Tdim, var_type, Tsol, Tres, Tmsh}
  comm::MPI.Comm
  commsize::Int
  myrank::Int

  # we include a ParamType object of all variable types, because occasionally
  # we need to do a calculation in  variables other than var_type
  # params (above) typically points to the same object as one of these
  params_conservative::ParamType{Tdim, :conservative, Tsol, Tres, Tmsh}
  params_entropy::ParamType{Tdim, :entropy, Tsol, Tres, Tmsh}

  # used for complex-stepping the boundary conditions
  # this should really be Complex{Tsol}, but that isn't allowed
  # once we switch to dual number this can be improved
  params_complex::ParamType{Tdim, :conservative, Complex128, Complex128, Tmsh}

  # the following arrays hold data for all nodes
  q::Array{Tsol,3}  # holds conservative variables for all nodes
  q_bar::Array{Tsol, 3}  # adjoint part of q
  q_face::Array{Tsol, 4}  # store solution values interpolated to faces
  q_face_bar::Array{Tsol, 4}  # adjoint part of q_face
  q_bndry::Array{Tsol, 3}  # store solution variables interpolated to
  q_bndry_bar::Array{Tsol, 3}  # adjoint part

  q_flux::Array{Tsol, 3}  # flux variable solution

  q_vec::Array{Tres,1}            # initial condition in vector form

  aux_vars::Array{Tres, 3}        # storage for auxiliary variables
  aux_vars_bar::Array{Tres, 3}    # adjoint part
  aux_vars_face::Array{Tres, 3}    # storage for aux variables interpolated
                                  # to interior faces
  aux_vars_face_bar::Array{Tres, 3}  # adjoint part
  aux_vars_sharedface::Array{Array{Tres, 3}, 1}  # storage for aux varables interpolate
                                       # to shared faces
  aux_vars_bndry::Array{Tres,3}   # storage for aux variables interpolated
                                  # to the boundaries
  aux_vars_bndry_bar::Array{Tres, 3}  # adjoint part

  # hold fluxes in all directions
  # [ndof per node by nnodes per element by num element by num dimensions]
  flux_parametric::Array{Tsol,4}  # flux in xi and eta direction
  shared_data::Array{SharedFaceData{Tsol}, 1}  # MPI send and receive buffers
  shared_data_bar::Array{SharedFaceData{Tres}, 1} # adjoint part (for eqn.q_bar
                                                  # or eqn.res_bar)

  flux_face::Array{Tres, 3}  # flux for each interface, scaled by jacobian
  res::Array{Tres, 3}             # result of computation
  res_bar::Array{Tres, 3}         # adjoint part

  res_vec::Array{Tres, 1}         # result of computation in vector form
  Axi::Array{Tsol,4}               # Flux Jacobian in the xi-direction
  Aeta::Array{Tsol,4}               # Flux Jacobian in the eta-direction
  res_edge::Array{Tres, 4}       # edge based residual used for stabilization
                           # numdof per node x nnodes per element x numEl x num edges per element

  edgestab_alpha::Array{Tmsh, 4}  # alpha needed by edgestabilization
                                  # Tdim x Tdim x nnodesPerElement x numEl
  bndryflux::Array{Tsol, 3}       # boundary flux
  stabscale::Array{Tsol, 2}       # stabilization scale factor

  # artificial dissipation operator:
  #   a square numnodes x numnodes matrix for every element
  dissipation_mat::Array{Tmsh, 3}

  Minv3D::Array{Tmsh, 3}       # inverse mass matrix for application to res, not res_vec
  Minv::Array{Tmsh, 1}         # inverse mass matrix
  M::Array{Tmsh, 1}            # mass matrix

  # TODO: consider overloading getField instead of having function as
  #       fields
  multiplyA0inv::Function         # multiply an array by inv(A0), where A0
                                  # is the coefficient matrix of the time derivative
  majorIterationCallback::Function # called before every major (Newton/RK) itr

  src_func::SRCType  # functor for the source term
  flux_func::FluxType  # functor for the face flux
  flux_func_revm::FluxType_revm # Functor for the reverse mode of face flux
  flux_func_revq::FluxType_revq
  flux_func_diff::FluxType_diff
  volume_flux_func::FluxType  # functor for the volume flux numerical flux
                              # function
  volume_flux_func_diff::FluxType_diff
  volume_flux_func_revm::FluxType_revm
  volume_flux_func_revq::FluxType_revq

  viscous_flux_func::FluxType  # functor for the viscous flux numerical flux function
  face_element_integral_func::FaceElementIntegralType  # function for face
                                                       # integrals that use
                                                       # volume data
# minorIterationCallback::Function # called before every residual evaluation
  shock_capturing::AbstractShockCapturing

  assembler::AssembleElementData  # temporary place to stash the assembler

  file_dict::Dict{String, IO}  # dictionary of all files used for logging

 # inner constructor
  function EulerData_{Tsol, Tres, Tdim, Tmsh, var_type}(mesh::AbstractMesh, sbp::AbstractOperator, opts; open_files=true) where {Tsol, Tres, Tdim, Tmsh, var_type}

    myrank = mesh.myrank
    @mpi_master begin
      println("\nConstruction EulerData object")
      println("  Tsol = ", Tsol)
      println("  Tres = ", Tres)
      println("  Tdim = ", Tdim)
      println("  Tmsh = ", Tmsh)
    end

    eqn = new()  # incomplete initialization

    eqn.comm = mesh.comm
    eqn.commsize = mesh.commsize
    eqn.myrank = mesh.myrank

    numfacenodes = mesh.numNodesPerFace

    vars_orig = opts["variable_type"]
    opts["variable_type"] = :conservative
    eqn.params_conservative = ParamType{Tdim, :conservative, Tsol, Tres, Tmsh}(
                                       mesh, sbp, opts, mesh.order)
    opts["variable_type"] = :entropy
    eqn.params_entropy = ParamType{Tdim, :entropy, Tsol, Tres, Tmsh}(
                                       mesh, sbp, opts, mesh.order;
                                       open_file=false)

    opts["variable_type"] = vars_orig
    if vars_orig == :conservative
      eqn.params = eqn.params_conservative
    elseif vars_orig == :entropy
      eqn.params = eqn.params_entropy
    else
      println(BSTDERR, "Warning: variable_type not recognized")
    end

    eqn.params_complex = ParamType{Tdim, :conservative, Complex128, Complex128, Tmsh}(
                                       mesh, sbp, opts, mesh.order;
                                       open_file=false)

    eqn.multiplyA0inv = matVecA0inv
    eqn.majorIterationCallback = majorIterationCallback

    eqn.Minv = calcMassMatrixInverse(mesh, sbp, eqn)
    eqn.Minv3D = calcMassMatrixInverse3D(mesh, sbp, eqn)
    eqn.M = calcMassMatrix(mesh, sbp, eqn)

    jac_type = opts["jac_type"]::Int
    if opts["use_dissipation"] || opts["use_dissipation_prec"]
      dissipation_name = opts["dissipation_name"]
      eqn.dissipation_mat = calcDissipationOperator(mesh, sbp, eqn, opts,
                                                    dissipation_name)
    else
      eqn.dissipation_mat = zeros(Tmsh, 0, 0, 0)
    end

    # Must initialize them because some datatypes (BigFloat)
    #   don't automatically initialize them
    # Taking a sview(A,...) of undefined values is illegal
    # I think its a bug that Array{Float64}(...) initializes values
    eqn.q = zeros(Tsol, mesh.numDofPerNode, sbp.numnodes, mesh.numEl)

    if opts["use_staggered_grid"]
      eqn.q_flux = zeros(Tsol, mesh.numDofPerNode, mesh.mesh2.numNodesPerElement, mesh.numEl)
    else
      eqn.q_flux = zeros(Tsol, 0, 0, 0)
    end

    #TODO: don't store these, recalculate as needed
    eqn.Axi = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, sbp.numnodes,
                    mesh.numEl)
    eqn.Aeta = zeros(eqn.Axi)
    eqn.aux_vars = zeros(Tsol, 1, sbp.numnodes, mesh.numEl)

    if opts["precompute_volume_flux"]
      eqn.flux_parametric = zeros(Tsol, mesh.numDofPerNode, sbp.numnodes,
                                  mesh.numEl, Tdim)
    else
      eqn.flux_parametric = zeros(Tsol, 0, 0, 0, 0)
    end

    eqn.res = zeros(Tres, mesh.numDofPerNode, sbp.numnodes, mesh.numEl)

    if opts["use_edge_res"]
      eqn.res_edge = zeros(Tres, mesh.numDofPerNode, sbp.numnodes, mesh.numEl,
                           mesh.numTypePerElement[2])
    else
      eqn.res_edge = zeros(Tres, 0, 0, 0, 0)
    end

    if mesh.isDG
      eqn.q_vec = reshape(eqn.q, mesh.numDof)
      eqn.res_vec = reshape(eqn.res, mesh.numDof)
    else
      eqn.q_vec = zeros(Tres, mesh.numDof)
      eqn.res_vec = zeros(Tres, mesh.numDof)
    end

    if opts["precompute_q_bndry"]
      eqn.q_bndry = zeros(Tsol, mesh.numDofPerNode, numfacenodes,
                                mesh.numBoundaryFaces)
    else
      eqn.q_bndry = zeros(Tsol, 0, 0, 0)
    end

    if opts["precompute_q_face"]
      eqn.q_face = zeros(Tsol, mesh.numDofPerNode, 2, numfacenodes, mesh.numInterfaces)
    else
      eqn.q_face = zeros(Tsol, 0, 0, 0, 0)
    end

    #TODO: why are there 2 if mesh.isDG blocks?
    if mesh.isDG
     if opts["precompute_face_flux"]
        eqn.flux_face = zeros(Tres, mesh.numDofPerNode, numfacenodes,
                                    mesh.numInterfaces)
      else
        eqn.flux_face = zeros(Tres, 0, 0, 0)
      end

      eqn.aux_vars_face = zeros(Tres, 1, numfacenodes, mesh.numInterfaces)
      eqn.aux_vars_bndry = zeros(Tres, 1, numfacenodes, mesh.numBoundaryFaces)
    else
      eqn.q_face = zeros(Tres, 0, 0, 0, 0)
      eqn.flux_face = zeros(Tres, 0, 0, 0)
      eqn.aux_vars_face = zeros(Tres, 0, 0, 0)
      eqn.aux_vars_bndry = zeros(Tres, 0, 0, 0)
    end

    if opts["precompute_boundary_flux"]
      eqn.bndryflux = zeros(Tsol, mesh.numDofPerNode, numfacenodes,
                            mesh.numBoundaryFaces)
    else
      eqn.bndryflux = zeros(Tsol, 0, 0, 0)
    end

    # send and receive buffers
    eqn.aux_vars_sharedface = Array{Array{Tres, 3}}(mesh.npeers)
    if mesh.isDG
      for i=1:mesh.npeers
        eqn.aux_vars_sharedface[i] = zeros(Tres, mesh.numDofPerNode,
                                        numfacenodes, mesh.peer_face_counts[i])
      end
      eqn.shared_data = getSharedFaceData(Tsol, mesh, sbp, opts,
                                          opts["parallel_data"])

    else
      eqn.shared_data = Array{SharedFaceData{Tsol}}(0)
    end

    if eqn.params.use_edgestab
      eqn.stabscale = zeros(Tres, sbp.numnodes, mesh.numInterfaces)
      eqn.edgestab_alpha = zeros(Tmsh,Tdim,Tdim,sbp.numnodes, mesh.numEl)
      calcEdgeStabAlpha(mesh, sbp, eqn)
    else
      eqn.stabscale = zeros(Tres, 0, 0)
      eqn.edgestab_alpha = zeros(Tmsh, 0, 0, 0, 0)
    end

    # functor defaults. functorThatErrors() is defined in ODLCommonTools
    eqn.flux_func = functorThatErrors()
    eqn.flux_func_revm = functorThatErrors_revm()
    eqn.volume_flux_func = functorThatErrors()
    eqn.viscous_flux_func = functorThatErrors()

    if opts["need_adjoint"]
      eqn.q_bar = zeros(eqn.q)
      eqn.q_face_bar = zeros(eqn.q_face)
      eqn.q_bndry_bar = zeros(eqn.q_bndry)

      eqn.aux_vars_bar = zeros(eqn.aux_vars)
      eqn.aux_vars_face_bar = zeros(eqn.aux_vars_face)
      eqn.aux_vars_bndry_bar = zeros(eqn.aux_vars_bndry)

      if mesh.isDG
        eqn.shared_data_bar = getSharedFaceData(Tsol, mesh, sbp, opts, PARALLEL_DATA_ELEMENT)
        #eqn.shared_data_res_bar = getSharedFaceData(Tres, mesh, sbp, opts, PARALLEL_DATA_ELEMENT)
      else
        # eqn.shared_data_res_bar = zeros(SharedFaceData, 0)
        eqn.shared_data_bar = Array{SharedFaceData{Tres}}(0)
      end

      eqn.res_bar = zeros(eqn.res)
    else  # don't allocate arrays if they are not needed
      eqn.q_bar = zeros(Tsol, 0, 0, 0)
      eqn.q_face_bar = zeros(Tsol, 0, 0, 0, 0)
      eqn.q_bndry_bar = zeros(Tsol, 0, 0, 0)

      eqn.aux_vars_bar = zeros(Tres, 0, 0, 0)
      eqn.aux_vars_face_bar = zeros(Tres, 0, 0, 0)
      eqn.aux_vars_bndry_bar = zeros(Tres, 0, 0, 0)

      #eqn.shared_data_res_bar = Array{SharedFaceData}(0)
      eqn.shared_data_bar = Array{SharedFaceData{Tres}}(0)
      eqn.res_bar = zeros(Tres, 0, 0, 0)
   end

   sensor = createShockSensor(mesh, sbp, eqn, opts)
   getShockCapturing(mesh, sbp, eqn, opts, sensor)
   eqn.assembler = NullAssembleElementData

   if open_files
     eqn.file_dict = openLoggingFiles(opts, mesh.myrank)
   else
     eqn.file_dict = Dict{String, IO}()
   end

   return eqn

  end  # end of constructor

end  # end of type declaration

"""
  Useful alias for 2D ParamType
"""
const ParamType2 = ParamType{2}

"""
  Useful alias for 3D ParamType
"""
const ParamType3 = ParamType{3}

import PDESolver: openLoggingFiles, closeLoggingFiles

"""
  This function opens all used for logging data.  In particular, every data
  file that has data appended to it in majorIterationCallback should be
  opened here.  Most files are of type BufferedIO, so they must be flushed
  periodically.

  This function requires each output to have two keys: "write_outname"
  and "write_outname_fname", where the first has a boolean value that
  controls whether or not to write the output, and the second is the
  file name (including extension) to write.

  This function contains a list of all possible log files.  Every new
  log file must be added to the list

  **Inputs**:

   * opts: options dictionary
   * myrank: MPI rank

  **Outputs**:

   * file_dict: dictionary mapping names of files to the file object
                 ie. opts["write_entropy_fname"] => f

  Exceptions: this function will throw an exception if any two file names
              are the same

  Implementation notes:
    When restarting, all files must be appended to.  Currently, files
    are appended to in all cases.
"""
function openLoggingFiles(opts, myrank::Integer)

  # output dictionary
  file_dict = Dict{AbstractString, IO}()

  # map output file names to the key name that specified them
  used_names = Dict{AbstractString, AbstractString}()

  # use the fact that the key names are formulaic
  names = ["entropy", "integralq", "kinetic_energy", "kinetic_energydt", "enstrophy", "drag"]
  @mpi_master for name in names  # only open files on the master process
    keyname = string("write_", name)
    if opts[keyname]  # if this file is being written
      fname_key = string("write_", name, "_fname")
      fname = opts[fname_key]

      if fname_key in keys(used_names)
        other_keyname = used_names[fname]
        throw(ErrorException("data file name $fname used for key $keyname is already used for key $other_keyname"))
      end

      used_names[fname] = keyname  # record this fname as used

      f = BufferedIO(opts[fname_key], "a")  # append to files (safe default)

      file_dict[fname] = f

    end  # end if
  end  # end

  return file_dict
end

function openLoggingFiles(eqn::EulerData, opts)

  # verify all existing files are closed
  for f in values(eqn.file_dict)
    @assert !isopen(f)
  end

  eqn.file_dict = openLoggingFiles(opts, eqn.myrank)

  # make all params object use the same file
  eqn.params_entropy.f = eqn.params_conservative.f
  eqn.params_complex.f = eqn.param.f

  return nothing
end

function closeLoggingFiles(eqn::EulerData, opts)

  cleanup(eqn, opts)
end

"""
  This function performs all cleanup activities before the run_physics()
  function returns.  The mesh, sbp, eqn, opts are returned by run_physics()
  so there is not much cleanup that needs to be done, mostly closing files.

  **Inputs/Outputs**:

   * eqn: the EulerData object
   * opts: the options dictionary

"""
function cleanup(eqn::EulerData, opts)

  for f in values(eqn.file_dict)
    if isopen(f)
      close(f)
    end
  end

  return nothing
end

@doc """
### EulerEquationMod.getTypeParameters

Gets the type parameters for mesh and equation objects.

**Input**

* `mesh` : Object of abstract meshing type.
* `eqn`  : Euler Equation object.

**Output**

* `Tmsh` : Type parameter of the mesh.
* `Tsol` : Type parameter of the solution array.
* `Tres` : Type parameter of the residual array.
"""->

function getTypeParameters(mesh::AbstractMesh{Tmsh}, eqn::EulerData{Tsol, Tres}) where {Tmsh, Tsol, Tres}
  return Tmsh, Tsol, Tres
end

import ODLCommonTools.getAllTypeParams

@doc """
### EulerEquationMod.getAllTypeParameters

Gets the type parameters for mesh and equation objects.

**Input**

* `mesh` : Object of abstract meshing type.
* `eqn`  : Euler Equation object.
* `opts` : Options dictionary

**Output**

* `tuple` : Tuple of type parameters. Ordering is same as that of the concrete eqn object within this physics module.

"""->
function getAllTypeParams(mesh::AbstractMesh{Tmsh}, eqn::EulerData_{Tsol, Tres, Tdim, Tmsh, var_type}, opts) where {Tmsh, Tsol, Tres, Tdim, var_type}

  tuple = (Tsol, Tres, Tdim, Tmsh, var_type)

  return tuple
end

import PDESolver.updateMetricDependents

function updateMetricDependents(mesh::AbstractMesh, sbp::AbstractOperator,
                                 eqn::EulerData, opts)

  #TODO: don't reallocate the arrays, update in place
  eqn.Minv = calcMassMatrixInverse(mesh, sbp, eqn)
  eqn.Minv3D = calcMassMatrixInverse3D(mesh, sbp, eqn)
  eqn.M = calcMassMatrix(mesh, sbp, eqn)

  if eqn.params.use_edgestab
    calcEdgeStabAlpha(mesh, sbp, eqn)
  end

  sensor = getShockSensor(eqn.shock_capturing)
  updateMetrics(mesh, sbp, opts, sensor)

  return nothing
end

import PDESolver: copyParameters

function copyParameters(eqn_old::EulerData, eqn_new::EulerData)

  copyParameters(eqn_old.params, eqn_new.params)
  copyParameters(eqn_old.params_complex, eqn_new.params_complex)
  copyParameters(eqn_old.params_conservative, eqn_new.params_conservative)
  copyParameters(eqn_old.params_entropy, eqn_new.params_entropy)
end

function copyParameters(params_old::ParamType, params_new::ParamType)

   params_new.h = params_old.h
  params_new.cv = params_old.cv
  params_new.R = params_old.R
  params_new.R_ND = params_old.R_ND
  params_new.gamma = params_old.gamma
  params_new.gamma_1 = params_old.gamma_1

  params_new.Ma = params_old.Ma

  params_new.aoa = params_old.aoa
  params_new.sideslip_angle = params_old.aoa
  params_new.rho_free = params_old.rho_free
  params_new.p_free = params_old.p_free
  params_new.T_free = params_old.T_free
  params_new.E_free = params_old.E_free
  params_new.a_free = params_old.a_free

  params_new.edgestab_gamma = params_old.edgestab_gamma

  # debugging options
  params_new.writeflux = params_old.writeflux
  params_new.writeboundary = params_old.writeboundary
  params_new.writeq = params_old.writeq
  params_new.use_edgestab = params_old.use_edgestab
  params_new.use_filter = params_old.use_filter
  params_new.use_res_filter = params_old.use_res_filter

  params_new.use_dissipation = params_old.use_dissipation
  params_new.dissipation_const = params_old.dissipation_const

  params_new.tau_type = params_old.tau_type

  params_new.use_Minv = params_old.use_Minv
  params_new.vortex_x0 = params_old.vortex_x0
  params_new.vortex_strength = params_old.vortex_strength

  resize!(params_new.x_design, length(params_old.x_design))
  copy!(params_new.x_design, params_old.x_design)

  return nothing
end


text:
typealias CudaDims Union{Integer, Tuple{Integer}, Tuple{Integer, Integer}, Tuple{Integer, Integer, Integer}}

const CudaExtent = rt.cudaExtent
CudaExtent(dims::CudaDims) = CudaExtent(get_size_x(dims), get_size_y(dims), get_size_z(dims))
CudaExtent{T}(::Type{T}, dims::CudaDims) = CudaExtent(get_size_x(dims)*sizeof(T), get_size_y(dims), get_size_z(dims))

const CudaPos = rt.cudaPos
CudaPos(dims::CudaDims) = CudaPos(get_size_x(dims)-1, get_size_y(dims)-1, get_size_z(dims)-1)
CudaPos{T}(::Type{T}, dims::CudaDims) = CudaPos((get_size_x(dims)-1)*sizeof(T), get_size_y(dims)-1, get_size_z(dims)-1)

get_size_x(g::Integer) = g
get_size_x(g::Tuple{Integer}) = g[1]
get_size_x(g::Tuple{Integer, Integer}) = g[1]
get_size_x(g::Tuple{Integer, Integer, Integer}) = g[1]

get_size_y(g::Integer) = 1
get_size_y(g::Tuple{Integer}) = 1
get_size_y(g::Tuple{Integer, Integer}) = g[2]
get_size_y(g::Tuple{Integer, Integer, Integer}) = g[2]

get_size_z(g::Integer) = 1
get_size_z(g::Tuple{Integer}) = 1
get_size_z(g::Tuple{Integer, Integer}) = 1
get_size_z(g::Tuple{Integer, Integer, Integer}) = g[3]

# dim3 and utility functions
const Dim3 = rt.dim3

Dim3(dim::CudaDims) = Dim3(get_dim_x(dim), get_dim_y(dim), get_dim_z(dim))
Dim3(dims::Integer...) = Dim3(dims)

get_dim_x(g::CudaDims) = UInt32(get_size_x(g))
get_dim_y(g::CudaDims) = UInt32(get_size_y(g))
get_dim_z(g::CudaDims) = UInt32(get_size_z(g))


text:
function applywarptoshifts(p::Vector{T},
    Bs,
    st_ind::Int,
    itp_a,
    itp_b)::Vector{T} where T <: Real

    x = similar(p)
    applywarptoshifts!(x, Bs, p, st_ind, itp_a, itp_b)

    return x
end

function applywarptoshifts!(x::Vector{T},
    Bs,
    p::Vector{T},
    st_ind::Int,
    itp_a,
    itp_b) where T <: Real

    j = st_ind - 1

    for n = 1:length(Bs)

        # first.
        i = 1
        j += 1
        x[j] = p[j]

        # itp.
        target = convertcompactdomain(p[j], -one(T), one(T), zero(T), one(T))
        a = itp_a(target)
        b = itp_b(target)

        for i = 2:length(Bs[n].ss_params.d)
            j += 1

            x[j] = MonotoneMaps.evalcompositelogisticprobit(p[j], a, b, -one(T), one(T))
        end

        for i = 1:length(Bs[n].d_singlets)
            j += 1

            x[j] = p[j]
        end
    end

    return j
end

# assume shifts always between [-1,1]
function updatemixturedwarp!(Bs::Vector{NMRSpectraSimulator.FIDModelType{T,SST}},
    p::Vector{T},
    st_ind::Int,
    fs::T,
    SW::T,
    #warp_param_set::Vector{Piecewise2DLineType{T}},
    Δsys_cs::Vector{Vector{T}},
    itp_a,
    itp_b)::Int where {T <: Real, SST}

    #@assert length(warp_param_set) == length(Δ_shifts)

    j1 = applywarptoshifts!(p, Bs, p, st_ind, itp_a, itp_b)
    j = updatemixtured!(Bs, p, st_ind, fs, SW, Δsys_cs)

    @assert j1 == j # debug.

    return j
end

function updatemixtured!(As::Vector{NMRSpectraSimulator.FIDModelType{T,NMRSpectraSimulator.SpinSysParamsType1{T}}},
    p::Vector{T},
    st_ind::Int,
    fs::T,
    SW::T,
    Δsys_cs::Vector{Vector{T}})::Int where T <: Real

    #@assert length(warp_param_set) == length(Δ_shifts)

    j = st_ind - 1

    for n = 1:length(As)

        N_spins_sys = length(As[n].ss_params.d)
        @assert length(Δsys_cs[n]) == N_spins_sys + length(As[n].d_singlets)

        for i = 1:N_spins_sys
            j += 1

            #p2 = p[j]*0.05 #*Δ_shifts[j]
            p2 = p[j]*Δsys_cs[n][i]
            As[n].ss_params.d[i] = convertΔcstoΔω0(p2, fs, SW)
        end

        for i = 1:length(As[n].d_singlets)
            j += 1

            p2 = p[j]*Δsys_cs[n][i+N_spins_sys]
            As[n].d_singlets[i] = convertΔcstoΔω0(p2, fs, SW)
        end
    end

    return j
end

function convertΔcstoΔω0(x::T, fs::T, SW::T)::T where T
    return x*2*π*fs/SW
end

function updatemixtured!(As::Vector{NMRSpectraSimulator.FIDModelType{T,NMRSpectraSimulator.SpinSysParamsType2{T}}},
    p::Vector{T},
    st_ind::Int,
    fs::T,
    SW::T,
    #warp_param_set::Vector{Piecewise2DLineType{T}},
    shift_constants::Tuple{Vector{Vector{T}},T})::Int where T <: Real

    Δsys_cs, γ = shift_constants
    @assert length(Δsys_cs) == length(As)

    j = st_ind - 1

    d_i::T = NaN

    for n = 1:length(As)

        for i = 1:length(As[n].ss_params.d)

            # first.
            j += 1
            p2 = p[j]*Δsys_cs[n][i]
            d_i = convertΔcstoΔω0(p2, fs, SW)

            As[n].ss_params.d[i][1] = d_i

            # rest.
            for k = 2:length(As[n].ss_params.d[i])

                j += 1
                As[n].ss_params.d[i][k] = p[j]*γ + d_i
            end
        end

        for i = 1:length(As[n].d_singlets)
            j += 1

            p2 = p[j]*Δ_shifts[j]
            As[n].d_singlets[i] = convertΔcstoΔω0(p2, fs, SW)
        end
    end

    return j
end

function getNd(A)::Int

    counter_sys = 0
    for i = 1:length(A.ss_params.d)
        counter_sys += length(A.ss_params.d[i])
    end

    return counter_sys + length(A.d_singlets)
end

# function getNd(A::NMRSpectraSimulator.CompoundFIDType{T,NMRSpectraSimulator.SpinSysFIDType2{T}}) where T
#
#     counter_sys = 0
#     for i = 1:length(A.ss_params.d)
#         for k = 1:length(A.ss_params.d[i])
#             counter_sys += length(A.ss_params.d[i][k])
#         end
#     end
#
#     return counter_sys + length(A.d_singlets)
# end

function updateβ!(Bs::Vector{NMRSpectraSimulator.FIDModelType{T,SST}},
    κs_β_orderings::Vector{Vector{Vector{Int}}},
    κs_β_DOFs::Vector{Vector{Int}},
    p::Vector{T},
    st_ind::Int)::Int where {T <: Real, SST}

    j = st_ind - 1

    for n = 1:length(Bs)
        for i = 1:length(Bs[n].ss_params.κs_β)

            @assert length(κs_β_orderings[n][i]) == length(Bs[n].ss_params.κs_β[i])

            for l = 1:length(Bs[n].ss_params.κs_β[i])

                idx = κs_β_orderings[n][i][l] + j # use the idx-th unique equivalent value.

                Bs[n].ss_params.κs_β[i][l] = p[idx]
            end
            j += κs_β_DOFs[n][i] # add the number of unique values for this spin system we just processed.

        end

        for i = 1:length(Bs[n].β_singlets)
            j += 1

            Bs[n].β_singlets[i] = p[j]
        end
    end

    return j
end

function updateβ!(Bs::Vector{NMRSpectraSimulator.FIDModelType{T,SST}},
    p::Vector{T},
    st_ind::Int)::Int where {T <: Real, SST}

    j = st_ind - 1

    for n = 1:length(Bs)
        for i = 1:length(Bs[n].ss_params.κs_β)
            for l = 1:length(Bs[n].ss_params.κs_β[i])

                j += 1
                Bs[n].ss_params.κs_β[i][l] = p[j]
            end
        end

        for i = 1:length(Bs[n].β_singlets)
            j += 1
            Bs[n].β_singlets[i] = p[j]
        end
    end

    return j
end

function getNβ(A::NMRSpectraSimulator.FIDModelType{T,SST}) where {T,SST}

    counter_sys = 0
    for i = 1:length(A.ss_params.κs_β)
        counter_sys += length(A.ss_params.κs_β[i])
    end

    return counter_sys + length(A.β_singlets)
end

function getNβ(κs_β_DOF::Vector{Int}, B::NMRSpectraSimulator.FIDModelType{T,SST}) where {T,SST}

    counter_sys = sum( κs_β_DOF[i] for i = 1:length(κs_β_DOF) )

    return counter_sys + length(B.β_singlets)
end

function updateλ!(Bs::Vector{NMRSpectraSimulator.FIDModelType{T,NMRSpectraSimulator.SpinSysParamsType1{T}}},
    p::Vector{T},
    st_ind::Int)::Int where T <: Real

    j = st_ind - 1

    for n = 1:length(Bs)
        for i = 1:length(Bs[n].ss_params.κs_λ)

            j += 1
            Bs[n].ss_params.κs_λ[i] = p[j]
        end

        for i = 1:length(Bs[n].κs_λ_singlets)

            j += 1
            Bs[n].κs_λ_singlets[i] = p[j]
        end
    end

    return j
end

"""
use common λ for type2, just to get code going.
"""
function updateλ!(As::Vector{NMRSpectraSimulator.FIDModelType{T,NMRSpectraSimulator.SpinSysParamsType2{T}}},
    p::Vector{T},
    st_ind::Int)::Int where T <: Real

    j = st_ind - 1

    κ_λ_i::T = NaN

    for n = 1:length(As)

        for i = 1:length(As[n].ss_params.κs_λ)

            # first.
            j += 1
            κ_λ_i = p[j]

            As[n].ss_params.κs_λ[i][1] = κ_λ_i

            # rest.
            for k = 2:length(As[n].ss_params.κs_λ[i])

                ### uncommon λ.
                #j += 1
                #As[n].ss_params.κs_λ[i][k] = κ_λ_i + p[j]*γ_λ
                ### end uncommon λ.

                ### common λ
                As[n].ss_params.κs_λ[i][k] = κ_λ_i
                ### end common λ
            end
        end

        for i = 1:length(As[n].κs_λ_singlets)

            j += 1
            As[n].κs_λ_singlets[i] = p[j]
        end
    end

    return j
end

function getNλ(A::NMRSpectraSimulator.FIDModelType{T,NMRSpectraSimulator.SpinSysParamsType1{T}}) where T

    counter_sys = 0
    for i = 1:length(A.ss_params.κs_λ)
        counter_sys += length(A.ss_params.κs_λ[i])
    end

    return counter_sys + length(A.κs_λ_singlets)
end

function getNλ(A::NMRSpectraSimulator.CompoundFIDType{T,NMRSpectraSimulator.SpinSysFIDType2{T}}) where T

    counter_sys = 0

    ### uncommon λ
    # for i = 1:length(A.ss_params.κs_λ)
    #     for k = 1:length(A.ss_params.κs_λ[i])
    #         counter_sys += length(A.ss_params.κs_λ[i])
    #     end
    # end

    ### common λ.
    for i = 1:length(A.ss_params.κs_λ)
        counter_sys += length(A.ss_params.κs_λ[i])
    end

    return counter_sys + length(A.κs_λ_singlets)
end


text:
module Clarus

#SUB-MODULES
include("Services.jl")
include("Resource.jl")
include(joinpath("categories","Ccp.jl"))
include(joinpath("categories","Compliance.jl"))
include(joinpath("categories","Credit.jl"))
include(joinpath("categories","Dates.jl"))
include(joinpath("categories","Etd.jl"))
include(joinpath("categories","Frtb.jl"))
include(joinpath("categories","Hedge.jl"))
include(joinpath("categories","Load.jl"))
include(joinpath("categories","Margin.jl"))
include(joinpath("categories","Market.jl"))
include(joinpath("categories","Mifid.jl"))
include(joinpath("categories","Portfolio.jl"))
include(joinpath("categories","ProfitLoss.jl"))
include(joinpath("categories","Risk.jl"))
include(joinpath("categories","Sdr.jl"))
include(joinpath("categories","Simm.jl"))
include(joinpath("categories","Trade.jl"))
include(joinpath("categories","Util.jl"))
include(joinpath("categories","Xva.jl"))
#

using Clarus.Services: api_key, api_request, api_secret, Response, dataframe!,drilldown,pivot
using Clarus.Resource: read,write

end #MODULE-END


text:
using CUDAnative, GPUArrays, CuArrays
using CUDAnative: device!, devices, CuDevice
using CuYao
import CuYao: cu
CuArrays.allowscalar(false)

function cu(chem::QuantumMPS)
    QuantumMPS(chem.nbit_measure, chem.nbit_virtual, chem.nbit_ancilla, chem.circuit, chem.initial_reg |> cu, chem.input_state)
end

max_stars_repo_name: Roger-luo/QuantumMPS

max_issues_repo_name: Roger-luo/QuantumMPS

max_forks_repo_name: Roger-luo/QuantumMPS


text:
function catchfire(i, j, p, forest)
    countlateral = 0 # Count the buring trees around point i, j
    countdiagonal = 0
    for k in (i-1):(i+1), l in (j-1):(j+1)
        if forest[k, l] == 2
            if k != i && l != j
                countdiagonal += 1
            else
                countlateral += 1
            end
        end
    end
    fire = rand() >= ((1-p)^countlateral)*(1-0.7*p)^countdiagonal
    return fire
end

function burn(forest, p)
    nx = size(forest, 1)
    ny = size(forest, 2)
    new_forest = fill(0, (nx, ny))

    for i in 1:nx, j in 1:ny
        if forest[i, j] == 3
            new_forest[i, j] = 3 - catchfire(i, j, p, forest)
        elseif forest[i, j] == 2
            new_forest[i, j] = 1
        else
            new_forest[i, j] = forest[i, j]
        end
    end
    return new_forest
end

max_stars_repo_name: TheSquareRoot/ForestFire

max_issues_repo_name: TheSquareRoot/ForestFire

max_forks_repo_name: TheSquareRoot/ForestFire


text:
@enum TerminationStatus OPTIMAL INFEASIBLE INF_OR_UNBOUNDED EARLY_TERMINATION NOT_OPTIMIZED

mutable struct Result
    primal_bound::Float64
    dual_bound::Float64
    best_solution::Vector{Float64}
    termination_status::TerminationStatus
    total_node_count::Int
    total_simplex_iters::Int
    total_elapsed_time_sec::Float64
    total_model_builds::Int
    total_warm_starts::Int

    function Result()
        return new(Inf, -Inf, Float64[], NOT_OPTIMIZED, 0, 0, 0, 0, 0)
    end
end

function Result(state::CurrentState, config::AlgorithmConfig)
    update_dual_bound!(state)
    result = Result()
    result.primal_bound = state.primal_bound
    result.dual_bound = state.dual_bound
    copy!(result.best_solution, state.best_solution)
    if state.primal_bound == state.dual_bound == Inf
        result.termination_status = INFEASIBLE
    elseif state.primal_bound == -Inf
        result.termination_status = INF_OR_UNBOUNDED
    elseif _optimality_gap(state.primal_bound, state.dual_bound) <=
           config.gap_tol
        result.termination_status = OPTIMAL
    else
        result.termination_status = EARLY_TERMINATION
    end
    result.total_node_count = state.total_node_count
    result.total_simplex_iters = state.total_simplex_iters
    result.total_elapsed_time_sec = state.total_elapsed_time_sec
    result.total_model_builds = state.total_model_builds
    result.total_warm_starts = state.total_warm_starts
    return result
end

@enum NodeResultStatus NOT_SOLVED PRUNED_BY_PARENT_BOUND OPTIMAL_LP INFEASIBLE_LP UNBOUNDED_LP

mutable struct NodeResult
    status::NodeResultStatus
    cost::Float64
    x::Vector{Float64}
    simplex_iters::Int
    depth::Int
    int_infeas::Int

    function NodeResult(node::Node)
        node_result = new()
        node_result.status = NOT_SOLVED
        node_result.cost = NaN
        node_result.simplex_iters = 0
        node_result.depth = node.depth
        node_result.int_infeas = 0
        return node_result
    end
    function NodeResult(
        status::NodeResultStatus,
        cost::Float64,
        x::Vector{Float64},
        simplex_iters::Int,
        depth::Int,
        int_infeas::Int,
    )
        return new(status, cost, x, simplex_iters, depth, int_infeas)
    end
end


text:
simple = @stochastic_model begin
    @stage 1 begin
        @decision(model, x₁ >= 40)
        @decision(model, x₂ >= 20)
        objective = @expression(model, 100*x₁ + 150*x₂)
        @objective(model, Min, objective)
        @constraint(model, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @uncertain q₁ q₂ d₁ d₂
        @variable(model, 0 <= y₁ <= d₁)
        @variable(model, 0 <= y₂ <= d₂)
        @objective(model, Max, q₁*y₁ + q₂*y₂)
        @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end

ξ₁ = Scenario(q₁ = 24.0, q₂ = 28.0, d₁ = 500.0, d₂ = 100.0, probability = 0.4)
ξ₂ = Scenario(q₁ = 28.0, q₂ = 32.0, d₁ = 300.0, d₂ = 300.0, probability = 0.6)

simple_res = SPResult([46.67,36.25], -855.83, -1518.75, 662.92, 286.92, -1445.92, -568.92)
push!(problems, (simple, [ξ₁,ξ₂], simple_res, "Simple"))


text:
import Base.length
Base.length(ps::Particles) = length(ps.calls)

function initialize_filter(observations::K,
                           num_particles::Int,
                           fn::Function,
                           args::Tuple) where K <: ConstrainedSelection
    ps, _ = importance_sampling(observations, num_particles, fn, args)
    return ps
end

function filter_step!(observations::K,
                      ps::Particles,
                      argdiffs::D,
                      new_args::Tuple) where {K <: ConstrainedSelection, D <: Diff}
    num_particles = length(ps)
    Threads.@threads for i in 1:num_particles
        _, ps.calls[i], uw, _, _ = update(observations, ps.calls[i], argdiffs, new_args...)
        ps.lws[i] += uw
    end
end

function filter_step!(observations::K,
                      ps::Particles,
                      argdiffs::D,
                      new_args::Tuple,
                      proposal::Function,
                      proposal_args::Tuple) where {K <: ConstrainedSelection, D <: Diff}
    num_particles = length(ps)
    Threads.@threads for i in 1:num_particles
        _, p_cl, p_w = propose(proposal, ps.calls[i], proposal_args...)
        sel = selection(p_cl)
        merge!(sel, observations)
        _, ps.calls[i], u_w, _, _ = update(sel, ps.calls[i], argdiffs, new_args...)
        ps.lws[i] += u_w - p_w
    end
end

function check_ess_resample!(ps::Particles)
    num_particles = length(ps.calls)
    ltw, lnw = nw(ps.lws)
    if ess(lnw) < length(ps) / 2
        weights = exp.(lnw)
        ps.lmle += ltw - log(num_particles)
        selections = rand(Categorical(weights/sum(weights)), num_particles)
        calls = map(selections) do ind
            ps.calls[ind]
        end
        ps.calls = calls
        ps.lws = zeros(num_particles)
        return true
    end
    return false
end

function get_lmle(ps::Particles)
    return ps.lmle + lse(ps.lws) - log(length(ps))
end


text:
"""
    AbstractDVec{K,V}

Abstract type for data structures that behave similar to sparse vectors, but are indexed
by an arbitrary type `V` (could be non-integers) similarly to dictionaries. `AbstractDVec`s
are  designed to work well with [`lomc!`](@ref) and
[KrylovKit](https://github.com/Jutho/KrylovKit.jl).

Concrete implementations are available as [`DVec`](@ref)
and [`InitiatorDVec`](@ref).

`AbstractDvec`s lie somewhere between `AbstractDict`s and sparse `AbstractVector`s, while
being subtyped to neither.
Generally they behave
like a dictionary, while supportting various linear algebra functionality. Indexing with a
value not stored in the dictionary returns `zero(V)`. Setting a stored value to 0 or below
`eps(V::AbstractFloat)` removes the value from the dictionary. Their `length` signals the
number of stored elements, not the size of the vector space.

They have a [`StochasticStyle`](@ref) which selects the spawning algorithm in `FCIQMC`.

To iterate over an `AbstractDVec`, use `keys`, `pairs`, or `values`.

# Interface

The interface is similar to the `AbstractDict` interface, but with the changed behaviour
as noted above.
Implement what would be needed for the `AbstractDict` interface (`pairs`, `keys`, `values`,
`setindex!, getindex, delete!, length, haskey, empty!, isempty`) and, in addition:
* [`StochasticStyle`](@ref)
* [`storage`](@ref) returns an `AbstractDict` storing the raw data with possibly
  different `valtype` than `V`.
* [`deposit`](@ref)

See also [`DictVectors`](@ref), [`Interfaces`](@ref).
"""
abstract type AbstractDVec{K,V} end

"""
    deposit!(w::AbstractDVec, add, val, parent::Pair)

Add `val` into `w` at address `add`, taking into account initiator rules if applicable.
`parent` contains the `address => value` pair from which the pair `add => val`
was created. [`InitiatorDVec`](@ref) can intercept this and add its own functionality.
"""
function deposit!(w, add, val, _)
    w[add] += convert(valtype(w), val)
end

"""
    localpart(dv) -> AbstractDVec

Get the part of `dv` that is located on this MPI rank. Returns `dv` itself for
`AbstractDVec`s.
"""
localpart(dv) = dv # default for local data

"""
    storage(dvec) -> AbstractDict

Return the raw storage associated with `dvec` as an `AbstractDict`. Used in MPI
communication.
"""
storage(v::AbstractVector) = v

StochasticStyle(::AbstractArray{T}) where {T} = default_style(T)

"""
    freeze(dv)

Create a "frozen" version of `dv` which can no longer be modified or used in the
conventional manner, but supports faster dot products.

If `dv` is an [`MPIData`](@ref), synchronize its contents among the ranks first.
"""
freeze(v::AbstractVector) = copy(v)


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 04/16/2014
# Renamed from DeterministicPR and added changes to match LL pilot response/dynamics

module LLDetPRImpl

export
    initialize,
    update,

    updatePilotResponse,

    LLDetPR,
    LLDetPRAircraft,
    LLDetPRCommand,
    LLDetPRRA,
    LLDetPROutput

using AbstractPilotResponseImpl
using AbstractPilotResponseInterfaces
using CommonInterfaces

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractPilotResponseInterfaces.updatePilotResponse

using DataStructures

import Base: isequal, ==, empty!

type LLDetPRAircraft
    vh::Float64             #vertical rate of aircraft (ft/s)
end

type LLDetPRCommand
    t::Float64
    v_d::Float64             #commanded acceleration (ft/s^2)
    h_d::Float64             #commanded vertical rate (ft/s)
    psi_d::Float64           #commanded turn rate (deg/s)
end

type LLDetPROutput
    t::Float64
    v_d::Float64             #commanded acceleration (ft/s^2)
    h_d::Float64             #commanded vertical rate (ft/s)
    psi_d::Float64           #commanded turn rate (deg/s)
    dh_min::Float64          #dh_min of advisory being followed (ft/s)
    dh_max::Float64          #dh_max of advisory being followed (ft/s)
    target_rate::Float64     #target_rate of advisory being followed (ft/s)
    ddh::Float64             #ddh of advisory being followed (ft/s^2)
    logProb::Float64 #log probability of generating this command
end

LLDetPROutput() = LLDetPROutput(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

type LLDetPRRA
    dh_min::Float64
    dh_max::Float64
    target_rate::Float64
    ddh::Float64
end

type QueueEntry
    t::Int64          #time left before reaching front of queue (s)
    RA::LLDetPRRA     #associated RA
end

type LLDetPR <: AbstractPilotResponse
    initial_resp_time::Int64 #initial response delay
    subsequent_resp_time::Int64 #subsequent response delay
    limit_hd::Bool #limit pilot command rate to (dh_min,dh_max)
    state::Symbol #{none,follow,stay}
    queue::Vector{QueueEntry} #TODO: replace this with Queue
    timer::Int64
    COC_RA::LLDetPRRA
    output::LLDetPROutput #preallocate
    vh_history::Queue{Float64}
    history_length::Int64
    accel_thresh::Float64 #accel above which counts as too much change, in ft/s^2
    accel_near_RA::Bool #flag to indicate that accel thresh was exceeded near initial RA

    function LLDetPR(initial_resp_time::Int64, subsequent_resp_time::Int64; limit_hd::Bool=true,
        history_length::Int64=5, accel_thresh::Float64=5.0)
        obj = new()
        obj.initial_resp_time = initial_resp_time
        obj.subsequent_resp_time = subsequent_resp_time
        obj.limit_hd = limit_hd;
        obj.state = :none
        obj.queue = QueueEntry[]
        obj.timer = 0
        obj.COC_RA = LLDetPRRA(-9999.0, 9999.0, 0.0, 0.0) #COC TODO: Remove hardcoding of 9999.0, or at least centralize it
        obj.output = LLDetPROutput()
        obj.vh_history = Queue(Float64)
        obj.history_length = history_length
        obj.accel_thresh = accel_thresh
        obj.accel_near_RA = false
        obj
    end
end

==(x::LLDetPRRA, y::LLDetPRRA) = isequal(x, y)
function isequal(x::LLDetPRRA, y::LLDetPRRA)
  x.dh_min == y.dh_min && x.dh_max == y.dh_max && x.target_rate == y.target_rate
end

function add_to_queue!(q::Vector{QueueEntry}, RA::LLDetPRRA,
                       queuetime::Int64)
  if queuetime < 0
    return
  end
  filter!(x -> x.t < queuetime, q) #all elements should have a time smaller than what's being added
  el = QueueEntry(queuetime, RA)
  push!(q, el)
end

isfirstRA(pr::LLDetPR) = isequal(pr.queue[1].RA, pr.COC_RA) && length(pr.queue) == 1

function updatePilotResponse(pr::LLDetPR, command::LLDetPRCommand, RA::LLDetPRRA,
    ac_state::LLDetPRAircraft)
    t, v_d, h_d, psi_d = command.t, command.v_d, command.h_d, command.psi_d

    @assert RA.dh_min <= RA.target_rate <= RA.dh_max #this should always hold

    #decrement timer for all ra's in queue
    for i = 1:endof(pr.queue)
        pr.queue[i].t = max(0, pr.queue[i].t - 1) #TODO: remove hardcoding of t-1 for step size
    end

    #incorporate the new RA if it's new.  Could be coc
    if !isequal(pr.queue[end].RA, RA)
        if isequal(RA, pr.COC_RA)
            queuetime = 0 #COC's happen immediately
        elseif isfirstRA(pr)
            queuetime = pr.initial_resp_time
        else
            queuetime = pr.subsequent_resp_time
        end
        add_to_queue!(pr.queue, RA, queuetime)
    end

    #shift items to the next one with time = 0
    last_index = findlast(map(x->x.t, pr.queue), 0)
    splice!(pr.queue, 1:(last_index - 1))

    @assert !isempty(pr.queue) #queue should never be empty
    @assert pr.queue[1].t == 0 #first item should always have its time expired

    dh_min      = pr.queue[1].RA.dh_min
    dh_max      = pr.queue[1].RA.dh_max
    target_rate = pr.queue[1].RA.target_rate
    ddh         = pr.queue[1].RA.ddh

    #keep a sliding history of vh
    enqueue!(pr.vh_history, ac_state.vh)
    while length(pr.vh_history) > pr.history_length
        dequeue!(pr.vh_history)
    end
    #reset at each time step so that flag is only true once if true at all
    pr.accel_near_RA = false

    if isequal(pr.queue[1].RA, pr.COC_RA) #currently COC
        pr.state = :none

        #restrict pilot to not accelerate against any pending RA's
        if pr.limit_hd && length(pr.queue) > 1 #currently COC, but RA is queued
            h_d = pr.queue[2].RA.target_rate > ac_state.vh ?
                max(ac_state.vh, h_d) : min(ac_state.vh, h_d)
        end
        #Check history for accelerations exceeding thresh
        #currently COC, but RA is queued
        #only check at time of initial RA
        if pr.history_length > 0 && length(pr.queue) > 1 && pr.queue[2].t == pr.initial_resp_time
            pr.accel_near_RA = exceeds_thresh(pr.vh_history, pr.accel_thresh)
        end
    else #RA
        #perfect compliance
        #let LLADM (scripted_dynamics) handle the RA
        #following based on dh_min, dh_max, ddh
        #don't interfere with resolve_TCAS_and_script in LLADM
        h_d = ac_state.vh

        if length(pr.queue) > 1
            #currently following an RA and there is another RA queued
            pr.state = :stay
        else
            #currently following an RA and there are no other RA's queued
            pr.state = :follow
        end
    end

    # if there is a next item, show its timer, otherwise 0
    pr.timer = length(pr.queue) > 1 ? pr.queue[2].t : 0

    pr.output.t = t
    pr.output.v_d = v_d
    pr.output.h_d = h_d
    pr.output.psi_d = psi_d
    pr.output.dh_min = dh_min
    pr.output.dh_max = dh_max
    pr.output.target_rate = target_rate
    pr.output.ddh = ddh
    pr.output.logProb = 0.0 #probability = 1

    pr.output
end

function exceeds_thresh(q::Queue{Float64}, thresh::Float64)
    x = first(q)
    for y in q
        x = abs(x-y) <= thresh ? y : return true
    end
    false
end

function update(pr::LLDetPR, command, RA, ac_state)
    update(pr, convert(LLDetPRCommand, command), convert(LLDetPRRA, RA),
        convert(LLDetPRAircraft, ac_state))
end
function update(pr::LLDetPR, command::LLDetPRCommand, RA::LLDetPRRA, ac_state::LLDetPRAircraft)
    updatePilotResponse(pr, command, RA, ac_state)
end

function empty!{T}(q::Queue{T})
    while length(q) > 0
        dequeue!(q)
    end
end

function initialize(pr::LLDetPR)
    pr.state = :none
    empty!(pr.queue)
    add_to_queue!(pr.queue, pr.COC_RA, 0)
    empty!(pr.vh_history)

    pr.output.t = 0.0
    pr.output.v_d = 0.0
    pr.output.h_d = 0.0
    pr.output.psi_d = 0.0
    pr.output.dh_min = -9999.0
    pr.output.dh_max = 9999.0
    pr.output.target_rate = 0.0
    pr.output.ddh = 0.0
    pr.output.logProb = 0.0

    nothing
end

end #module

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
"""
    gridwalk(nrows::Integer, ncols::Integer, ntimesteps::Integer,
             initialvertex=rand(1:nrows*ncols)) -> Vector{Int}
"""
function gridwalk(nrows::Integer, ncols::Integer, ntimesteps::Integer,
                  initialvertex=rand(1:nrows*ncols))
    rw = RandomWalk(grid((nrows, ncols)), initialvertex)
    return rand(rw, ntimesteps)
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
include("../CovarianceFunctions/CovarianceFunctions.jl")

# For linear algebra
using LinearAlgebra
# logging
using Logging
using PyPlot
## PDEs type
abstract type AbstractPDEs end
struct NonlinElliptic2d{Tα,Tm,TΩ} <: AbstractPDEs
    # eqn: -Δu + α*u^m = f in [Ω[1,1],Ω[2,1]]*[Ω[1,2],Ω[2,2]]
    α::Tα
    m::Tm
    Ω::TΩ
    bdy::Function
    rhs::Function
end

## sample points
function sample_points_rdm(eqn::NonlinElliptic2d, N_domain, N_boundary)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]

    X_domain = hcat(rand(Float64,(N_domain, 1))*(x1r-x1l).+x1l,rand(Float64,(N_domain, 1))*(x2r-x2l).+x2l)

    N_bd_each=convert(Int64, N_boundary/4)
    if N_boundary != 4 * N_bd_each
        println("[sample points] N_boundary not divided by 4, replaced by ", 4 * N_bd_each)
        N_boundary = 4 * N_bd_each
    end

    X_boundary = zeros((N_boundary, 2))
    # bottom face
    X_boundary[1:N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2l*ones(N_bd_each))
    # right face
    X_boundary[N_bd_each+1:2*N_bd_each, :] = hcat(x1r*ones(N_bd_each),(x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    # top face
    X_boundary[2*N_bd_each+1:3*N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2r*ones(N_bd_each))

    # left face
    X_boundary[3*N_bd_each+1:N_boundary, :] = hcat(x1l*ones(N_bd_each), (x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    return X_domain', X_boundary'
end

function sample_points_grid(eqn::NonlinElliptic2d, h_in, h_bd)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]
    x = x1l + h_in:h_in:x1r-h_in
    y = x2l + h_in:h_in:x2r-h_in
    X_domain = reduce(hcat,[[x[i], y[j]] for i in 1:length(x) for j in 1:length(x)])

    l = length(x1l:h_bd:x1r-h_bd)
    X_boundary = vcat([x1l:h_bd:x1r-h_bd x2l*ones(l)], [x1r*ones(l) x2l:h_bd:x2r-h_bd], [x1r:-h_bd:x1l+h_bd x2r*ones(l)], [x1l*ones(l) x2r:-h_bd:x1l+h_bd])
    return X_domain, X_boundary'
end

## exact algorithms
# assemby Gram matrices
function get_Gram_matrices(eqn::NonlinElliptic2d, cov::AbstractCovarianceFunction, X_domain, X_boundary, sol_now)
    d = 2
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    Δδ_coefs = -1.0
    δ_coefs_int = eqn.α*eqn.m*(sol_now.^(eqn.m-1))

    # get linearized PDEs correponding measurements
    meas_δ = [PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_Δδ = [ΔδPointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), Δδ_coefs, δ_coefs_int[i]) for i = 1:N_domain]
    meas_test_int = [PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]

    Theta_train = zeros(N_domain+N_boundary,N_domain+N_boundary)

    measurements = Vector{Vector{<:AbstractPointMeasurement}}(undef,2)
    measurements[1] = meas_δ; measurements[2] = meas_Δδ
    cov(Theta_train, reduce(vcat,measurements))

    Theta_test = zeros(N_domain+N_boundary,N_domain+N_boundary)
    measurements_δ = Vector{Vector{<:AbstractPointMeasurement}}(undef,2)
    measurements_δ[1] = meas_δ; measurements_δ[2] = meas_test_int

    cov(Theta_test, reduce(vcat,measurements_δ), reduce(vcat,measurements))
    # cov(view(Theta_test,1:N_domain,1:N_boundary), meas_test_int, meas_δ)
    # cov(view(Theta_test,1:N_domain,N_boundary+1:N_domain+N_boundary), meas_test_int, meas_Δδ)
    return Theta_train, Theta_test
end

function  get_initial_covariance(cov::AbstractCovarianceFunction, X_domain, X_boundary)
    d = 2
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    meas_δ_bd = [PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_bd_int = [PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]
    Theta_initial = zeros(N_domain+N_boundary,N_domain+N_boundary)
    measurements = Vector{Vector{<:AbstractPointMeasurement}}(undef,2)
    measurements[1] = meas_δ_bd; measurements[2] = meas_bd_int
    cov(Theta_initial, reduce(vcat,measurements))
    return Theta_initial
end

# iterative GPR
function iterGPR_exact(eqn, cov, X_domain, X_boundary, sol_init, noise, GNsteps)
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    # get the rhs and bdy data
    rhs = [eqn.rhs(X_domain[:,i]) for i in 1:N_domain] .+ sqrt(noise)*randn(N_domain)
    bdy = [eqn.bdy(X_boundary[:,i]) for i in 1:N_boundary]
    v = sol_init

    Theta_train = zeros(N_domain+N_boundary, N_domain+N_boundary)
    Theta_test = zeros(N_domain,N_domain+N_boundary)
    for _ in 1:GNsteps
        Theta_train, Theta_test = get_Gram_matrices(eqn, cov, X_domain, X_boundary, v)
        rhs_now = vcat(bdy, rhs.+eqn.α*(eqn.m-1)*v[N_boundary+1:end].^eqn.m)

        v = Theta_test*(((Theta_train+noise*LinearAlgebra.I))\rhs_now)
    end

    Cov = get_initial_covariance(cov, X_domain, X_boundary)
    Cov_posterior = Cov .- Theta_test*(((Theta_train+noise*LinearAlgebra.I))\Theta_test')
    return v[N_boundary+1:end], Cov_posterior
end

α = 10.0
m = 3
Ω = [[0,1] [0,1]]
h_in = 0.02
h_bd = 0.02
lengthscale = 0.3
kernel = "Matern5half"
cov = MaternCovariance5_2(lengthscale)
noise = 0.0001
GNsteps = 3

# ground truth solution
freq = 1000
s = 4
function fun_u(x)
    ans = 0
    @inbounds for k = 1:freq
        ans += sin(pi*k*x[1])*sin(pi*k*x[2])/k^s
        # H^t norm squared is sum 1/k^{2s-2t}, so in H^{s-1/2}
    end
    return ans
end

# right hand side
function fun_rhs(x)
    ans = 0
    @inbounds for k = 1:freq
        ans += (2*k^2*pi^2)*sin(pi*k*x[1])*sin(pi*k*x[2])/k^s
    end
    return ans + α*fun_u(x)^m
end

# boundary value
function fun_bdy(x)
    return fun_u(x)
end

@info "[solver started] NonlinElliptic2d"
@info "[equation] -Δu + $α u^$m = f"
eqn = NonlinElliptic2d(α,m,Ω,fun_bdy,fun_rhs)

# N_domain = 2000
# N_boundary = 400
# X_domain, X_boundary = sample_points_rdm(eqn,N_domain, N_boundary)

X_domain, X_boundary = sample_points_grid(eqn, h_in, h_bd)
N_domain = size(X_domain,2)
N_boundary = size(X_boundary,2)
@info "[sample points] grid size $h_in"
@info "[sample points] N_domain is $N_domain, N_boundary is $N_boundary"
@info "[kernel] choose $kernel, lengthscale $lengthscale\n"
@info "[noise] $noise"
@info "[GNsteps] $GNsteps"

sol_init = zeros(N_domain+N_boundary) # initial solution
truth = [fun_u(X_domain[:,i]) for i in 1:N_domain]
@time sol_exact, sol_postvar = iterGPR_exact(eqn, cov, X_domain, X_boundary, sol_init, noise, GNsteps)
pts_accuracy_exact = sqrt(sum((truth-sol_exact).^2)/sum(truth.^2))
@info "[L2 accuracy: exact method] $pts_accuracy_exact"
pts_max_accuracy_exact = maximum(abs.(truth-sol_exact))/maximum(abs.(truth))
@info "[Linf accuracy: exact method] $pts_max_accuracy_exact"

sol_std = [sqrt(sol_postvar[i,i]) for i in N_boundary+1:N_boundary+N_domain]

Nh = convert(Int,sqrt(N_domain))
figure()
plot_surface(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(truth,Nh,Nh), label = "Reference")
plot_surface(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(sol_exact - sol_std,Nh,Nh), color="C1", label = "Lower")
plot_surface(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(sol_exact + sol_std,Nh,Nh), color="C1", label = "Upper")
display(gcf())

figure()
contourf(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(abs.(truth-sol_exact),Nh,Nh))
colorbar()
display(gcf())

figure()
contourf(reshape(X_domain[1,:],Nh,Nh), reshape(X_domain[2,:],Nh,Nh), reshape(sol_std,Nh,Nh))
colorbar()
display(gcf())


text:
using DelimitedFiles, FileIO

function save_training(filename, qopt, loss::Vector, params::Vector)
    save(filename, "qopt", qopt, "loss", loss, "params", params)
end

function load_training(filename)
    res = load(filename)
    #res["qopt"], res["loss"], res["params"]
end

function decode(nx::Int, ny::Int, depth::Int=5, nv::Int=1)
    suff = "data/j1j2-nx$nx-ny$ny-nv$nv-d$depth"
    res = load_training("$suff.jld2")
    writedlm("$suff-params.dat", res["params"])
    writedlm("$suff-loss.dat", res["loss"])
end


text:
export
    # cudnnSoftmaxMode_t
    CUDNN_SOFTMAX_MODE_INSTANCE,
    CUDNN_SOFTMAX_MODE_CHANNEL,

    # cudnnSoftmaxAlgorithm_t
    CUDNN_SOFTMAX_FAST,
    CUDNN_SOFTMAX_ACCURATE,
    CUDNN_SOFTMAX_LOG

function softmax{T}(x::CuArray{T}; algo=CUDNN_SOFTMAX_ACCURATE, mode=CUDNN_SOFTMAX_MODE_CHANNEL)
    h = handle(x)
    dims1 = ntuple(_ -> 1, 4-ndims(x))
    x4d = reshape(x, dims1..., size(x)...)
    xdesc = TensorDesc(x4d)
    y = similar(x)
    cudnnSoftmaxForward(h, algo, mode, T[1], xdesc, x, T[0], xdesc, y)

    function backward!(dy, dx)
        isvoid(dx) && return
        cudnnSoftmaxBackward(h, algo, mode, T[1], xdesc, y, xdesc, dy, T[1], xdesc, dx)
    end
    y, backward!
end

logsoftmax(x) = softmax(x, algo=CUDNN_SOFTMAX_LOG)


text:
using PyCall
using Cambrian
import Distances
import Random
import Formatting
import Base.GC

export GymEnv, fitness, new_gen!, gym_envs

import RDatasets

"Gym BERL environment"
mutable struct GymEnv <: BERLenv
    name::String
    memory::Array # Array of storable data
    gym
    env_name
    gen
    n_steps::Int64
end

environments["gym"]=GymEnv
gym_envs = ["Acrobot-v1", "CartPole-v1", "MountainCar-v0", "MountainCarContinuous-v0", "Pendulum-v0"]

function GymEnv(cfg::Dict)
    env_name = get!(cfg, "gym_env", "MountainCar-v0") # Default mountain car
    cfg["env"]="Gym-$env_name"

    gym = pyimport("gym")
    game = gym.make(env_name)
    cfg["n_out"] = game.action_space.n # length(env.action_space.sample())
    cfg["n_in"] = length(game.observation_space.sample())

    gen = get!(cfg, "seed", 0)
    Random.seed!(gen)

    GymEnv(cfg["env"], [], gym, env_name, gen, 0)
end

"Returns the fitness of one Cambrian Individual in this environment"
function fitness(indiv::Cambrian.Individual, env::GymEnv)
    game = env.gym.make(env.env_name)
    game.seed(env.gen)
    obs = game.reset()
    total_reward = 0.0
    done = false

    max_obs = Float64(max(-minimum(game.observation_space.low),
                          maximum(game.observation_space.high)))

    while ~done
        action_vect = process(indiv, obs ./ max_obs)
        # println(action_vect)
        action = argmax(action_vect)-1
        obs, reward, done, _ = game.step(action)
        total_reward += reward
        env.n_steps += 1
    end

    game.close()
    game = nothing
    Base.GC.gc()
    [total_reward]
end

"Sets up the env for the next generation"
function new_gen!(env::GymEnv, cfg::Dict=nothing)
    env.gen += 1
end


text:
module Freestreams

export Freestream

using ..Elements
using ..Motions

struct Freestream <: Element
    U::ComplexF64
end

Elements.@kind Freestream Singleton

Motions.induce_velocity(::ComplexF64, f::Freestream, t) = f.U
Motions.induce_velocity(f::Freestream, src, t) = nothing
Motions.induce_velocity!(vel, f::Freestream, src, t) = nothing

Elements.streamfunction(z::ComplexF64, f::Freestream) = -imag(f.U*z')
Elements.complexpotential(z::ComplexF64, f::Freestream) = f.U'*z

Elements.circulation(::Freestream) = 0.0
Elements.flux(::Freestream) = 0.0

Motions.self_induce_velocity!(vel, f::Freestream, t) = nothing

Motions.allocate_velocity(::Freestream) = ComplexF64

end


text:
try Pkg.clone("https://github.com/zsunberg/LaserTag.jl.git") end
Pkg.build("LaserTag")
Pkg.add("Reel")

using POMDPs
POMDPs.add("QMDP")

using LaserTag
using POMDPToolbox
using ParticleFilters
using QMDP

using Reel

rng = MersenneTwister(7)
pomdp = gen_lasertag(rng=rng, robot_position_known=true)
policy = solve(QMDPSolver(), pomdp)
filter = SIRParticleFilter(pomdp, 10000)

frames = Frames(MIME("image/png"), fps=2)

print("Simulating and generating LaserTag gif")
for step in stepthrough(pomdp, policy, filter, "a,r,sp,o,bp", rng=rng)
    push!(frames, LaserTagVis(pomdp, step...))
    print('.')
end
println(" Done.")
write("out.gif", frames)


text:
"""
    non_dominated_sort(population)

Return a vector of integers `r` containing in `r[i]` the rank for `population[i]`.
"""
function non_dominated_sort(pop)
    n = length(pop)

    dom_list = [ Int[] for i in 1:n ]
    rank = zeros(Int, n)
    dom_count = zeros(Int, n)

    for i in 1:n
        for j in i+1:n
            comparison = compare(pop[i], pop[j])
            if comparison == 1 #is_better(pop[i], pop[j])
                push!(dom_list[i], j)
                dom_count[j] += 1
            elseif comparison == 2 #is_better(pop[j], pop[i])
                push!(dom_list[j], i)
                dom_count[i] += 1
            end
        end
        if dom_count[i] == 0
            rank[i] = 1
        end
    end

    k = UInt16(2)
    while any(==(k-one(UInt16)), (rank[p] for p in 1:n)) #ugly workaround for #15276
        for p in 1:n
            if rank[p] == k-one(UInt16)
                for q in dom_list[p]
                    dom_count[q] -= one(UInt16)
                    if dom_count[q] == zero(UInt16)
                        rank[q] = k
                    end
                end
            end
        end
        k += one(UInt16)
    end

    return rank
end

"""
    get_fronts(population, computed_ranks = true)

Return each sub-front in an array. If `computed_ranks == true`, this method assumes that
`fast_non_dominated_sort!(population)` has been called before.
"""
function get_fronts(population, computed_ranks = true)
    if computed_ranks
        ranks = [s.rank for s in population]
    else
        ranks = non_dominated_sort(population)
    end

    fronts = [Int[] for i in eachindex(unique(ranks))]

    for (i, r) in enumerate(ranks)
        push!(fronts[r], i)
    end

    return fronts
end

"""
    fast_non_dominated_sort!(population)

Sort population using the fast non dominated sorting algorithm. Note that `s.rank` is
updated for each solution `s ∈ population`.
"""
function fast_non_dominated_sort!(pop)
    rank = non_dominated_sort(pop)

    for i in 1:length(pop)
        pop[i].rank = rank[i]
    end

    sort!(pop, by = s -> s.rank, alg = Base.Sort.QuickSort)
    pop
end

"""
    get_non_dominated_solutions_perm(population)

Return a vector of integers `v` such that `population[v]` are the non dominated
solutions contained in `population`.
"""
function get_non_dominated_solutions_perm(population)
    ids = Int[1]
    n = length(population)

    for i in 2:n
        j = 1
        while j <= length(ids)
            jj = ids[j]
            relation = compare(population[i], population[jj])

            if relation == 2 # j dominates i
                break
            elseif relation == 1 # i dominates j
                deleteat!(ids, j)
                continue
            end

            j += 1
        end

        if j > length(ids)
            push!(ids, i)
        end

    end

    return ids
end

"""
    get_non_dominated_solutions(population)

Return the non dominated solutions contained in `population`.
"""
function get_non_dominated_solutions(population)
    return population[get_non_dominated_solutions_perm(population)]
end


text:
module JupyterFormatter

using IJulia
using JuliaFormatter

function format_current_cell()
    current_cell = IJulia.In[IJulia.n]
    formatted = JuliaFormatter.format_text(current_cell)
    if formatted != current_cell
        IJulia.load_string(formatted, true)
    end
end

enable_autoformat() = IJulia.push_preexecute_hook(format_current_cell)
disable_autoformat() = IJulia.pop_preexecute_hook(format_current_cell)

export enable_autoformat, disable_autoformat
end


text:
using ExpmV
using Test
using LinearAlgebra
using SparseArrays
using LinearMaps

@testset "Real matrices" begin
    @testset "Positive: $positive" for positive in [true, false]
        @testset "Size: $d" for d in [10, 100]
            for i = 1:10 # Just repeat it a few times
                r = sprandn(d, d, .1)

                if positive
                    r += abs.(r)
                end

                rv = randn(d)
                rv ./= norm(rv, 2)

                rt = randn()

                @testset "Against expm" begin
                    @test expmv(rt,r,rv;check_positive=positive) ≈ exp(Matrix(rt*r))*rv
                end

                # Test the StepRangeLen version against the normal version
                @testset "Timespan $nt timesteps" for nt in [5 11 51]
                    t = range(0, stop=rt, length=nt)
                    x = expmv(t,r,rv;check_positive=positive)
                    y = hcat([expmv(ti,r,rv;check_positive=positive) for ti in t]...)
                    @test x ≈ y
                end

                # Test that it works with matrix instead of vectors
                @testset "Second dim: $d2" for d2 in 2:4
                    rv = randn(d, d2)
                    rv = rv * diagm(0 => [1.0/norm(rv[:,j],2) for j in 1:d2])

                    x = expmv(rt,r,rv;check_positive=positive)
                    @test x ≈ exp(Matrix(rt*r))*rv
                end
            end
        end
    end
end

@testset "Complex matrices" begin
    @testset "Hermitian: $herm"  for herm in [true, false]
        @testset "Size: $d" for d in [10, 100]
            for i = 1:10
                r = sprandn(d,d,.1)+1im*sprandn(d,d,.1)
                if herm
                    r = (r-r')/2
                end

                rv = randn(Complex{Float64}, d)
                rv ./= norm(rv, 2)

                rt = randn()

                x = expmv(rt,r,rv)
                @testset "Against expm" begin
                    @test x ≈ exp(Matrix(rt*r))*rv
                end

                # Test the StepRangeLen version against the normal version
                @testset "Timespan $nt timesteps" for nt in [5 11 51]
                    t = range(0, stop=rt, length=nt)
                    x = expmv(t,r,rv)
                    y = hcat([expmv(ti,r,rv) for ti in t]...)
                    @test x ≈ y
                end

                @testset "Second dim: $d2" for d2 in 2:4
                    rv = randn(Complex{Float64}, d, d2)
                    rv = rv * diagm(0 => [1.0/norm(rv[:,j],2) for j in 1:d2])

                    x = expmv(rt,r,rv)
                    @test x ≈ exp(Matrix(rt*r))*rv
                end
            end
        end
    end
end

@testset "Linear Operator" begin
    d = 100
    A = sprandn(ComplexF64,d,d,.1)
    L = LinearMap(A)

    rv = randn(ComplexF64, d)
    rv ./= norm(rv, 2)

    rt = randn()
    @test expmv(rt, L, rv; shift = false) ≈ expmv(rt, A, rv)

    t = range(0, stop=rt, length=20)
    @test expmv(t, L, rv; shift = false) ≈ expmv(t, A, rv)
end


text:
using MPI
MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
size = MPI.Comm_size(comm)

dst = mod(rank+1, size)
src = mod(rank-1, size)

N = 4

send_mesg = Array{Float64}(undef, N)
recv_mesg = Array{Float64}(undef, N)

fill!(send_mesg, Float64(rank))

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)

print("$rank: Sending   $rank -> $dst = $send_mesg\n")
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

stats = MPI.Waitall!([rreq, sreq])

print("$rank: Received $src -> $rank = $recv_mesg\n")

MPI.Barrier(comm)


text:
@testset "793.preimage-size-of-factorial-zeroes-function.jl" begin
    @test preimage_size_fzf(0) == 5
    @test preimage_size_fzf(1000000) == 5
    @test preimage_size_fzf(5) == 0
end