text:
@cache struct SymplecticEulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::SymplecticEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  SymplecticEulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct SymplecticEulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::SymplecticEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SymplecticEulerConstantCache()

@cache struct VelocityVerletCache{uType,rateType,uEltypeNoUnits} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  half::uEltypeNoUnits
end

struct VelocityVerletConstantCache{uEltypeNoUnits} <: OrdinaryDiffEqConstantCache
    half::uEltypeNoUnits
end

function alg_cache(alg::VelocityVerlet,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = zero(rate_prototype)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  half = uEltypeNoUnits(1//2)
  VelocityVerletCache(u,uprev,k,tmp,fsalfirst,half)
end

alg_cache(alg::VelocityVerlet,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = VelocityVerletConstantCache(uEltypeNoUnits(1//2))

@cache struct Symplectic2Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::VerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = VerletLeapfrogConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic2Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::VerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
      VerletLeapfrogConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

function alg_cache(alg::PseudoVerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = PseudoVerletLeapfrogConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic2Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::PseudoVerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
PseudoVerletLeapfrogConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

function alg_cache(alg::McAte2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte2ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic2Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
McAte2ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Symplectic3Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::Ruth3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = Ruth3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic3Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::Ruth3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
Ruth3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

function alg_cache(alg::McAte3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic3Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
McAte3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Symplectic4Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::McAte4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic4Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

function alg_cache(alg::CandyRoz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = CandyRoz4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic4Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::CandyRoz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Symplectic45Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::CalvoSanz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = CalvoSanz4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic45Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::CalvoSanz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = CalvoSanz4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

function alg_cache(alg::McAte42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte42ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic45Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte42ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Symplectic5Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::McAte5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic5Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Symplectic6Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::Yoshida6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = Yoshida6ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic6Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::Yoshida6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Yoshida6ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Symplectic62Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::KahanLi6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = KahanLi6ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  Symplectic62Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::KahanLi6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = KahanLi6ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct McAte8Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::McAte8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  McAte8Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct KahanLi8Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::KahanLi8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = KahanLi8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  KahanLi8Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::KahanLi8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = KahanLi8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct SofSpa10Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::SofSpa10,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = SofSpa10ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  SofSpa10Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::SofSpa10,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SofSpa10ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))


text:
#####################################
#DiracMatrix#########################
#####################################

type DiracMatrix{K<:Ket, B<:Bra, T<:Union(Number, ScalarExpr)} <: AbstractOperator
	coeffs::SparseMatrixCSC{T, Int}
	rowb::AbstractBasis{K}
	colb::AbstractBasis{B}
	function DiracMatrix(coeffs, rowb, colb)
		@assert size(coeffs)==(length(rowb), length(colb)) "Dimensions of coefficient array do not match bases"
		return new(coeffs, rowb, colb)
	end
end
DiracMatrix{K<:Ket, B<:Bra, T<:Union(Number, ScalarExpr)}(coeffs::SparseMatrixCSC{T}, rowb::AbstractBasis{K}, colb::AbstractBasis{B}) = DiracMatrix{K,B,T}(coeffs, rowb, colb)
DiracMatrix{K<:Ket, B<:Bra}(coeffs::SparseMatrixCSC, rowb::AbstractBasis{K}, colb::AbstractBasis{B}) = DiracMatrix{K,B,ScalarExpr}(convert(SparseMatrixCSC{ScalarExpr},coeffs), rowb, colb)
DiracMatrix{K<:Ket, B<:Bra}(coeffs::Array, rowb::AbstractBasis{K}, colb::AbstractBasis{B}) = DiracMatrix(qsparse(coeffs), rowb, colb)
DiracMatrix{K<:Ket}(coeffs::AbstractArray, b::AbstractBasis{K}) = dmat(coeffs, b, b')
DiracMatrix{B<:Bra}(coeffs::AbstractArray, b::AbstractBasis{B}) = dmat(coeffs, b', b)

function constructop!(coeffs, fcoeff, fstate, b) #basically never use this function except for in the below context
	@inbounds begin
		for i=1:length(b), j=1:length(b)
			coeffs[i,j] = fcoeff(b[j])*inner(b[i]',fstate(b[j]))
		end
	end
end

function DiracMatrix{K<:Ket}(fcoeff::Function, fstate::Function, b::AbstractBasis{K}, t::DataType=typeof(fcoeff(b[1])))
	coeffs = ones(t,length(b),length(b))
	constructop!(coeffs, fcoeff, fstate, b)
	return dmat(coeffs, b)
end

dmat=DiracMatrix

copy(dm::DiracMatrix) = dmat(copy(dm.coeffs), copy(dm.rowb), copy(dm.colb))

#####################################
#Getter-style Functions##############
#####################################
bsym(dm::DiracMatrix) = [bsym(dm.rowb), bsym(dm.colb)]

#####################################
#Boolean Functions###################
#####################################
isequal(a::DiracMatrix, b::DiracMatrix) = isequal(a.coeffs,b.coeffs) && a.rowb==b.rowb && a.colb==b.colb
==(a::DiracMatrix, b::DiracMatrix) = a.coeffs==b.coeffs && a.rowb==b.rowb && a.colb==b.colb
isdual(a::DiracMatrix, b::DiracMatrix) = a.coeffs'==b.coeffs && isdual(a.rowb,b.colb) && isdual(b.rowb,a.colb)
#####################################
#Array-like Functions################
#####################################
size(dm::DiracMatrix, args...) = size(dm.coeffs, args...)
for op=(:endof, :eltype, :length, :find, :findn, :findnz, :nnz, :ndims, :countnz)
	@eval ($op)(dm::DiracMatrix) = ($op)(dm.coeffs)
end

in(o::OuterProduct, dm::DiracMatrix) = in(o.ket, dm.rowb) && in(o.bra, dm.colb)
in{K<:Ket}(s::State{K}, dm::DiracMatrix) = in(s, dm.rowb)
in{B<:Bra}(s::State{B}, dm::DiracMatrix) = in(s, dm.colb)
in(n, dm::DiracMatrix) = in(n, dm.coeffs)

ctranspose(dm::DiracMatrix) = dmat(dm.coeffs', dm.colb', dm.rowb')
getindex(dm::DiracMatrix, x...) = dm.coeffs[x...]
setindex!(dm::DiracMatrix, y, x...) = setindex!(dm.coeffs,y,x...)

#####################################
#Dict-like Functions#################
#####################################
getpos{K<:Ket}(dm::DiracMatrix{K}, k::State{K}) = get(dm.rowb, k)
getpos{K<:Ket, B<:Bra}(dm::DiracMatrix{K,B}, b::State{B}) = get(dm.colb, b)
getpos{K<:Ket, B<:Bra}(dm::DiracMatrix{K,B}, k::State{K}, b::State{B}) = (getpos(dm, k), getpos(dm, b))
getpos{K<:Ket, B<:Bra}(dm::DiracMatrix{K,B}, o::OuterProduct{K,B}) = getpos(dm, o.ket, o.bra)
getpos(dm::DiracMatrix, arg) = throw(KeyError(arg))
getpos(dm::DiracMatrix, args...) = throw(KeyError(args))

get{K<:Ket}(dm::DiracMatrix, s::State{K}) = dvec(dm[getpos(dm, s), :], dm.colb)
get{B<:Bra}(dm::DiracMatrix, s::State{B}) = dvec(dm[:, getpos(dm, s)], dm.rowb)
get{K<:Ket, B<:Bra}(dm::DiracMatrix, k::State{K}, b::State{B}) = dm[getpos(dm, k), getpos(dm, b)]
get(dm::DiracMatrix, o::OuterProduct) = get(dm, o.ket, o.bra)

get{K<:Ket}(dm::DiracMatrix, s::State{K}, notfound) = in(k, dm.rowb) ? get(dm, k) : notfound
get{B<:Bra}(dm::DiracMatrix, s::State{B}, notfound) = in(b, dm.colb) ? get(dm, b) : notfound

function get{K<:Ket, B<:Bra}(dm::DiracMatrix, k::State{K}, b::State{B}, notfound)
	if in(k, dm) && in(b, dm)
		return get(dm, k, b)
	else
		return notfound
	end
end

get(dm::DiracMatrix, o::OuterProduct, notfound) = get(dm, o.ket, o.bra, notfound)

#####################################
#Show Functions######################
#####################################
summary(dm::DiracMatrix) = "$(size(dm,1))x$(size(dm,2)) $(typeof(dm))"

function showcompact(io::IO, dm::DiracMatrix)
	if length(dm.coeffs)==0
		print(io, "$(typeof(dm))[]")
	else
		tempio = IOBuffer()
		nz = hcat(findn(dm)...)
		print(tempio, [" + $(dm.coeffs[nz[i,1],nz[i,2]])$(dm.rowb[nz[i,1]])$(dm.colb[nz[i,2]])" for i=1:size(nz,1)]...)
		print(io, takebuf_string(tempio)[3:end])
	end
end

function show(io::IO, dm::DiracMatrix)
	if length(dm) > 1000000
		showsp(io, dm)
	else
		#the below is horribly inefficient for large cases,
		#hence the if statement
		println(io, summary(dm))
		table = cell(length(dm.rowb)+1, length(dm.colb)+1)
		for i = 1:length(dm.rowb)
			table[i+1,1] = dm.rowb[i]
		end
		for j = 1:length(dm.colb)
			table[1,j+1] = dm.colb[j]
		end
		table[1,1] = 0
		table[2:end, 2:end] = dm.coeffs
		temp_io = IOBuffer()
		show(temp_io, table)
		io_str = takebuf_string(temp_io)
		print(io, io_str[searchindex(io_str, "\n")+3:end])
	end
end

function showsp{K<:Ket,B<:Bra,T<:Union(Float64,Int)}(io::IO, dm::DiracMatrix{K,B,T})
	numpad = length("$(round(maximum(dm.coeffs.nzval), 4))")
	showsp_dmat(io, dm, print_el = ((io,el)->print(io, rpad("$(round(el, 4))", numpad))))
end

showsp(io::IO, dm::DiracMatrix) = showsp_dmat(io, dm)

function showsp_dmat(io::IO, dm::DiracMatrix; print_el=((io,el)->showcompact(io, el)))
 	print(io, "$(size(dm,1))x$(size(dm,2)) DiracMatrix with $(nnz(dm)) $(eltype(dm)) entries: ")
	rows = Base.tty_rows()
	S = dm.coeffs
    half_screen_rows = div(rows - 8, 2)
    pad = ndigits(max(S.m,S.n))
    k = 0
    sep = "\n\t"
    for col = 1:S.n, k = S.colptr[col] : (S.colptr[col+1]-1)
        if k < half_screen_rows || k > nnz(S)-half_screen_rows
            print(io, sep, '[', rpad(S.rowval[k], pad), ", ", lpad(col, pad), "]  =  ")
            print_el(io, S.nzval[k])
            print(io, "  $(dm.rowb[S.rowval[k]])$(dm.colb[col])")
        elseif k == half_screen_rows
            print(io, sep, '\u22ee')
        end
        k += 1
    end
end

#####################################
#Function-passing Functions##########
#####################################
find(f::Function, dm::DiracMatrix) = find(f, dm.coeffs)
findstates(f::Function, dm::DiracMatrix) = find(f, [dm.rowb[i]*dm.colb[j] for i=1:length(dm.rowb), j=1:length(dm.colb)]) #f takes OuterProduct as argument
map(f::Function, dm::DiracMatrix) = dmat(map(f, dm.coeffs), dm.rowb, dm.colb)
qeval(f::Function, dm::DiracMatrix) = map(x->qeval(f, x), dm)

#####################################
#Arithmetic Operations###############
#####################################
for op=(:.*,:.-,:.+,:./,:.^)
	@eval ($op)(a::DiracMatrix, b::DiracVector) = dmat(($op)(a.coeffs,b.coeffs), a.rowb, a.colb)
	@eval ($op)(a::DiracVector, b::DiracMatrix) = dmat(($op)(a.coeffs,b.coeffs), b.rowb, b.colb)
	@eval ($op)(a::DiracMatrix, b::DiracMatrix) = dmat(($op)(a.coeffs,b.coeffs), a.rowb, a.colb)
	@eval ($op)(n::DiracCoeff, d::DiracMatrix) = dmat(($op)(n,d.coeffs), d.rowb, d.colb)
	@eval ($op)(d::DiracMatrix, n::DiracCoeff) = dmat(($op)(d.coeffs,n), d.rowb, d.colb)
end

/(dm::DiracMatrix, d::DiracCoeff) = dmat(dm.coeffs/d, dm.rowb, dm.colb)
*(dm::DiracMatrix, d::DiracCoeff) = dmat(dm.coeffs*d, dm.rowb, dm.colb)
*(d::DiracCoeff, dm::DiracMatrix) = dmat(d*dm.coeffs, dm.rowb, dm.colb)

kron(c::DiracCoeff, dm::DiracMatrix) = c*dm
kron(dm::DiracMatrix, c::DiracCoeff) = dm*c

function multbystate!{K<:Ket}(dm::DiracMatrix, s::State{K}, dv::AbstractArray)
	for i in findn(dm)[1]
		coeff = 0
		for j in findn(dm)[2]
			coeff+=dm[i,j]*(dm.colb[j]*s)
		end
		dv[i]=dv[i]+coeff
	end
	return dvec(dv, dm.rowb)
end

function multbystate!{B<:Bra}(dm::DiracMatrix, s::State{B}, dv::AbstractArray)
	for j in findn(dm)[2]
		coeff = 0
		for i in findn(dm)[1]
			coeff+=dm[i,j]*(s*dm.rowb[i])
		end
		dv[j]=dv[j]+coeff
	end
	return dvec(dv, dm.colb)
end

consmult{K<:Ket}(dm::DiracMatrix, s::State{K}) = multbystate!(dm, s, samebasis(dm.colb,s) ? zeros(length(dm.rowb)) : convert(Array{Any}, zeros(length(dm.rowb))))
consmult{B<:Bra}(dm::DiracMatrix, s::State{B}) = multbystate!(dm, s, samebasis(dm.rowb,s) ? zeros(1,length(dm.colb)) : convert(Array{Any}, zeros(1,length(dm.rowb))))

#currently the below consmults are slow;
#should be replaced by a something
#more in line with the above single
#state versions, which preallocates
#memory and then fills the array using
#a separate function
consmult{K<:Ket}(dm::DiracMatrix, d::DiracVector{K}) = reduce(+,[dm.rowb[i]*reduce(+,[dm[i,j]*(dm.colb[j]*d) for j=1:length(dm.colb)]) for i=1:length(dm.rowb)])
consmult{B<:Bra}(dm::DiracMatrix, d::DiracVector{B}) = reduce(+,[dm.colb[j]*reduce(+,[dm[i,j]*(d*op.rowb[i]) for j=1:length(dm.rowb)]) for i=1:length(dm.colb)])
consmult(a::DiracMatrix, b::DiracMatrix) = reduce(+,[(a[i,j]*b[m,n]*(a.colb[j]*b.rowb[m]))*(a.rowb[i]*b.colb[n]) for n=1:size(b,2), m=1:size(b,1), j=1:size(a,2), i=1:size(a,1)])

*{K<:Ket}(dm::DiracMatrix, s::State{K}) = in(s', dm.colb) ? get(dm, s') : consmult(dm,s)
*{B<:Bra}(s::State{B}, dm::DiracMatrix) = in(s', dm.rowb) ? get(dm, s') : consmult(dm,s)
*{K<:Ket}(dm::DiracMatrix, dv::DiracVector{K}) = isdual(dm.colb, dv.basis) ? dvec(dm.coeffs*dv.coeffs, dm.rowb) : consmult(dm,dv)
*{B<:Bra}(dv::DiracVector{B},dm::DiracMatrix) = isdual(dm.rowb, dv.basis) ? dvec(dm.coeffs*dv.coeffs, dm.colb) : consmult(dm,dv)
*(a::DiracMatrix, b::DiracMatrix) = isdual(a.colb, b.rowb) ? dmat(a.coeffs*b.coeffs, a.rowb, b.colb) : consmult(a,b)
*(dm::DiracMatrix, o::OuterProduct) = (dm*o.ket)*o.bra
*(o::OuterProduct, dm::DiracMatrix) = o.ket*(o.bra*dm)

kron(a::DiracMatrix, b::DiracMatrix) = dmat(kron(a.coeffs, b.coeffs), tensor(a.rowb, b.rowb), tensor(a.colb, b.colb))
kron{K<:Ket}(dm::DiracMatrix, dv::DiracVector{K}) = dmat(kron(dm.coeffs, dv.coeffs), tensor(dm.rowb, dv.basis), dm.colb)
kron{K<:Ket}(dv::DiracVector{K}, dm::DiracMatrix) = dmat(kron(dv.coeffs, dm.coeffs), tensor(dv.basis, dm.rowb), dm.colb)
kron{B<:Bra}(dm::DiracMatrix, dv::DiracVector{B}) = dmat(kron(dm.coeffs, dv.coeffs), dm.rowb, tensor(dm.colb, dv.basis))
kron{B<:Bra}(dv::DiracVector{B}, dm::DiracMatrix) = dmat(kron(dv.coeffs, dm.coeffs), dm.rowb, tensor(dv.basis, dm.colb))
kron{K<:Ket}(s::State{K}, dm::DiracMatrix) = dmat(dm.coeffs, tensor(s, dm.rowb), dm.colb)
kron{K<:Ket}(dm::DiracMatrix, s::State{K}) = dmat(dm.coeffs, tensor(s, dm.rowb), dm.colb)
kron{B<:Bra}(s::State{B}, dm::DiracMatrix) = dmat(dm.coeffs, dm.rowb, tensor(s, dm.colb))
kron{B<:Bra}(dm::DiracMatrix, s::State{B}) = dmat(dm.coeffs, dm.rowb, tensor(s, dm.colb))
kron(o::OuterProduct, dm::DiracMatrix) = dmat(dm.coeffs, tensor(o.ket, dm.rowb), tensor(o.bra, dm.colb))
kron(dm::DiracMatrix, o::OuterProduct) = dmat(dm.coeffs, tensor(dm.rowb, o.ket), tensor(dm.colb, o.bra))

function add_out!(dm, o)
	dm[getpos(dm, o)...] = 1+get(dm, o)
	return dm
end

add_out(dm, o) = add_out!(copy(dm), o)
add_out{K,B,T<:InnerProduct}(dm::DiracMatrix{K,B,T}, o) = add_out!(1*dm, o) #converts res to DiracMatrix{K,B,ScalarExpr}

function add_bra(dm::DiracMatrix, o::OuterProduct)
	newcol = zeros(size(dm,1))
	newcol[get(dm.rowb, o.ket)] = 1
	return dmat(hcat(dm.coeffs,newcol), dm.rowb, bcat(dm.colb, o.bra))
end

function add_bra(o::OuterProduct, dm::DiracMatrix)
	newcol = zeros(size(dm,1))
	newcol[get(dm.rowb, o.ket)] = 1
	return dmat(hcat(newcol, dm.coeffs), dm.rowb, bcat(o.bra, dm.colb))
end

function add_ket(dm::DiracMatrix, o::OuterProduct)
	newrow = zeros(1,size(dm,2))
	newrow[get(dm.colb, o.bra)] = 1
	return dmat(vcat(dm.coeffs,newrow), bcat(dm.rowb, o.ket), dm.colb)
end

function add_ket(o::OuterProduct, dm::DiracMatrix)
	newrow = zeros(1,size(dm,2))
	newrow[get(dm.colb, o.bra)] = 1
	return dmat(vcat(newrow,dm.coeffs), bcat(o.ket,dm.rowb), dm.colb)
end

function add_both!(dm::DiracMatrix, o::OuterProduct)
	dm[getpos(dm, o)...] = 1+get(dm, o)
	return dm
end
function add_both!(dm::DiracMatrix, o::OuterProduct,rb,cb,res)
	res[1:size(dm,1), 1:size(dm,2)] = dm.coeffs
	return add_both!(dmat(res, rb, cb), o)
end

function add_both!(o::OuterProduct,dm::DiracMatrix,rb,cb,res)
	res[2:size(dm,1)+1, 2:size(dm,2)+1] = dm.coeffs
	return add_both!(dmat(res, rb, cb), o)
end

add_both(dm::DiracMatrix, o::OuterProduct, rb, cb) = add_both!(dm, o, rb, cb, zeros(length(rb), length(cb)))
add_both{K,B,T}(dm::DiracMatrix{K,B,T}, o::OuterProduct) = add_both(dm, o, bcat(dm.rowb, o.ket), bcat(dm.colb, o.bra))
add_both{K,B,T<:InnerProduct}(dm::DiracMatrix{K,B,T}, o::OuterProduct) = add_both(1*dm, o)

add_both(o::OuterProduct, dm::DiracMatrix, rb, cb) = add_both!(o, dm, rb, cb, zeros(length(rb), length(cb)))
add_both{K,B,T}(o::OuterProduct, dm::DiracMatrix{K,B,T}) = add_both(o, dm, bcat(o.ket, dm.rowb), bcat(o.bra, dm.colb))
add_both{K,B,T<:InnerProduct}(o::OuterProduct, dm::DiracMatrix{K,B,T}) = add_both(o,1*dm)

function +{K<:Ket, B<:Bra}(op::DiracMatrix{K,B}, o::OuterProduct{K,B})
	if in(o.ket, op.rowb)
		if in(o.bra, op.colb)
			return add_out(op, o)
		else
			return add_bra(op, o)
		end
	elseif in(o.bra, op.colb)
		return add_ket(op, o)
	else
		return add_both(op,o)
	end
end

function +{K<:Ket, B<:Bra}(o::OuterProduct{K,B}, op::DiracMatrix{K,B})
	if in(o.ket, op.rowb)
		if in(o.bra, op.colb)
			return add_out(op,o)
		else
			return add_bra(o,op)
		end
	elseif in(o.bra, op.colb)
		return add_ket(o,op)
	else
		return add_both(o,op)
	end
end

function add_dmat!(a::DiracMatrix, b::DiracMatrix)
	for i in findn(b)[1]
		for j in findn(b)[2]
			a=a+(b.rowb[i]*b.colb[j])
			a[getpos(a, b.rowb[i], b.colb[j])...] = get(a, b.rowb[i], b.colb[j])+b[i,j]-1
		end
	end
	return a
end

function +{K<:Ket, B<:Bra}(a::DiracMatrix{K,B}, b::DiracMatrix{K,B})
	if a.rowb==b.rowb && a.colb==b.colb
		return dmat(a.coeffs+b.coeffs, a.rowb, a.colb)
	else
		return add_dmat!(copy(a), b)
	end
end

-(op::DiracMatrix) = -1*op
-(a::DiracMatrix, b::DiracMatrix) = a+(-b)
-(a::DiracMatrix, b::OuterProduct) = a+(-b)
-(a::OuterProduct, b::DiracMatrix) = a+(-b)

exp(op::DiracMatrix) = dmat(exp(op.coeffs), op.rowb, op.colb)
^(op::DiracMatrix, i::Integer) = dmat(^(op.coeffs, i), op.rowb, op.colb)
^(op::DiracMatrix, n::Number) = dmat(^(op.coeffs, n), op.rowb, op.colb)
trace(op::DiracMatrix) = trace(op.coeffs)

commutator(a::DiracMatrix, b::DiracMatrix) = (a*b) - (b*a)

function ptrace(op::DiracMatrix, ind::Int)
	@assert isdual(op.colb, op.rowb) "BasisMismatch"
	trrow = tensor(vcat(separate(op.rowb)[1:ind-1], separate(op.rowb)[ind+1:end])...)
	trcol = trrow'
	len = length(trcol)
	#the line below should be optimized by pre-allocating the array and
	#writing an in-place function to fill it
	coeffs = [reduce(+,[op[(((i-1)*len)+k), (((i-1)*len)+j)] for i=1:length(separate(op.rowb)[ind])]) for j=1:len, k=1:len]
	# return dmat(vcat([hcat(coeffs[i, :]...) for i=1:size(coeffs, 1)]...), trrow, trcol) #use vcat()/hcat() trick to convert to most primitive common type
	return dmat(coeffs, trrow, trcol)
end

actop{K<:Ket}(op::AbstractOperator, s::Tensor{K}, i::Integer) = kron(vcat(s[1:i-1], op*s[i], s[i+1:end])...)
actop{B<:Bra}(op::AbstractOperator, s::Tensor{B}, i::Integer) = kron(vcat(s[1:i-1], s[i]*op, s[i+1:end])...)

ext: jl


text:
export construct_univariate_relaxation!, construct_bilinear_relaxation!

"""
    construct_univariate_relaxation!(m,f,x,y,x_partition;f_dash=x->ForwardDiff.derivative(f,x),error_tolerance=NaN64,length_tolerance=1e-6,derivative_tolerance=1e-6,num_additional_partitions=0)

Add MILP relaxation of `y=f(x)` to given JuMP model and return an object with
new variables and constraints.

# Mandatory Arguments
- `m::Jump.Model`: model to which relaxation is to be added.
- `f::Function`: function or oracle for which a polyhedral relaxation is
    required, usually non-linear.
- `x::Jump.VariableRef`: JuMP variable for domain of `f`.
- `y::JuMP.VariableRef`: JuMP variable for evaluation of `f`.
- `x_partition::Vector{<:Real}`: partition of the domain of `f`.
- `milp::Bool`: build MILP relaxation if true, LP relaxation otherwise. Note
    that the MILP relaxation uses the incremental formulation presented in the
    paper, but the LP relaxation uses a lambda form that builds a formulation
    as the convex combination of triangle vertices that are at the intersection
    of tangents to the curve.

# Optional Arguments
- `f_dash::Function`: function or oracle for derivative of `f`, defaults to
    the `derivative` function from the `ForwardDiff` package.
- ` error_tolerance::Float64`: Maximum allowed vertical distance between over
    and under estimators of `f`, defaults to NaN64.
- `length_tolerance::Float64`: maximum length of a sub-interval in a partition,
    defaults to ``1 \\times 10^{-6}``.
- `derivative_tolerance::Float64`: minimum absolute difference between
    derivaties at successive elements of a partition for them to be considered
    different, defaults to ``1 \\times 10^{-6}``. If the difference of a partition sub-interval
    is smaller than this value, that sub-interval will be refined.
- `num_additional_partitions::Int64`: budget on number of sub-intervals in
    partition, defaults to 0. Note that if the number of partitions is `n` and
    the number of additional partitions is `m`, then the function will return a
    relaxation with at most `n+m` partitions.
- `variable_pre_base_name::AbstractString`: base_name that needs to be added to the auxiliary
    variables for meaningful LP files

Assume that:
- `f` is a bounded function of 1 variable.
- `x_partition` is a partition of the domain of `f` such that `f` is either
    convex or concave each sub-interval of the partition.
- `f_dash` is not equal at two consecutive elements of `x_partition`.

This function builds an incremental formulation, which is the formulation
presented in the paper.
"""
function construct_univariate_relaxation!(
    m::JuMP.Model,
    f::Function,
    x::JuMP.VariableRef,
    y::JuMP.VariableRef,
    x_partition::Vector{<:Real},
    milp::Bool;
    f_dash::Function = x -> ForwardDiff.derivative(f, x),
    error_tolerance::Float64 = NaN64,
    length_tolerance::Float64 = EPS,
    derivative_tolerance::Float64 = EPS,
    num_additional_partitions::Int64 = 0,
    variable_pre_base_name::AbstractString = ""
)::FormulationInfo
    univariate_function_data = UnivariateFunctionData(
        f,
        f_dash,
        x_partition,
        error_tolerance,
        length_tolerance,
        derivative_tolerance,
        num_additional_partitions,
        length(x_partition),
    )
    _validate(univariate_function_data)
    _validate(x, x_partition)
    _refine_partition!(univariate_function_data)
    func = milp ? _build_univariate_milp_relaxation! : _build_univariate_lp_relaxation!
    return func(m, x, y, univariate_function_data, variable_pre_base_name)
end

"""
    construct_bilinear_relaxation!(m,x,y,z,x_partition,y_partition)

Add polyhedral relaxation of `z = xy` to given JuMP model and return an object with
new variables and constraints.

# Mandatory Arguments
- `m::Jump.Model`: model to which relaxation is to be added.
- `x::Jump.VariableRef`: JuMP variable `x`.
- `y::JuMP.VariableRef`: JuMP variable `y`.
- `z::JuMP.VariableRef`: JuMP variable `z`.
- `x_partition::Vector{<:Real}`: partition of the domain of `x`.
- `y_partition::Vector{<:Real}`: partition of the domain of `y`.

# Optional Arguments
- `variable_pre_base_name::AbstractString`: base_name that needs to be added to the auxiliary
    variables for meaningful LP files

This function builds an incremental formulation, and currently supports more than
one partition only on one of the variables `x` or `y` and not on both. It will
throw an error when more than one partitions are provided on both variables.
When exactly one partition is input for both variables, it populates the model
with the McCormick relaxation. The incremental formulation is similar to the triangle
chain relaxation in the manuscript with the triangles replaced with tetrahedrons.
Adjacent tetrahedrons share an edge (instead of a point in the triangles case).
"""
function construct_bilinear_relaxation!(
    m::JuMP.Model,
    x::JuMP.VariableRef,
    y::JuMP.VariableRef,
    z::JuMP.VariableRef,
    x_partition::Vector{<:Real},
    y_partition::Vector{<:Real};
    variable_pre_base_name::AbstractString = ""
)::FormulationInfo
    _validate(x, y, x_partition, y_partition)
    if length(x_partition) == 2 && length(y_partition) == 2
        return _build_mccormick_relaxation!(m, x, y, z)
    end
    return _build_bilinear_milp_relaxation!(m, x, y, z, x_partition, y_partition, variable_pre_base_name)
end


text:
# high level functions used in the project
using PyPlot
using DelimitedFiles

function run_fit_procedure(description::String, data::Array{Int,2}, algo_fit_spec,
                            nr_repetitions::Int, mcmc_data_file::String,
                            top_n_fraction::Float64, rl_module_path::String, opts...)
    start_time = now()
    println("start: $start_time")
    println(description)
    data_string = RL_Fit.to_short_string(RL_Fit.get_data_stats(data))
    println(data_string)
        flush(stdout)
    println("run_fit_procedure")

    mcmc_samples = RL_Fit.run_parallel_Metropolis(data, algo_fit_spec, rl_module_path, opts...; nr_of_repetitions=nr_repetitions)
    nr_map_samples = Int(ceil(size(mcmc_samples)[1] * top_n_fraction))
    println("sampling data size: $(size(mcmc_samples)), top_N: $nr_map_samples")
    top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples = RL_Fit.get_MAP(mcmc_samples, nr_map_samples)
    println("result timestamp: $(now())")
    println("top_MAP: $(round.(top_MAP,digits=4)) (Max Likelihood estimate (top 1 sample). LL includes Prior LL)")
    println("top_n_avgMAP: $(round.(top_n_avgMAP,digits=4)) . (Avg over top $nr_map_samples samples. LL includes Prior LL)")
    println("top_n_min: $(round.(top_n_min,digits=4))")
    println("top_n_max: $(round.(top_n_max,digits=4))")
    # evaluate using top_n_avgMAP parameter. get the LL without the prior (penalty)
    LL_noPrior_avgMAP, all_signals_avg = algo_fit_spec.rl_impl(data, top_n_avgMAP[1:end-1], opts...)
    println("LL top_n_avgMAP=$(round(LL_noPrior_avgMAP, digits=4))")
    # evaluate using top_MAP parameter. get the LL without the prior (penalty)
    LL_noPrior_MAP, all_signals = algo_fit_spec.rl_impl(data, top_MAP[1:end-1], opts...)
    println("LL top_MAP=$(round(LL_noPrior_MAP, digits=4))")

    flush(stdout)
    end_time = now()
    println("end of fitting: $end_time, duration: $(end_time-start_time)")
    top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples, mcmc_samples, LL_noPrior_MAP, LL_noPrior_avgMAP, all_signals
end

function _write_subj_data_file(subj_data_path::String , data_stats::Data_stats ,subj_signals::Step_by_step_signals, pop_signals::Step_by_step_signals, subj_fit_mcmc_figure)
    # writes the data to the file system.
    # the directory subj_data_path is created if it does not yet exist.
    if !ispath(subj_data_path)
        mkpath(subj_data_path)
    end

    #learning signals when per-subject (posterior) params are used
    path = "$(subj_data_path)params_subj.csv"
    writedlm(path, subj_signals.param_vector)

    # always write the action_selection_prob
    writedlm("$(subj_data_path)action_selection_prob_subj.csv", subj_signals.action_selection_prob)

    if subj_fit_mcmc_figure != nothing
    # keep track of the MCMC samples for later visualizations
    # writedlm("$(subj_data_path)/mcmc_samples_subj.csv", subj_fit_mcmc_samples)
    # plot
        subj_fit_mcmc_figure.savefig("$(subj_data_path)MCMC.png")   # save the figure to file
    end

    if length(subj_signals.RPE)>0
      path = "$(subj_data_path)RPE_subj.csv"
        d = hcat(subj_signals.RPE, subj_signals.action_selection_prob)
        writedlm(path, d)
    end
    if length(subj_signals.SPE)>0
        path = "$(subj_data_path)SPE_subj.csv"
        d = hcat(subj_signals.SPE, subj_signals.action_selection_prob)
        writedlm(path, d)
    end

    # learning signals when population (prior) params are used:
    path = "$(subj_data_path)params_pop.csv"
    writedlm(path, pop_signals.param_vector)

    # always write the action_selection_prob
    writedlm("$(subj_data_path)action_selection_prob_pop.csv", pop_signals.action_selection_prob)

    if length(pop_signals.RPE)>0
        path = "$(subj_data_path)RPE_pop.csv"
          d = hcat(pop_signals.RPE, pop_signals.action_selection_prob)
          writedlm(path, d)

    end
    if length(pop_signals.SPE)>0
        path = "$(subj_data_path)SPE_pop.csv"
        writedlm(path, hcat(pop_signals.SPE, pop_signals.action_selection_prob))
    end

    if haskey(pop_signals.results_dict, "RPE_QV")
        path = "$(subj_data_path)RPE_QV_pop.csv"
        writedlm(path, pop_signals.results_dict["RPE_QV"])
    end
    if haskey(subj_signals.results_dict, "RPE_QV")
        path = "$(subj_data_path)RPE_QV_subj.csv"
        writedlm(path, subj_signals.results_dict["RPE_QV"])
    end

    if haskey(pop_signals.results_dict, "trialByTrial_Vs")
        path = "$(subj_data_path)StateValue_pop.csv"
        writedlm(path, pop_signals.results_dict["trialByTrial_Vs"])
    end
    if haskey(subj_signals.results_dict, "trialByTrial_Vs")
        path = "$(subj_data_path)StateValue_subj.csv"
        writedlm(path, subj_signals.results_dict["trialByTrial_Vs"])
    end

    ###############################
    # Store V-values (if there are any in the results_dict)
    if haskey(pop_signals.results_dict, "StateValue_V1_record")
        path = "$(subj_data_path)StateValue_V1_pop.csv"
        writedlm(path, pop_signals.results_dict["StateValue_V1_record"])
    end
    if haskey(subj_signals.results_dict, "StateValue_V1_record")
        path = "$(subj_data_path)StateValue_V1_subj.csv"
        writedlm(path, subj_signals.results_dict["StateValue_V1_record"])
    end
    if haskey(pop_signals.results_dict, "StateValue_V2_record")
        path = "$(subj_data_path)StateValue_V2_pop.csv"
        writedlm(path, pop_signals.results_dict["StateValue_V2_record"])
    end
    if haskey(subj_signals.results_dict, "StateValue_V2_record")
        path = "$(subj_data_path)StateValue_V2_subj.csv"
        writedlm(path, subj_signals.results_dict["StateValue_V2_record"])
    end
    if haskey(pop_signals.results_dict, "StateValue_V3_record")
        path = "$(subj_data_path)StateValue_V3_pop.csv"
        writedlm(path, pop_signals.results_dict["StateValue_V3_record"])
    end
    if haskey(subj_signals.results_dict, "StateValue_V3_record")
        path = "$(subj_data_path)StateValue_V3_subj.csv"
        writedlm(path, subj_signals.results_dict["StateValue_V3_record"])
    end
    ##############################

    if haskey(pop_signals.results_dict, "Surprise")
        path = "$(subj_data_path)Surprise_pop.csv"
        writedlm(path, pop_signals.results_dict["Surprise"])
    end
    if haskey(subj_signals.results_dict, "Surprise")
        path = "$(subj_data_path)Surprise_subj.csv"
        writedlm(path, subj_signals.results_dict["Surprise"])
    end
    if haskey(pop_signals.results_dict, "gamma_surprise")
        path = "$(subj_data_path)gamma_surprise_pop.csv"
        writedlm(path, pop_signals.results_dict["gamma_surprise"])
    end
    if haskey(subj_signals.results_dict, "gamma_surprise")
        path = "$(subj_data_path)gamma_surprise_subj.csv"
        writedlm(path, subj_signals.results_dict["gamma_surprise"])
    end
    ##############################
    if haskey(pop_signals.results_dict, "DeltaPolicy")
        path = "$(subj_data_path)DeltaPolicy_pop.csv"
        writedlm(path, pop_signals.results_dict["DeltaPolicy"])
    end
    if haskey(subj_signals.results_dict, "DeltaPolicy")
        path = "$(subj_data_path)DeltaPolicy_subj.csv"
        writedlm(path, subj_signals.results_dict["DeltaPolicy"])
    end
    #
    if haskey(pop_signals.results_dict, "PolicyParametersImbalance")
        path = "$(subj_data_path)PolicyParametersImbalance_pop.csv"
        writedlm(path, pop_signals.results_dict["PolicyParametersImbalance"])
    end
    if haskey(subj_signals.results_dict, "PolicyParametersImbalance")
        path = "$(subj_data_path)PolicyParametersImbalance_subj.csv"
        writedlm(path, subj_signals.results_dict["PolicyParametersImbalance"])
    end
    #
    if haskey(pop_signals.results_dict, "PolicyImbalance")
        path = "$(subj_data_path)PolicyImbalance_pop.csv"
        writedlm(path, pop_signals.results_dict["PolicyImbalance"])
    end
    if haskey(subj_signals.results_dict, "PolicyImbalance")
        path = "$(subj_data_path)PolicyImbalance_subj.csv"
        writedlm(path, subj_signals.results_dict["PolicyImbalance"])
    end
    #
    if haskey(pop_signals.results_dict, "PolicyParametersImbalanceAllstates")
        path = "$(subj_data_path)PolicyParametersImbalanceAllstates_pop.csv"
        writedlm(path, pop_signals.results_dict["PolicyParametersImbalanceAllstates"])
    end
    if haskey(subj_signals.results_dict, "PolicyParametersImbalanceAllstates")
        path = "$(subj_data_path)PolicyParametersImbalanceAllstates_subj.csv"
        writedlm(path, subj_signals.results_dict["PolicyParametersImbalanceAllstates"])
    end
    #
    if haskey(pop_signals.results_dict, "trial_by_trial_p_sa_max")
        path = "$(subj_data_path)action_selection_prob_MAX_pop.csv"
        writedlm(path, pop_signals.results_dict["trial_by_trial_p_sa_max"])
    end
    if haskey(subj_signals.results_dict, "trial_by_trial_p_sa_max")
        path = "$(subj_data_path)action_selection_prob_MAX_subj.csv"
        writedlm(path, subj_signals.results_dict["trial_by_trial_p_sa_max"])
    end

    # LL of subject
    path = "$(subj_data_path)LL_pop.csv"
    writedlm(path, pop_signals.LL)

    path = "$(subj_data_path)LL_subj.csv"
    writedlm(path, subj_signals.LL)

    # BIC of subject
    path = "$(subj_data_path)BIC_pop.csv"
    BIC_pop = get_BIC(pop_signals.LL, length(pop_signals.param_vector), data_stats.nr_actions)
    writedlm(path, BIC_pop)

    path = "$(subj_data_path)BIC_subj.csv"
    BIC_subj = get_BIC(subj_signals.LL, length(subj_signals.param_vector), data_stats.nr_actions)
    writedlm(path, BIC_subj)

    #
    if haskey(pop_signals.results_dict, "PolicyParametersDiffChosenVsUnchosen")
        path = "$(subj_data_path)PolicyParametersDiffChosenVsUnchosen_pop.csv"
        writedlm(path, pop_signals.results_dict["PolicyParametersDiffChosenVsUnchosen"])
    end
    if haskey(subj_signals.results_dict, "PolicyParametersDiffChosenVsUnchosen")
        path = "$(subj_data_path)PolicyParametersDiffChosenVsUnchosen_subj.csv"
        writedlm(path, subj_signals.results_dict["PolicyParametersDiffChosenVsUnchosen"])
    end
    ##############################

    # textfile describing unstructured data for info only
    path = "$(subj_data_path)description.txt"
    f = open(path, "w")
    write(f,"DateTime of fitting: $(now())\n")
    write(f,"DATA_STATS\n")
    write(f,to_string(data_stats))
    write(f,"\nSUBJECT_SIGNALS\n")
    write(f, to_string(subj_signals))
    write(f,"\nPOPULATION_SIGNALS\n")
    write(f, to_string(pop_signals))
    close(f)
end

function run_single_subj_fit_write_csv(
    data_path::String,
    base_path::String,
    id_subject::Int64,
    fit_max_episodes_per_subj::Int64,
    learner_function_V::Function,
    population_prior,
    population_map_params,
    burn_in::Int64,
    nr_mcmc_samples::Int64,
    inter_collect_interval::Int64,
    nr_of_repetitions::Int64,
    top_n_fraction::Float64,
    descritpion::String,
    rl_module_path::String;
    write_mcmc_picture_file = true)
    # a higher level script that writes all the information we currently
    # want to the file system.

    data_single_subj, data_stats = load_SARSPEI(data_path; subject_id = id_subject,
                                            max_episodes_per_subj=fit_max_episodes_per_subj)

    mSpec = Metropolis_Spec(learner_function_V, population_prior, burn_in,nr_mcmc_samples, inter_collect_interval)
    mcmc_samples = run_parallel_Metropolis(data_single_subj, mSpec, rl_module_path; nr_of_repetitions=nr_of_repetitions)

    nr_map_samples = Int(ceil(size(mcmc_samples)[1] * top_n_fraction))
    top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples = get_MAP(mcmc_samples, nr_map_samples)

    map_val_subj = deepcopy(top_MAP)

    println("subj: $id_subject, population prior, map_val_subj: $(round.(map_val_subj,digits=4))")
    flush(stdout)
    # for evaluation (unlike fitting) load ALL Episodes
    data_single_subj, data_stats = load_SARSPEI(data_path; subject_id = id_subject)
    LL, subj_signals = learner_function_V(data_single_subj, map_val_subj)
    subj_signals = subj_signals[id_subject]

    LL_pop, pop_signals = learner_function_V(data_single_subj, population_map_params)
    pop_signals = pop_signals[id_subject]

    subj_fit_mcmc_figure = nothing
    if write_mcmc_picture_file
        fig_title = "Subj:$id_subject, $(RL_Fit.to_short_string(data_stats)), $learner_function_V, LL_sub: $(round(LL,digits=2)), LL_pop: $(round(LL_pop,digits=2))"
        subj_fit_mcmc_figure = RL_Fit.plot_result(mcmc_samples, population_prior, fig_title)#, top_n_fraction = 1.)
    end

    subj_data_path = "$(base_path)subj_$(id_subject)/"
    _write_subj_data_file(subj_data_path, data_stats,subj_signals, pop_signals, subj_fit_mcmc_figure)
    plt.close("all")
    subj_signals, pop_signals, mcmc_samples, data_stats
end

function get_BIC(LL, nr_params, nr_actions)
    BIC = -2. * LL + nr_params * log(nr_actions)
end

function get_BIClogevidence(LL, nr_params, nr_actions)
    BIClogevidence = LL - nr_params * log(nr_actions) / 2.
    # or
    # BIClogevidence = - get_BIC(LL, nr_params, nr_actions) / 2.
end


text:
# QR-decomposition using the MSG algorithm

R = zeros(n,n)

for k=1:n
    R[k,k] = norm(A[:,k])
    A[:,k] /= R[k,k]
    R[k:k,k+1:n] = A[:,k]’*A[:,k+1:n]
    A[:,k+1:n]  -= A[:,k]*R[k:k,k+1:n]
end

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Numerical Linear Algebra/QRDecomposition.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
function soyster_load_parameters_chain!(
    ata_model::SoysterMaximinModel;
    items_file = "items.jld2",
    items::Vector{Psychometrics.Item} = Psychometrics.Item[],
    kwargs...,
)

    try
        T = ata_model.settings.T
        n_items = ata_model.settings.n_items
        irt_parameters = ata_model.settings.irt.parameters
        irt_model = ata_model.settings.irt.model
        irt_D = ata_model.settings.irt.D
        irt_parametrization = ata_model.settings.irt.parametrization
        R = ata_model.obj.R
        n_items = ata_model.settings.n_items
        if !isfile(items_file) || (size(items, 1) > 0)
            error!(
                ata_model,
                string(
                    "",
                    items_file,
                    " does not exist.\nProvide a `Psychometrics.Item` vector through the `items` argument, or the name of a file which can be loaded by `FileIO` that contains that vector through the `items_file` argument.",
                ),
            )
        elseif isfile(items_file)
            items = FileIO.load(items_file)
            items = items[collect(keys(items))[1]]
        end
        IIF = Vector{Array{Float64,2}}(undef, T)
        # ICF = Vector{Array{Float64,2}}(undef, T)
        K = zeros(Int, T)
        for t = 1:T
            K[t] = size(ata_model.obj.cores[t].points, 1)
            IIF[t] = zeros(K[t], n_items)
            # ICF[t] = zeros(K[t], n_items)
        end
        for i = 1:n_items
            item = items[i]
            if size(item.parameters.chain, 1) != R
                error!(ata_model, string("Item 1 does not have R=", R, " chains."))
            end
            if irt_model == "1PL"
                if !(item.parameters isa Psychometrics.Parameters1PL)
                    error!(ata_model, string("Item 1 is not of type ", irt_model, "."))
                end
                df = DataFrames.DataFrame(b = [item.parameters.chain[r][1]])
            elseif irt_model == "2PL"
                if !(item.parameters isa Psychometrics.Parameters2PL)
                    error!(ata_model, string("Item 1 is not of type ", irt_model, "."))
                end
                df = DataFrames.DataFrame(
                    a = [item.parameters.chain[r][1] for r = 1:R],
                    b = [item.parameters.chain[r][2] for r = 1:R],
                )
            elseif irt_model == "3PL"
                if !(item.parameters isa Psychometrics.Parameters3PL)
                    error!(ata_model, string("Item 1 is not of type ", irt_model, "."))
                end
                df = DataFrames.DataFrame(
                    a = [item.parameters.chain[r][1] for r = 1:R],
                    b = [item.parameters.chain[r][2] for r = 1:R],
                    c = [item.parameters.chain[r][3] for r = 1:R],
                )
            end
            for t = 1:T
                IIF_i = zeros(R, K[t])
                IIF_i_t_min = Inf .* ones(K[t])
                #check if chain has length R
                for k = 1:K[t]
                    IIF_i[:, k] = item_info(
                        df,
                        ata_model.obj.cores[t].points[k];
                        model = irt_model,
                        parametrization = irt_parametrization,
                        D = irt_D,
                    )
                    # ICF_i[:, k] = item_char(
                    #     df,
                    #     ata_model.obj.cores[t].points[k];
                    #     model = irt_model,
                    #     parametrization = irt_parametrization,
                    #     D = irt_D,
                    # )[1][
                    #     :,
                    #     :,
                    #     1,
                    # ]
                    IIF[t][k, i] = max(
                        0.01,
                        StatsBase.mean(IIF_i[:, k]) - 3 * StatsBase.std(IIF_i[:, k]),
                    )
                    # ICF[t][k, i] = StatsBase.mean(IIF_i[:, k]) - StatsBase.std(IIF_i[:, k])
                end
            end
        end
        for t = 1:T
            ata_model.obj.cores[t].IIF = IIF[t]
        end
    catch e
        error!(ata_model, string(sprint(showerror, e)))

    end
    return nothing
end


text:
using JSON
using Base.Test
using Compat
import DataStructures

include(joinpath(dirname(@__FILE__),"json_samples.jl"))

@test JSON.parse("{\"x\": 3}"; ordered = true) == DataStructures.OrderedDict{String,Any}([("x",3)])

# Test definitions -------
validate_c(c) = begin
                    j = JSON.parse(c);
                    @test j != nothing
                    @test typeof(j["widget"]["image"]["hOffset"]) == Int
                    @test j["widget"]["image"]["hOffset"] == 250
                    @test typeof(j["widget"]["text"]["size"]) == Float64
                    @test j["widget"]["text"]["size"] == 36.5
                end

validate_e(e) = begin
                    j=JSON.parse(e)
                    @test j != nothing
                    @test typeof(j) == Dict{AbstractString, Any}
                    @test length(j) == 1
                    @test typeof(j["menu"]) == Dict{AbstractString, Any}
                    @test length(j["menu"]) == 2
                    @test j["menu"]["header"] == "SVG\tViewerα"
                    @test isa(j["menu"]["items"], Array)
                    @test length(j["menu"]["items"]) == 22
                    @test j["menu"]["items"][3] == nothing
                    @test j["menu"]["items"][2]["id"] == "OpenNew"
                    @test j["menu"]["items"][2]["label"] == "Open New"
                end

validate_flickr(flickr) = begin
                              k = JSON.parse(flickr)
                              @test k != nothing
                              @test k["totalItems"] == 222
                              @test k["items"][1]["description"][12] == '\"'
                          end

validate_unicode(unicode) = begin
                                u = JSON.parse(unicode)
                                @test u != nothing
                                @test u["অলিম্পিকস"]["রেকর্ড"][2]["Marathon"] == "জনি হেইস"
                            end
# -------

finished_async_tests = RemoteRef()

@async begin
    s = listen(7777)
    s = accept(s)

    s.line_buffered = false
    Base.start_reading(s)

    @test JSON.parse(s) != nothing  # a
    @test JSON.parse(s) != nothing  # b
    validate_c(s)                     # c
    @test JSON.parse(s) != nothing  # d
    validate_e(s)                     # e
    @test JSON.parse(s) != nothing  # gmaps
    @test JSON.parse(s) != nothing  # colors1
    @test JSON.parse(s) != nothing  # colors2
    @test JSON.parse(s) != nothing  # colors3
    @test JSON.parse(s) != nothing  # twitter
    @test JSON.parse(s) != nothing  # facebook
    validate_flickr(s)                # flickr
    @test JSON.parse(s) != nothing  # youtube
    @test JSON.parse(s) != nothing  # iphone
    @test JSON.parse(s) != nothing  # customer
    @test JSON.parse(s) != nothing  # product
    @test JSON.parse(s) != nothing  # interop
    validate_unicode(s)               # unicode
    @test JSON.parse(s) != nothing  # issue5
    @test JSON.parse(s) != nothing  # dollars
    @test JSON.parse(s) != nothing  # brackets

    put!(finished_async_tests, nothing)
end

w = connect("localhost", 7777)

@test JSON.parse(a) != nothing
write(w, a)

@test JSON.parse(b) != nothing
write(w, b)

validate_c(c)
write(w, c)

@test JSON.parse(d) != nothing
write(w, d)

validate_e(e)
write(w, e)

@test JSON.parse(gmaps) != nothing
write(w, gmaps)

@test JSON.parse(colors1) != nothing
write(w, colors1)

@test JSON.parse(colors2) != nothing
write(w, colors2)

@test JSON.parse(colors3) != nothing
write(w, colors3)

@test JSON.parse(twitter) != nothing
write(w, twitter)

@test JSON.parse(facebook) != nothing
write(w, facebook)

validate_flickr(flickr)
write(w, flickr)

@test JSON.parse(youtube) != nothing
write(w, youtube)

@test JSON.parse(iphone) != nothing
write(w, iphone)

@test JSON.parse(customer) != nothing
write(w, customer)

@test JSON.parse(product) != nothing
write(w, product)

@test JSON.parse(interop) != nothing
write(w, interop)

validate_unicode(unicode)
write(w, unicode)

#Issue 5 on Github
issue5 = "[\"A\",\"B\",\"C\\n\"]"
JSON.parse(issue5)
write(w, issue5)

# $ escaping issue
dollars = ["all of the \$s", "µniçø∂\$"]
json_dollars = json(dollars)
@test JSON.parse(json_dollars) != nothing
write(w, json_dollars)

# unmatched brackets
brackets = @compat Dict("foo"=>"ba}r", "be}e]p"=>"boo{p")
json_brackets = json(brackets)
@test JSON.parse(json_brackets) != nothing
write(w, json_dollars)

fetch(finished_async_tests)

zeros = @compat Dict("\0"=>"\0")
json_zeros = json(zeros)
@test contains(json_zeros,"\\u0000")
@test !contains(json_zeros,"\\0")
@test JSON.parse(json_zeros) == zeros

#Uncomment while doing timing tests
#@time for i=1:100 ; JSON.parse(d) ; end

# Printing an empty array or Dict shouldn't cause a BoundsError
@test json(ASCIIString[]) == "[]"
@test json(Dict()) == "{}"

#test for issue 26
obj = JSON.parse("{\"a\":2e10}")
@test(obj["a"] == 2e10)

#test for issue 21
a=JSON.parse(test21)
@test isa(a, Array{Any})
@test length(a) == 2
#Multidimensional arrays
@test json([0 1; 2 0]) == "[[0,2],[1,0]]"

# ::Nothing values should be encoded as null
testDict = @compat Dict("a" => nothing)
nothingJson = JSON.json(testDict)
nothingDict = JSON.parse(nothingJson)
@test testDict == nothingDict

# test for issue #57
obj = JSON.parse("{\"\U0001d712\":\"\\ud835\\udf12\"}")
@test(obj["𝜒"] == "𝜒")

# test for single values
@test JSON.parse("true") == true
@test JSON.parse("null") == nothing
@test JSON.parse("\"hello\"") == "hello"
@test JSON.parse("\"a\"") == "a"
@test JSON.parse("1") == 1
@test JSON.parse("1.5") == 1.5

# test parsefile
tmppath, io = mktemp()
write(io, facebook)
close(io)
@unix_only @test haskey(JSON.parsefile(tmppath), "data")
# don't use mmap on Windows, to avoid ERROR: unlink: operation not permitted (EPERM)
@windows_only @test haskey(JSON.parsefile(tmppath; use_mmap=false), "data")
rm(tmppath)

# check indented json has same final value as non indented

fb = JSON.parse(facebook)
fbjson1 = json(fb, 2)
fbjson2 = json(fb)
@test JSON.parse(fbjson1) == JSON.parse(fbjson2)

ev = JSON.parse(e)
ejson1 = json(ev, 2)
ejson2 = json(ev)
@test JSON.parse(ejson1) == JSON.parse(ejson2)

# test symbols are treated as strings
symtest = @compat Dict(:symbolarray => [:apple, :pear], :symbolsingleton => :hello)
@test (JSON.json(symtest) == "{\"symbolarray\":[\"apple\",\"pear\"],\"symbolsingleton\":\"hello\"}"
         || JSON.json(symtest) == "{\"symbolsingleton\":\"hello\",\"symbolarray\":[\"apple\",\"pear\"]}")


text:
@testset "Test analyzer" begin
    @testset "Semantic analysis" begin
        for (code, expected) in [
            ("foobar", "identifier not found: foobar"),
            ("x = 2;", "identifier not found: x"),
            ("let a = 2; let a = 4;", "a is already defined"),
            ("let f = fn(x) { let x = 3; }", "x is already defined"),
            ("break;", "syntax error: break outside loop"),
            ("continue;", "syntax error: continue outside loop"),
            ("fn () { break; } ", "syntax error: break outside loop"),
            ("fn () { continue; } ", "syntax error: continue outside loop"),
            ("while (x == 2) {}", "identifier not found: x"),
            ("while (true) { x = 2; }", "identifier not found: x"),
            ("if (x == 2) {}", "identifier not found: x"),
            ("if (true) { x }", "identifier not found: x"),
            ("if (true) { 2 } else { x }", "identifier not found: x"),
            ("-x", "identifier not found: x"),
            ("x + 2", "identifier not found: x"),
            ("x(3)", "identifier not found: x"),
            ("x[3]", "identifier not found: x"),
            ("[1, 2][x]", "identifier not found: x"),
            ("len(x)", "identifier not found: x"),
            ("[x]", "identifier not found: x"),
            ("{x: 2}", "identifier not found: x"),
            ("{2: x}", "identifier not found: x"),
            (
                "let f = fn(x) { if (x > 0) { f(x - 1); f = 2; } }",
                "cannot reassign the current function being defined: f",
            ),
        ]
            test_object(m.analyze(code), expected)
        end
    end
end


text:
__precompile__()

module StataCall

    ##############################################################################
    ##
    #   TODO:
    #
    #   FUNCTIONALITY:
    #   -
    #
    #   TECHNICAL:
    #   - distinguish different versions of Stata, save in a particular version of
    #   the dta format to ensure readability
    #   - read/write of files to/from stata could be better (to ensure variable type)
    #   - use tempname() and tempdir() for temporary files
    #
    #   NOTE:
    #   This package is not written for speed -- if you need to write fast code, you
    #   should probably not use Stata in the first place. Instead, it should come in
    #   handy if you would want to use some of Stata's cleaning routines on fairly
    #   small DataFrames.
    #
    ##
    ##############################################################################

    ##############################################################################
    ##
    ## Dependencies
    ##
    ##############################################################################

    using FileIO, StatFiles, DataFrames, CSV, Dates

    ##############################################################################
    ##
    ## Exported methods and types
    ##
    ##############################################################################

    export stataCall

    ##############################################################################
    ##
    ## Load files
    ##
    ##############################################################################

    include("setup.jl")
    include("call.jl")

    executable  = locate_executable()

end # module


text:
using GeometricStructures
using Base.Test

import GeometricStructures.intersection
import GeometricStructures.convert
import GeometricStructures.convexhull

include("points.jl")
include("linesegments.jl")
include("convexhull.jl")

max_stars_repo_name: helgith/GeometricStructures

max_issues_repo_name: helgith/GeometricStructures

max_forks_repo_name: helgith/GeometricStructures


text:
## activate project environment
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

## load packages
using DelimitedFiles, Statistics
using OrdinaryDiffEq
using Trixi

## set up benchmark code
"""
    initial_condition_isentropic_vortex(x, t, equations)

The classical isentropic vortex test case of
- Chi-Wang Shu (1997)
  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
  Schemes for Hyperbolic Conservation Laws.
  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)
"""
function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
  ϱ0 = 1.0               # background density
  v0 = SVector(1.0, 1.0) # background velocity
  p0 = 10.0              # background pressure
  ε = 20.0               # vortex strength
  L = 10.0               # size of the domain per coordinate direction

  T0 = p0 / ϱ0           # background temperature
  γ = equations.gamma    # ideal gas constant

  x0 = v0 * t            # current center of the vortex
  dx = vortex_center.(x - x0, L)
  r2 = sum(abs2, dx)

  # perturbed primitive variables
  T = T0 - (γ - 1) * ε^2 / (8 * γ * π^2) * exp(1 - r2)
  v = v0 + ε / (2 * π) * exp(0.5 * (1 - r2)) * SVector(-dx[2], dx[1])
  ϱ = ϱ0 * (T / T0)^(1 / (γ - 1))
  p = ϱ * T

  return prim2cons(SVector(ϱ, v..., p), equations)
end

function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations3D)
  ϱ0 = 1.0               # background density
  v0 = SVector(1.0, 1.0, 0.0) # background velocity
  p0 = 10.0              # background pressure
  ε = 20.0               # vortex strength
  L = 10.0               # size of the domain per coordinate direction

  T0 = p0 / ϱ0           # background temperature
  γ = equations.gamma    # ideal gas constant

  x0 = v0 * t            # current center of the vortex
  dx = vortex_center.(x - x0, L)
  r2 = sum(abs2, dx)

  # perturbed primitive variables
  T = T0 - (γ - 1) * ε^2 / (8 * γ * π^2) * exp(1 - r2)
  v = v0 + ε / (2 * π) * exp(0.5 * (1 - r2)) * SVector(-dx[2], dx[1], 0.0)
  ϱ = ϱ0 * (T / T0)^(1 / (γ - 1))
  p = ϱ * T

  return prim2cons(SVector(ϱ, v..., p), equations)
end

vortex_center(x, L) = mod(x + L/2, L) - L/2

@noinline function flux_shima_etal_notinlined(u_ll, u_rrl, orientation_or_normal_direction, equations)
  flux_shima_etal(u_ll, u_rrl, orientation_or_normal_direction, equations)
end

@noinline function flux_ranocha_notinlined(u_ll, u_rrl, orientation_or_normal_direction, equations)
  flux_ranocha(u_ll, u_rrl, orientation_or_normal_direction, equations)
end

function run_benchmarks(polydeg, volume_flux, surface_flux, mesh)
  println("Benchmarking (", volume_flux, ", ", surface_flux, ") in ",
          ndims(mesh), "D with polydeg=", polydeg, " on ", nameof(typeof(mesh)))

  if ndims(mesh) == 2
    equations = CompressibleEulerEquations2D(1.4)
  else
    equations = CompressibleEulerEquations3D(1.4)
  end

  solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux,
                 volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

  semi = SemidiscretizationHyperbolic(mesh, equations,
                                      initial_condition_isentropic_vortex, solver)

  # compile
  ode = semidiscretize(semi, (0.0, 1.0))
  solve(ode, RDPK3SpFSAL49(), abstol=1.0e-7, reltol=1.0e-7,
        save_everystep=false, maxiters=1)

  pids = zeros(5)
  for i in eachindex(pids)
    sleep(0.5)
    take!(semi.performance_counter)
    solve(ode, RDPK3SpFSAL49(), abstol=1.0e-8, reltol=1.0e-8,
          save_everystep=false, maxiters=50)
    pid = 1.0e-9 * take!(semi.performance_counter) / Trixi.ndofs(semi)
    pids[i] = pid
  end

  return mean(pids), std(pids)
end

function run_benchmarks(polydegs, ndims::Integer, volume_flux, surface_flux)
  header = "# Polydeg"
  header = header * "\tTreeMesh-mean\tstd"
  header = header * "\tStructuredMesh-mean\tstd"
  header = header * "\tP4estMesh-mean\tstd"
  filename = joinpath(@__DIR__, "pids_$(ndims)D_$(volume_flux)_$(first(polydegs))_$(last(polydegs)).dat")
  open(filename, "w") do io
    println(io, header)
  end

  if ndims == 2
    coordinates_min = (-5.0, -5.0)
    coordinates_max = ( 5.0,  5.0)
    initial_refinement_level = 3
    cells_per_dimension = 2^initial_refinement_level .* (1, 1)
  else
    coordinates_min = (-5.0, -5.0, -5.0)
    coordinates_max = ( 5.0,  5.0,  5.0)
    initial_refinement_level = 3
    cells_per_dimension = 2^initial_refinement_level .* (1, 1, 1)
  end

  for polydeg in polydegs
    results_line = Float64[polydeg]

    # TreeMesh
    let mesh = TreeMesh(coordinates_min, coordinates_max,
                        initial_refinement_level=initial_refinement_level,
                        n_cells_max=100_000, periodicity=true)

      pid_mean, pid_std = run_benchmarks(polydeg, volume_flux, surface_flux, mesh)
      push!(results_line, pid_mean)
      push!(results_line, pid_std)
    end

    # StructuredMesh
    let mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
                              periodicity=true)

      pid_mean, pid_std = run_benchmarks(polydeg, volume_flux, surface_flux, mesh)
      push!(results_line, pid_mean)
      push!(results_line, pid_std)
    end

    # P4estMesh
    let mesh = P4estMesh(cells_per_dimension; coordinates_min, coordinates_max,
                         polydeg=1, periodicity=true)

      pid_mean, pid_std = run_benchmarks(polydeg, volume_flux, surface_flux, mesh)
      push!(results_line, pid_mean)
      push!(results_line, pid_std)
    end

    open(filename, "a") do io
      writedlm(io, results_line')
    end
  end

  return filename
end

## run benchmarks
polydegs = 3:15
# run_benchmarks(polydegs, 2, flux_shima_etal_notinlined, flux_shima_etal)
# run_benchmarks(polydegs, 2, flux_ranocha_notinlined, flux_ranocha)

# run_benchmarks(polydegs, 3, flux_shima_etal_notinlined, flux_shima_etal)
# run_benchmarks(polydegs, 3, flux_ranocha_notinlined, flux_ranocha)


text:
function trip_shutdown(path::String, day::Int) # day = 4 solo 4th, day = 1 tutti
    df = DataFrame()
    classes = Int[]
    attributes = CSV.File(path * "/Example_1.csv", drop = [1, 2]) |> DataFrame |> names

    for attr in attributes
        insertcols!(df, attr => Array{Float64, 1}[])
    end

    for row in CSV.File(path * "/DataInfo.csv")
        if row.Day >= day
            aux = CSV.File(path * "/Example_" * string(row.ExampleID) * ".csv", drop = [1, 2]) |> Tables.matrix

            push!(df, [aux[:, j] for j in 1:length(attributes)])

            push!(classes, row.Class)
        end
    end

    ds = ClassificationDataset([ModalFrame(df)], CategoricalArray(classes))

    return ds
end

function trip_no_trip(path::String; n_machine::Union{Missing, Int64} = missing)
    df = DataFrame()
    classes = Int[]
    attributes = CSV.File(path * "/Example_1.csv", drop = [1, 2]) |> DataFrame |> names

    for attr in attributes
        insertcols!(df, attr => Array{Float64, 1}[])
    end

    for row in CSV.File(path * "/DataInfo.csv")
        if row.Class == 0
            continue
        else
            if !ismissing(n_machine) && row.Datasource != n_machine
                continue
            end

            aux = CSV.File(path * "/Example_" * string(row.ExampleID) * ".csv", drop = [1, 2]) |> Tables.matrix

            push!(df, [aux[:, j] for j in 1:length(attributes)])

            if row.Day < 4
                push!(classes, 0)
            else
                push!(classes, 1)
            end
        end
    end

    ds = ClassificationDataset([ModalFrame(df)], CategoricalArray(classes))

    return ds
end