text:
function l1m1!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m+1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l+1,m+1,i-1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[3]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[1]+q[3]+1.0)
end

function lem1!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m+1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m+1,i+1]
		sum3 += p1[l-1,m+1,i-1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[3]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[2]+q[3]+1.0)
end

function l1me!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m-1,i+1]
		sum3 += p1[l+1,m-1,i-1]
		#ghost
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[4]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[1]+q[4]+1.0)
end

function leme!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m-1,i+1]
		sum3 += p1[l-1,m-1,i-1]
		#ghost
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[4]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[2]+q[4]+1.0)
end

function l1i1!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m-1,i+1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l+1,m-1,i+1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[5]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[1]+q[5]+1.0)
end

function lei1!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m-1,i+1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m+1,i+1]
		sum3 += p1[l-1,m-1,i+1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[5]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[2]+q[5]+1.0)
end

function leie!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m+1,i-1]
		sum3 += p1[l-1,m-1,i-1]
		#ghost
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[6]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[2]+q[6]+1.0)
end

function l1ie!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i-1]
		sum3 += p1[l+1,m-1,i-1]
		#ghost
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[6]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[1]+q[6]+1.0)
end

function m1i1!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l-1,m+1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l-1,m+1,i+1]
		#ghost
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[3]+q[5]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[3]+q[5]+1.0)
end

function m1ie!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i-1]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l-1,m+1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i-1]
		sum3 += p1[l-1,m+1,i-1]
		#ghost
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[3]+q[6]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[3]+q[6]+1.0)
end

function mei1!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l-1,m-1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m-1,i+1]
		sum3 += p1[l-1,m-1,i+1]
		#ghost
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[4]+q[5]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[4]+q[5]+1.0)
end

function meie!{T,N}(l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i-1]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m-1,i-1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l-1,m-1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m-1,i-1]
		sum3 += p1[l-1,m-1,i-1]
		#ghost
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[4]+q[6]-1.)*p0[l,m,i]
	p2[l,m,i] /= (q[4]+q[6]+1.0)
end

#adding source
function l1m1!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m+1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l+1,m+1,i-1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[3]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[1]+q[3]+1.0)
end

function lem1!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m+1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m+1,i+1]
		sum3 += p1[l-1,m+1,i-1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[3]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[2]+q[3]+1.0)
end

function l1me!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m-1,i+1]
		sum3 += p1[l+1,m-1,i-1]
		#ghost
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[4]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[1]+q[4]+1.0)
end

function leme!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l,m,i-1]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m-1,i+1]
		sum3 += p1[l-1,m-1,i-1]
		#ghost
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l,m,i+1]
		sum3 += p1[l,m,i-1]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[4]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[2]+q[4]+1.0)
end

function l1i1!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m-1,i+1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l+1,m-1,i+1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[5]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[1]+q[5]+1.0)
end

function lei1!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i+1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l,m-1,i+1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m+1,i+1]
		sum3 += p1[l-1,m-1,i+1]
		#ghost
		sum3 += p1[l,m+1,i+1]
		sum3 += p1[l,m-1,i+1]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[5]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[2]+q[5]+1.0)
end

function leie!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l-1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l-1,m+1,i]
		sum2 += p1[l-1,m-1,i]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l-1,m+1,i-1]
		sum3 += p1[l-1,m-1,i-1]
		#ghost
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[2]+q[6]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[2]+q[6]+1.0)
end

function l1ie!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l+1,m,i]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l,m,i-1]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l,m-1,i-1]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m-1,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i-1]
		sum3 += p1[l+1,m-1,i-1]
		#ghost
		sum3 += p1[l,m+1,i-1]
		sum3 += p1[l,m-1,i-1]
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l,m+1,i]
		sum3 += p1[l,m-1,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[1]+q[6]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[1]+q[6]+1.0)
end

function m1i1!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m+1,i+1]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l-1,m+1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i+1]
		sum3 += p1[l-1,m+1,i+1]
		#ghost
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[3]+q[5]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[3]+q[5]+1.0)
end

function m1ie!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i-1]
		sum1 += p1[l,m+1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m+1,i-1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l+1,m+1,i]
		sum2 += p1[l-1,m+1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m+1,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m+1,i-1]
		sum3 += p1[l-1,m+1,i-1]
		#ghost
		sum3 += p1[l+1,m+1,i]
		sum3 += p1[l-1,m+1,i]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[3]+q[6]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[3]+q[6]+1.0)
end

function mei1!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i+1]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m-1,i+1]
		sum2 += p1[l+1,m,i+1]
		sum2 += p1[l-1,m,i+1]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l-1,m-1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m,i+1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m-1,i+1]
		sum3 += p1[l-1,m-1,i+1]
		#ghost
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l+1,m,i+1]
		sum3 += p1[l-1,m,i+1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[4]+q[5]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[4]+q[5]+1.0)
end

function meie!{T,N}(n::Int,l::Int,m::Int,i::Int,
		   d1::T,d2::T,d3::T,d4::T,
		   p0::Array{T,N},p1::Array{T,N},p2::Array{T,N},q::Vector{T},s::Array{T,4})
	sum1,sum2,sum3,sum4 = 0.,0.,0.,0.
	if d1 != 0.
		sum1 += p1[l,m,i-1]
		sum1 += p1[l,m-1,i]
		sum1 += p1[l+1,m,i]
		sum1 += p1[l-1,m,i]
		#ghost
		sum1 += p1[l,m,i]
		sum1 += p1[l,m,i]
		sum1 *= d1
	end
	if d2 != 0.
		sum2 += p1[l,m-1,i-1]
		sum2 += p1[l+1,m,i-1]
		sum2 += p1[l-1,m,i-1]
		sum2 += p1[l+1,m-1,i]
		sum2 += p1[l-1,m-1,i]
		#ghost
		sum2 += p1[l,m,i]
		sum2 += p1[l,m-1,i]
		sum2 += p1[l,m,i-1]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 += p1[l+1,m,i]
		sum2 += p1[l-1,m,i]
		sum2 *= d2
	end
	if d3 != 0.
		sum3 += p1[l+1,m-1,i-1]
		sum3 += p1[l-1,m-1,i-1]
		#ghost
		sum3 += p1[l+1,m-1,i]
		sum3 += p1[l-1,m-1,i]
		sum3 += p1[l+1,m,i-1]
		sum3 += p1[l-1,m,i-1]
		sum3 += p1[l+1,m,i]
		sum3 += p1[l-1,m,i]
		sum3 *= d3
	end
	if d4 != 0.
		sum4 += d4*p1[l,m,i]
	end
	p2[l,m,i]  = sum1+sum2+sum3+sum4+(q[4]+q[6]-1.)*p0[l,m,i]+s[l,m,i,n]
	p2[l,m,i] /= (q[4]+q[6]+1.0)
end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
module PlotlyBase

using Base.Iterators
using JSON
using DocStringExtensions
using Requires
using UUIDs
using Dates
using Logging

import Base: ==

using Statistics: mean
using DelimitedFiles: readdlm

# import LaTeXStrings and export the handy macros
using LaTeXStrings
export @L_mstr, @L_str

# export some names from JSON
export json

_symbol_dict(x) = x
_symbol_dict(d::AbstractDict) =
    Dict{Symbol,Any}([(Symbol(k), _symbol_dict(v)) for (k, v) in d])

# include these here because they are used below
include("traces_layouts.jl")
include("styles.jl")

# core plot object
mutable struct Plot{TT<:AbstractTrace}
    data::Vector{TT}
    layout::AbstractLayout
    divid::UUID
    style::Style
end

# Default `convert` fallback constructor
Plot(p::Plot) = p

# include the rest of the core parts of the package
include("util.jl")
include("json.jl")
include("subplots.jl")
include("api.jl")
include("convenience_api.jl")
include("recession_bands.jl")
include("output.jl")

# Set some defaults for constructing `Plot`s
function Plot(;style::Style=CURRENT_STYLE[])
    Plot(GenericTrace{Dict{Symbol,Any}}[], Layout(), uuid4(), style)
end

function Plot(data::AbstractVector{T}, layout=Layout();
              style::Style=CURRENT_STYLE[]) where T<:AbstractTrace
    Plot(data, layout, uuid4(), style)
end

function Plot(data::AbstractTrace, layout=Layout();
              style::Style=CURRENT_STYLE[])
    Plot([data], layout; style=style)
end

# NOTE: we export trace constructing types from inside api.jl
# NOTE: we export names of shapes from traces_layouts.jl
export

    # core types
    Plot, GenericTrace, Layout, Shape, AbstractTrace, AbstractLayout,

    # plotly.js api methods
    restyle!, relayout!, update!, addtraces!, deletetraces!, movetraces!,
    redraw!, extendtraces!, prependtraces!, purge!, to_image, download_image,
    react!,

    # non-!-versions (forks, then applies, then returns fork)
    restyle, relayout, update, addtraces, deletetraces, movetraces, redraw,
    extendtraces, prependtraces, react,

    # helper methods
    plot, fork, vline, hline, attr,

    # new trace types
    stem,

    # convenience stuff
    add_recession_bands!,

    # styles
    use_style!, style, Style, Cycler, STYLES,

    # other
    savejson, savefig, savehtml, html_body

function __init__()
    env_style = Symbol(get(ENV, "PLOTLYJS_STYLE", ""))
    if env_style in STYLES
        global DEFAULT_STYLE
        DEFAULT_STYLE[] = Style(env_style)
    end
    @require IJulia="7073ff75-c697-5162-941a-fcdaad2a7d2a" begin
        function IJulia.display_dict(p::Plot)
            Dict(
                "application/vnd.plotly.v1+json" => JSON.lower(p),
                "text/plain" => sprint(show, "text/plain", p)
            )
        end
    end
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" include("dataframes_api.jl")
    @require Distributions="31c24e10-a181-5473-b8eb-7969acd0382f" include("distributions.jl")
    @require Colors="5ae59095-9a9b-59fe-a467-6f913c188581" JSON.lower(a::Colors.Colorant) = string("#", Colors.hex(a))
end

end # module

ext: jl


text: const GLOBAL_RNG = Random.MersenneTwister()


text:
function grad!(g, fset::FittingSet{T}, prob::NonlinearProblem, targ::TargetFunction{T}, μ, n) where {T}
    @unpack cache = fset
    @unpack gbuf = cache

    # Calculation of the gradient for each group of parameters
    fill!(gbuf, zero(T))
    # FIXED: Using 'map' instead of 'for' preserves 'type-stability'
    map(fset.params) do param
        grad!(cache, param, prob, targ, μ, n)
    end
    @turbo for i = 1:length(g)
        g[i] += gbuf[i]
    end

    return g
end

function grad!(cache::FittingCache, param::FittingParameter{:Tconn}, prob::NonlinearProblem, targ::TargetFunction, μ, n)
    @unpack C = prob
    @unpack gviews, bviews = param
    @unpack Pcalc = @inbounds prob.pviews[n]
    @unpack cbuf, cbuf2 = cache
    V = @inbounds param.vviews[n]

    mul!(cbuf, C, Pcalc)
    mul!(cbuf2, C, μ)
    @turbo for i = 1:length(cbuf)
        cbuf[i] *= cbuf2[i]
    end
    @inbounds @simd for i = 1:length(bviews)
        gviews[i] = bviews[i] * V[i]
    end

    return cache
end

function grad!(cache::FittingCache, param::FittingParameter{:Tconst}, prob::NonlinearProblem, targ::TargetFunction, μ, n)

    @unpack gviews, bviews = param
    @unpack Pcalc, Pi = @inbounds prob.pviews[n]
    @unpack tbuf = cache
    V = @inbounds param.vviews[n]

    @turbo for i = 1:length(tbuf)
        tbuf[i] = (Pcalc[i] - Pi[i]) * μ[i]
    end
    @inbounds @simd for i = 1:length(bviews)
        gviews[i] = bviews[i] * V[i]
    end

    return cache
end

function grad!(cache::FittingCache{T}, param::FittingParameter{:Vpi, T}, prob::NonlinearProblem{T}, targ::TargetFunction{T}, μ, n) where {T}

    @unpack tbuf, tbuf2 = cache
    @unpack "ⁿ", Pcalc, Pi, Swi, Bwi, Boi, cw, co, cf = @inbounds prob.pviews[n]
    @unpack "ⁿ⁻¹", Pcalc, Pi, Swi, Bwi, Boi, cw, co, cf, Δt = @inbounds prob.pviews[n-1]
    @unpack gviews, bviews, bviews2 = param
    Vⁿ = @inbounds param.vviews[n]
    Vⁿ⁻¹ = @inbounds param.vviews[n-1]

    @turbo for i = 1:length(tbuf)
        ΔP = Pcalcⁿ[i] - Piⁿ[i]
        tbuf[i] = Swiⁿ[i] / Bwiⁿ[i] * exp((cwⁿ[i] + cfⁿ[i]) * ΔP)
        tbuf[i] += (one(T) - Swiⁿ[i]) / Boiⁿ[i] * exp((coⁿ[i] + cfⁿ[i]) * ΔP)
        tbuf[i] *= μ[i] / Δtⁿ⁻¹[]
        ΔP = Pcalcⁿ⁻¹[i] - Piⁿ⁻¹[i]
        tbuf2[i] = Swiⁿ⁻¹[i] / Bwiⁿ⁻¹[i] * exp((cwⁿ⁻¹[i] + cfⁿ⁻¹[i]) * ΔP)
        tbuf2[i] += (one(T) - Swiⁿ⁻¹[i]) / Boiⁿ⁻¹[i] * exp((coⁿ⁻¹[i] + cfⁿ⁻¹[i]) * ΔP)
        tbuf2[i] *= μ[i] / Δtⁿ⁻¹[]
    end
    @inbounds @simd for i = 1:length(bviews)
        gviews[i] = bviews[i] * Vⁿ[i] - bviews2[i] * Vⁿ⁻¹[i]
    end

    return cache
end

function grad!(cache::FittingCache{T}, param::FittingParameter{:cf, T}, prob::NonlinearProblem{T}, targ::TargetFunction{T}, μ, n) where {T}

    @unpack tbuf, tbuf2 = cache
    @unpack "ⁿ", Pcalc, Vpi, Pi, Swi, Bwi, Boi, cw, co, cf = @inbounds prob.pviews[n]
    @unpack "ⁿ⁻¹", Pcalc, Vpi, Pi, Swi, Bwi, Boi, cw, co, cf, Δt = @inbounds prob.pviews[n-1]
    @unpack gviews, bviews, bviews2 = param
    Vⁿ = @inbounds param.vviews[n]
    Vⁿ⁻¹ = @inbounds param.vviews[n-1]

    @turbo for i = 1:length(tbuf)
        ΔP = Pcalcⁿ[i] - Piⁿ[i]
        tbuf[i] = Swiⁿ[i] * Vpiⁿ[i] / Bwiⁿ[i] * exp((cwⁿ[i] + cfⁿ[i]) * ΔP) * ΔP
        tbuf[i] += (one(T) - Swiⁿ[i]) * Vpiⁿ[i] / Boiⁿ[i] * exp((coⁿ[i] + cfⁿ[i]) * ΔP) * ΔP
        tbuf[i] *= μ[i] / Δtⁿ⁻¹[]
        ΔP = Pcalcⁿ⁻¹[i] - Piⁿ⁻¹[i]
        tbuf2[i] = Swiⁿ⁻¹[i] * Vpiⁿ⁻¹[i] / Bwiⁿ⁻¹[i] * exp((cwⁿ⁻¹[i] + cfⁿ⁻¹[i]) * ΔP) * ΔP
        tbuf2[i] += (one(T) - Swiⁿ⁻¹[i]) * Vpiⁿ⁻¹[i] / Boiⁿ⁻¹[i] * exp((coⁿ⁻¹[i] + cfⁿ⁻¹[i]) * ΔP) * ΔP
        tbuf2[i] *= μ[i] / Δtⁿ⁻¹[]
    end
    @inbounds @simd for i = 1:length(bviews)
        gviews[i] = bviews[i] * Vⁿ[i] - bviews2[i] * Vⁿ⁻¹[i]
    end

    return cache
end

function grad!(cache::FittingCache, param::FittingParameter{:λ}, prob::NonlinearProblem, targ::TargetFunction, μ, n)

    @unpack gviews, bviews = param
    @unpack Qinj = @inbounds prob.pviews[n]
    @unpack tbuf = cache
    V = @inbounds param.vviews[n]
    gλ = @inbounds targ.terms.Pinj.gλviews[n]

    @turbo for i = 1:length(tbuf)
        tbuf[i] = -Qinj[i] * μ[i] + gλ[i]
    end
    @inbounds @simd for i = 1:length(bviews)
        gviews[i] = bviews[i] * V[i]
    end

    return cache
end

function grad!(g, term::L2TargetTerm)
    @turbo for i = 1:length(g)
        g[i] += term.g[i]
    end
    return g
end


text:
struct URange{T<:Real} <: AbstractUnitRange{T}
    start::T
    stop::T
    URange{T}(start, stop) where {T} = new{T}(start, urange_last(start,stop))
end
URange(start::T, stop::T) where {T<:Real} = URange{T}(start, stop)

urange_last(::Bool, stop::Bool) = stop
urange_last(start::T, stop::T) where {T<:Integer} =
    ifelse(stop >= start, stop, convert(T,start-one(stop-start)))
urange_last(start::T, stop::T) where {T} =
    ifelse(stop >= start, convert(T,start+floor(stop-start)),
                          convert(T,start-one(stop-start)))

function Base.iterate(r::URange{T}) where {T}
    r.start > r.stop && return nothing
    x = oftype(r.start + one(T), r.start)
    return (x,x)
end
function Base.iterate(r::URange{T}, i) where T
    x = i + one(T)
    return (x > oftype(x, r.stop) ? nothing : (x,x))
end

Base.intersect(r::URange{T1}, s::URange{T2}) where {T1<:Integer,T2<:Integer} = URange(max(first(r),first(s)), min(last(r),last(s)))

@inline function Base.getindex(r::URange{R}, s::AbstractUnitRange{S}) where {R,S<:Integer}
    @boundscheck checkbounds(r, s)
    f = first(r)
    strt = f + first(s) - 1
    URange{R}(strt, strt+length(s)-1)
end

Base.promote_rule(::Type{URange{T1}},::Type{URange{T2}}) where {T1,T2} =
    URange{promote_type(T1,T2)}
Base.convert(::Type{URange{T}}, r::URange{T}) where {T<:Real} = r
Base.convert(::Type{URange{T}}, r::URange) where {T<:Real} = URange{T}(r.start, r.stop)

Base.promote_rule(::Type{URange{T1}}, ::Type{UR}) where {T1,UR<:AbstractUnitRange} =
    URange{promote_type(T1,eltype(UR))}
Base.promote_rule(::Type{UnitRange{T2}}, ::Type{URange{T1}}) where {T1,T2} =
    URange{promote_type(T1,T2)}
Base.convert(::Type{URange{T}}, r::AbstractUnitRange) where {T<:Real} = URange{T}(first(r), last(r))
URange{T}(r::AbstractUnitRange) where {T} = convert(Type{URange{T}}, r)

Base.show(io::IO, r::URange) = print(io, typeof(r).name, '(', repr(first(r)), ',', repr(last(r)), ')')


text:
using MCHammer
using Distributions
using Random

if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "Correlating Values" begin
    Random.seed!(1)
    TestMatrix = rand(Normal(),1000,5)
     test_r = cormat(TestMatrix)
     @test test_r[1,2] == 0.04501199701199701
 end

 @testset "Covariance Values" begin
     Random.seed!(1)
     TestMatrix = rand(Normal(),1000,5)
      test_r = covmat(TestMatrix)
      @test test_r[5,5] == 1.027673689237494
end

@testset "Correlating Distributions" begin
    Random.seed!(1)
    n_trials = 1000
    sample_data = [rand(LogNormal(0, 0.5),n_trials) rand(Normal(3,2),n_trials) rand(Gamma(1, 0.5),n_trials) rand(LogNormal(0, 0.5),n_trials) rand(Normal(3,2),n_trials) rand(Gamma(1, 0.5),n_trials)]

    test_cmatrix = [1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0;0 0 0 1 0.75 -0.7; 0 0 0 0.75 1 -0.95; 0 0 0 -0.7 -0.95 1 ]

    Random.seed!(1)
    test_r = cormat(corvar(sample_data, n_trials, test_cmatrix))
    @test test_r[3,4] == 0.01943957543957544

end

@testset "GetCertainty" begin
    Random.seed!(1)
    test = rand(Normal(),1000)
    @test GetCertainty(test, 0, 1) == 0.502
end

# @testset "fractiles" begin
#     Random.seed!(1)
#     test = rand(Normal(),1000)
#     @test fractiles(test) == "11×2 Array{Any,2}:
#  "P0.0"    -3.882
#  "P10.0"   -1.34325
#  "P20.0"   -0.860748
#  "P30.0"   -0.526089
#  "P40.0"   -0.274806
#  "P50.0"    0.00474446
#  "P60.0"    0.218685
#  "P70.0"    0.472055
#  "P80.0"    0.800966
#  "P90.0"    1.2504
#  "P100.0"   3.12432"
# end

@testset "truncate" begin
    Result_1 = truncate_digit(0.667)
    Result_2 = truncate_digit(0.661)
    @test  Result_1 == Result_2
end

@testset "GBMM" begin
    Random.seed!(1)
    test_r = GBMM(100000, 0.05,0.05,12)
    @test test_r[12,1] == 195258.14301210243
end

@testset "GBMMfit" begin
    Random.seed!(1)
    historical = rand(Normal(10,2.5),1000)
    test_r = GBMMfit(historical, 12)
    @test test_r[12,1] == 2.1886864479965875

end

@testset "RiskEvent" begin
    test_r =RiskEvent(0.3, Normal(), 10,seed=1)
    @test round(test_r[9,1];digits=6) ==  -0.693654

end

@testset "MarkovChains" begin
#Setup test data
    Marital_StatM = [0.85	0.12	0.02	0.01;
    0	0.88	0.08	0.04;
    0	0.13	0.45	0.42;
    0	0.09	0.02	0.89;
    ]

    initial_marital = [0.43, 0.29, 0.11, 0.17 ]

AnalyticSolution = markov_a(Marital_StatM)
@test round(AnalyticSolution[4,1]; digits=6) == 0.473998

TS_Solution = markov_ts(Marital_StatM, initial_marital, 5)

@test round(TS_Solution[4,5]; digits=6) == 0.305876
end


text:
#### Main

"""
    run(case)
Solve the Eddy-Diffusivity Mass-Flux (EDMF) equations for a
stand-alone `case`
"""
function run(param_set, case, output_dir)
    params = Params(param_set, case)

    tc = TurbConv(params, case, output_dir)

    grid = tc.grid
    q = tc.q
    q_new = tc.q_new
    aux = tc.aux
    q_tendencies = tc.tendencies
    tri_diag = tc.tri_diag
    aux_O2 = tc.aux_O2

    gm, en, ud, sd, al = allcombinations(q)

    # Initialize nc files
    nc_ics_q = NetCDFWriter(joinpath(output_dir, "ics_q"))
    nc_ics_q = init_data(nc_ics_q, grid, q)
    nc_ics_aux = NetCDFWriter(joinpath(output_dir, "ics_aux"))
    nc_ics_aux = init_data(nc_ics_aux, grid, aux)

    # Initialize reference state
    init_ref_state!(aux, grid, params, output_dir)
    append_data(nc_ics_aux, aux, 0)

    # Initialize forcing
    init_forcing!(q, aux, grid, params, output_dir, case)
    append_data(nc_ics_aux, aux, 0)

    # Initialize solution
    init_state_vecs!(q, aux, grid, params, output_dir, case)
    append_data(nc_ics_q, q, 0)

    params[:UpdVar] =
        [UpdraftVar(0, params[:SurfaceModel].area, length(ud)) for i in al]

    @unpack Δt, t_end = params

    i_Δt, i_export, t = [0], [0], [0.0]

    nc_q = NetCDFWriter(joinpath(output_dir, "prog_vs_time"))
    nc_q = init_data(nc_q, grid, q)
    nc_aux = NetCDFWriter(joinpath(output_dir, "aux_vs_time"))
    nc_aux = init_data(nc_aux, grid, aux)

    @show joinpath(pwd(), full_name(nc_q))

    ∑tendencies_params = (; params, grid, aux_O2, aux, case, tri_diag)
    while t[1] < t_end

        update_aux!(grid, q, aux, aux_O2, params, case)

        ∑tendencies!(q_tendencies, q, ∑tendencies_params, t[1])

        solve!(q_new, q, q_tendencies, ∑tendencies_params)

        update_dt!(grid, params, q, t)

        i_export[1] += 1
        if mod(i_export[1], params[:export_frequency]) == 0
            @info "Exporting file: $(nc_q.filename)"
            append_data(nc_q, q, t[1])
            append_data(nc_aux, aux, t[1])
        end
    end

    append_data(nc_q, q, t[1])
    append_data(nc_aux, aux, t[1])

    return (grid, q, aux, params)
end


text:
using PatModules
@mainmodule! subpackage ("one", "two") begin
    export say_hi, say_bye
end


text:
using Test
using RigidBodyDynamics
import Random

# A pathologically weird matrix which uses base -1 indexing
# for its first dimension and base 2 indexing for its second
struct NonOneBasedMatrix <: AbstractMatrix{Float64}
    m::Int
    n::Int
end

Base.size(m::NonOneBasedMatrix) = (m.m, m.n)
Base.axes(m::NonOneBasedMatrix) = ((1:m.m) .- 2, (1:m.n) .+ 1)

@testset "custom collections" begin
    @testset "nulldict" begin
        nd = RigidBodyDynamics.NullDict{Int, Int}()
        @test isempty(nd)
        @test length(nd) == 0
        for element in nd
            @test false # should never be reached, since the nulldict is always empty
        end
        show(IOBuffer(), nd)
    end

    @testset "IndexDict" begin
        Int32Dict{V} = RigidBodyDynamics.IndexDict{Int32, Base.OneTo{Int32}, V}
        dict = Dict(Int32(2) => 4., Int32(1) => 3.)
        expected = Int32Dict{Float64}(Base.OneTo(Int32(2)), [3., 4.])
        i32dict1 = @inferred Int32Dict(dict)
        @test i32dict1 == dict == expected
        @test keys(expected) === keys(i32dict1)
        @test values(expected) == values(i32dict1)
        i32dict2 = @inferred Int32Dict{Float64}(dict)
        @test i32dict2 == dict == expected
        @test keys(expected) === keys(i32dict2)
        @test values(expected) == values(i32dict2)
    end

    @testset "ConstDict" begin
        c = RigidBodyDynamics.ConstDict{Int}(2.0)
        @test c[1] == 2.0
        @test c[-1] == 2.0
        show(IOBuffer(), c)
    end

    @testset "SegmentedVector" begin
        x = [1., 2., 3., 4.]
        viewlength = i -> 2
        xseg = SegmentedVector{Int}(x, 1 : 2, viewlength)
        @test segments(xseg)[1] == [1., 2.]
        @test segments(xseg)[2] == [3., 4.]
        @test length(segments(xseg)) == 2
        yseg = similar(xseg, Int32)
        yseg .= 1 : 4
        @test segments(yseg)[1] == [1, 2]
        @test segments(yseg)[2] == [3, 4]
        @test length(segments(yseg)) == 2

        xseg2 = SegmentedVector(x, RigidBodyDynamics.IndexDict(Base.OneTo(2), [view(x, 1 : 3), view(x, 4 : 4)]))
        @test xseg2 isa SegmentedVector
        ranges2 = RigidBodyDynamics.CustomCollections.ranges(xseg2)
        @test ranges2[1] == 1 : 3
        @test ranges2[2] == 4 : 4
        @test xseg2 == xseg

        xseg3 = copy(xseg)
        @test xseg3 == xseg
        @test xseg3 isa SegmentedVector
    end

    @testset "SegmentedBlockDiagonalMatrix" begin
        Random.seed!(5)
        A = rand(10, 10)
        block_indices = [(1:1, 1:1),  # square
                         (2:4, 2:2),  # non-square
                         (5:4, 3:2),  # empty
                         (5:7, 3:6),  # 2x2
                         (8:10, 7:10)]
        RigidBodyDynamics.CustomCollections.check_contiguous_block_ranges(A, block_indices)
        S = RigidBodyDynamics.SegmentedBlockDiagonalMatrix(A, block_indices)
        for (i, block) in enumerate(block_indices)
            @test S[block...] == RigidBodyDynamics.CustomCollections.blocks(S)[i]
        end
        A .= 0
        for block in RigidBodyDynamics.CustomCollections.blocks(S)
            Random.rand!(block)
        end

        @testset "Malformed blocks" begin
            @testset "overlap" begin
            block_indices = [(1:1, 1:1),
                             (2:4, 2:3),
                             (5:4, 3:2),
                             (5:7, 3:6),
                             (8:10, 7:10)]
                @test_throws ArgumentError RigidBodyDynamics.CustomCollections.check_contiguous_block_ranges(A, block_indices)
            end

            @testset "out of bounds" begin
                block_indices = [(1:1, 0:1),
                                 (2:4, 2:2),
                                 (5:4, 3:2),
                                 (5:7, 3:6),
                                 (8:10, 7:10)]
                @test_throws ArgumentError RigidBodyDynamics.CustomCollections.check_contiguous_block_ranges(A, block_indices)

                block_indices = [(1:1, 1:1),
                                 (2:4, 2:2),
                                 (5:4, 3:2),
                                 (5:7, 3:6),
                                 (8:12, 7:10)]
                @test_throws ArgumentError RigidBodyDynamics.CustomCollections.check_contiguous_block_ranges(A, block_indices)
            end

            @testset "gap" begin
                block_indices = [(1:1, 1:1),
                                 (5:4, 3:2),
                                 (5:7, 3:6),
                                 (8:10, 7:10)]
                @test_throws ArgumentError RigidBodyDynamics.CustomCollections.check_contiguous_block_ranges(A, block_indices)
            end
        end

        @testset "Nonstandard indexing" begin
            M = NonOneBasedMatrix(5, 5)
            block_indices = [(-1:1, 2:3),
                             (2:3,  4:6)]
            RigidBodyDynamics.CustomCollections.check_contiguous_block_ranges(M, block_indices)
        end

        @testset "mul! specialization" begin
            M = rand(20, size(S, 1))
            C = M * S
            @test C ≈ M * parent(S) atol=1e-15

            M = rand(size(S, 2), 20)
            C = S * M
            @test C ≈ parent(S) * M atol=1e-15

            @test_throws DimensionMismatch rand(20, size(S, 1) + 1) * M
            @test_throws DimensionMismatch rand(size(A, 2) + 1, 20) * M
        end
    end

    @testset "UnorderedPair" begin
        UnorderedPair = RigidBodyDynamics.CustomCollections.UnorderedPair
        p1 = UnorderedPair(1, 2)
        p2 = UnorderedPair(2, 1)
        @test p1 == p2
        @test hash(p1) == hash(p2)
        @test p1 != UnorderedPair(3, 1)
        dict = Dict(p1 => 3)
        @test dict[p2] == 3
    end
end


text:
using Test, Random

@testset "basic" begin
    @testset "Chain" begin
        @test_nowarn Chain(Dense(10, 5, σ), Dense(5, 2))(randn(10))
        @test_throws DimensionMismatch Chain(Dense(10, 5, σ),Dense(2, 1))(randn(10))
        # numeric test should be put into testset of corresponding layer
    end

    @testset "Dense" begin
        @test  length(Dense(10, 5)(randn(10))) == 5
        @test_throws DimensionMismatch Dense(10, 5)(randn(1))
        @test_throws MethodError Dense(10, 5)(1) # avoid broadcasting
        @test_throws MethodError Dense(10, 5).(randn(10)) # avoid broadcasting

        @test Dense(10, 1, identity, initW = ones, initb = zeros)(ones(10,1)) == 10*ones(1, 1)
        @test Dense(10, 1, identity, initW = ones, initb = zeros)(ones(10,2)) == 10*ones(1, 2)
        @test Dense(10, 2, identity, initW = ones, initb = zeros)(ones(10,1)) == 10*ones(2, 1)
        @test Dense(10, 2, identity, initW = ones, initb = zeros)([ones(10,1) 2*ones(10,1)]) == [10 20; 10 20]

    end

    @testset "Diagonal" begin
        @test length(Flux.Diagonal(10)(randn(10))) == 10
        @test length(Flux.Diagonal(10)(1)) == 10
        @test length(Flux.Diagonal(10)(randn(1))) == 10
        @test_throws DimensionMismatch Flux.Diagonal(10)(randn(2))

        @test Flux.Diagonal(2)([1 2]) == [1 2; 1 2]
        @test Flux.Diagonal(2)([1,2]) == [1,2]
        @test Flux.Diagonal(2)([1 2; 3 4]) == [1 2; 3 4]
    end
end


text:
using Statistics,Serialization,StatsPlots

# parameter values
paras = ["beta_l5","alpha_o050","beta_l5_alpha_l025", "baseline_GER", "beta_l5_alpha_l025_alpha_o025"]

filename_prefix = "../../figures/working_paper_2/"

path_to_data_prefix = "../../data/working_paper_2/"

mkpath(filename_prefix)

include("$(path_to_data_prefix)/baseline_GER/covid_par_ini.jl")
datapoint = fld(T,datat)

for (i,p) in enumerate(paras)
    local path_to_data = string(path_to_data_prefix, p , "/")

    local worker_results = deserialize(open("$(path_to_data)batchdata.dat"))

    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    gdpstd = std(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    gdpall = hcat((results -> results[:gdppercaptraj]).(worker_results)...)
    local nn = size(gdpmean)[1]
    pl1a = plot(datat*collect(2:nn),[gdpall[2:nn,1],gdpall[2:nn,2],gdpall[2:nn,3],gdpall[2:nn,4],gdpall[2:nn,5],gdpall[2:nn,6]], linewidth = [2 2 2 2 2 2], linecolor = [:red :green :black :blue :cyan :grey],ylims = (32,48), xlabel = "Days", ylabel = "GDP", label = false)
    savefig(pl1a,"$(filename_prefix)gdp_dynamics_single_$(p).pdf")

    local hcapfrac = zeros(datapoint+1)
    for th = 1:datapoint+1
        hcapfrac[th] = hcap / (nhh * icufrac)
    end
    infmean = mean(cat((results -> results[:inftraj]).(worker_results)...,dims=3), dims=3)
    infstd = std(cat((results -> results[:inftraj]).(worker_results)...,dims=3), dims=3)
    infall = cat((results -> results[:inftraj]).(worker_results)...,dims=3)
    pl6a = plot(datat*collect(2:nn),[infall[2:nn,3,1],infall[2:nn,3,2],infall[2:nn,3,3],infall[2:nn,3,4],infall[2:nn,3,5],infall[2:nn,3,6]], linewidth = [2 2 2 2 2 2], linecolor = [:red :green :black :blue :cyan :grey], ylabel = "Inf [%]", xlabel = "Days", label = false)
    savefig(pl6a,"$(filename_prefix)infection_dynamics_single_$(p).pdf")

end


text:
using Revise
using Printf
using Random
using POMDPs
using POMDPPolicies
using POMDPSimulators
using BeliefUpdaters
using AutomotiveDrivingModels
using AutomotivePOMDPs
using AutoViz
using RLInterface
using DeepQLearning

rng = MersenneTwister(1)

const K = 4

pomdp = OCPOMDP(ΔT = 0.5, p_birth = 0.3, max_peds = 1)

env = KMarkovEnvironment(pomdp, k=K)

up = KMarkovUpdater(K)

policy = RandomPolicy(pomdp, rng=rng)

hr = HistoryRecorder(rng=rng, max_steps = 100)
s0 = rand(rng, initialstate(pomdp))
initial_observation = rand(rng, initialobs(pomdp, s0))
initial_obs_vec = fill(initial_observation, 4)
hist = POMDPs.simulate(hr, pomdp, policy, up, initial_obs_vec, s0)

# Visualize
using Reel
frames = Frames(MIME("image/png"), fps=4)
AutoViz.render(s0, pomdp.env, cam = StaticCamera(VecE2(25.0,0.0), 15.0))
for step in eachstep(hist, "s,a,r,sp")
    s, a, r, sp = step
    push!(frames, AutoViz.render(sp, pomdp.env, cam=StaticCamera(VecE2(25.0,0.0), 15.0)))
end
write("out.gif", frames)

pomdp.env.obstacles


text:
using JuMP
using Gurobi
using Plasmo
using PlasmoAlgorithms

mp = Model(solver = GurobiSolver())
sp1 = Model(solver = GurobiSolver())
sp2 = Model(solver = GurobiSolver())
sp3 = Model(solver = GurobiSolver())

@variable(mp,x[1:3]>=0)
@constraint(mp,x[1]+x[2]+x[3]<=500)
@objective(mp,Min,150x[1]+230x[2]+260x[3])

@variable(sp1, w[1:4]>=0)
@variable(sp1, y[1:2]>=0)
@variable(sp1, x[1:3]>=0)
@constraint(sp1, 3x[1]+y[1]-w[1]>=200)
@constraint(sp1, 3.6x[2]+y[2]-w[2]>= 240)
@constraint(sp1, w[3]+w[4]<=24x[3])
@constraint(sp1, w[3]<=6000)
@objective(sp1,Min,(-1/3)*(170w[1]-238y[1]+150w[2]-210y[2]+36w[3]+10w[4]))

@variable(sp2, w[5:8]>=0)
@variable(sp2, y[3:4]>=0)
@variable(sp2, x[1:3]>=0)
@constraint(sp2, 2.5x[1]+y[3]-w[5]>=200)
@constraint(sp2, 3x[2]+y[4]-w[6]>=240)
@constraint(sp2, w[7]+w[8]<=20x[3])
@constraint(sp2, w[7]<=6000)
@objective(sp2, Min, (-1/3)*(170w[5]-238y[3]+150w[6]-210y[4]+36w[7]+10w[8]))

@variable(sp3,w[9:12]>=0)
@variable(sp3,y[5:6]>=0)
@variable(sp3, x[1:3]>=0)
@constraint(sp3, 2x[1]+y[5]-w[9]>=200)
@constraint(sp3, 2.4x[2]+y[6]-w[10]>=240)
@constraint(sp3, w[11]+w[12]<=16x[3])
@constraint(sp3, w[11]<=6000)
@objective(sp3, Min, (-1/3)*(170w[9]-238y[5]+150w[10]-210y[6]+36w[11]+10w[12]))

g = ModelGraph()
setsolver(g, GurobiSolver())
n1 = add_node(g)
n2 = add_node(g)
n3 = add_node(g)
n4 = add_node(g)

setmodel(n1, mp)
setmodel(n2, sp1)
setmodel(n3, sp2)
setmodel(n4, sp3)

edge1 = Plasmo.add_edge(g,n1,n2)
edge2 = Plasmo.add_edge(g,n1,n3)
edge3 = Plasmo.add_edge(g,n1,n4)

@linkconstraint(g,[i in 1:3], n1[:x][i] == n2[:x][i])
@linkconstraint(g,[i in 1:3], n1[:x][i] == n3[:x][i])
@linkconstraint(g,[i in 1:3], n1[:x][i] == n4[:x][i])

r = bendersolve(g,max_iterations = 10)


text:
module GLMakie

using ModernGL, FixedPointNumbers, Colors, GeometryBasics, StaticArrays
using AbstractPlotting, FileIO

using AbstractPlotting: @key_str, Key, broadcast_foreach, to_ndim, NativeFont
using AbstractPlotting: Scene, Lines, Text, Image, Heatmap, Scatter
using AbstractPlotting: convert_attribute, @extractvalue, LineSegments
using AbstractPlotting: @get_attribute, to_value, to_colormap, extrema_nan
using AbstractPlotting: ClosedInterval, (..)
using ShaderAbstractions
using FreeTypeAbstraction

using Base: RefValue
import Base: push!, isopen, show
using Base.Iterators: repeated, drop

using LinearAlgebra

for name in names(AbstractPlotting)
    @eval import AbstractPlotting: $(name)
    @eval export $(name)
end

struct GLBackend <: AbstractPlotting.AbstractBackend
end

loadshader(name) = normpath(joinpath(@__DIR__, "..", "assets", "shader", name))

# don't put this into try catch, to not mess with normal errors
include("gl_backend.jl")

function activate!(use_display=true)
    b = GLBackend()
    AbstractPlotting.register_backend!(b)
    AbstractPlotting.set_glyph_resolution!(AbstractPlotting.High)
    AbstractPlotting.current_backend[] = b
    AbstractPlotting.inline!(!use_display)
end

function __init__()
    activate!()
end

export set_window_config!

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end


text:
# Script to create leontief inverse matrix, exports and value added

using DataFrames, XLSX, LinearAlgebra

# ---------- Notes -------------------------------------------------------------------------------------------------------

# Problems with value added. See comment for v_add_share.

# ---------- Load data -------------------------------------------------------------------------------------------------------

# Subsample of EU27+GB (code also works for original dataset too only need to adjust the row bounds from 6/5 to 7/6 and N = 43+1)
path = "clean/raw_wide_EU.xlsx"

raw = DataFrame(XLSX.readxlsx(path)["Sheet1"][:], :auto)

# ---------- Data wrangling -------------------------------------------------------------------------------------------------------

N = 27 + 1 # number of countries, EU27+GB
S = 56 # number of sectors

# cut out the additional variables in beginning/end (to have square matrix) (N*S×N*S rows/columns)
raw_D = raw[6:N*S+5, 5:N*S+4] # intermediary demand data
# cut sectoral input/output so only final demand variables remain (S*N rows and 5*N columns - 5 different measures of final demand:
# households/government/NGOs/inventory/capital_formation)
raw_FD = raw[6:N*S+5, 5+N*S:4+N*S+5*N] # final demand data

# convert into matrix (cannot make matrix calculations with DataFrames)
D = Matrix(convert.(Float64, raw_D)) # N*S×N*S
FD = Matrix(convert.(Float64, raw_FD)) # N*S×N*5

# ---------- Variable creation -------------------------------------------------------------------------------------------------------

# gross total inputs used by sector k and country i: sum rows of column k,i
inputs = [sum(D[1:end, i]) for i in 1:size(D, 2)] # N*S×1

# gross total output from sector k and country i: sum columns of row k,i
outputs = [sum(D[i, 1:end]) for i in 1:size(D, 1)] # N*S×1, without final demand

# gross final demand from sector k and country i per country: sum five final demand columns of row k, i
f_demand = [sum(FD[i, j:j+4]) for i in 1:size(FD, 1), j in 1:5:size(FD, 2)] # N*S×N, final demand

# gross output from sector k and country i per country = exports per country (without final demand => intermediate demand): sum all country columns of row k, i
i_demand = [sum(D[i, j:j+S-1]) for i in 1:size(D, 1), j in 1:S:size(D, 2)] # N*S×N, intermediate demand per country

# gross exports from sector k and country i per country: sum intermediate and final demand
exports = i_demand + f_demand #N*S×N, exports per country

# gross final demand from sector k and country i
final = [sum(f_demand[i, :]) for i in 1:size(f_demand, 1)] # N*S×1

# gross total demand from sector k and country i
outputs_total = outputs .+ final

# value added of sector k and country i (total): (output-input)/output
v_add_share = (outputs_total .- inputs) ./ outputs_total # N*S×1, with final demand

all_value_added = hcat(raw[6:end,[:x4, :x1]], v_add_share, makeunique=true)
rename!(all_value_added, [:iso3, :WIOD, :value_added])
transform!(all_value_added, [:iso3, :WIOD] .=> ByRow(string) .=> [:iso3, :WIOD], renamecols=false) # change type to string to use XLSX.writetable()

XLSX.writetable("clean/all_value_added.xlsx", all_value_added, overwrite=true)

# ---------- Leontief inverse matrix -------------------------------------------------------------------------------------------------------

# technical coefficients: "a dollar's value of inputs from country i into sector k per a dollar's worth of output of country j in sector s"
A = similar(D)

# divide by entire output of country-sector pair
for i in 1:size(A, 1)
    for j in 1:size(A, 2)
        A[i, j] = D[i, j] / outputs_total[j]
    end
end

replace!(A, NaN => 0.0) # remove all NaN, mostly sectors U, T

# Leontief inverse matrix
# Leontief coeffiecients measure the total of dollars worth of country j, sector k goods required to meet 1 dollar woth of country i, sector z 's final demand
L = inv(I - A)


text:
using Documenter, Extremes, Cairo, Fontconfig

CI = get(ENV, "CI", nothing) == "true"

makedocs(sitename = "Extremes.jl",
    format = Documenter.HTML(
    prettyurls = CI,
    ),
    pages = [
       "index.md",
       "Tutorial" =>["Getting started" => "tutorial/index.md",
            "Block maxima" => "tutorial/BlockMaxima.md",
            "Threshold exceedance" => "tutorial/ThresholdExceedance.md",
            "Declustering" => "tutorial/Declustering.md",
            "Non-stationary block maxima" => "tutorial/blockmaxima_ns.md",
            "Non-stationary threshold exceedance" => "tutorial/thresholdexceedance_ns.md"],
       "contributing.md",
       "functions.md"
       ]

)

# if CI
#     deploydocs(
#     repo   = "github.com/jojal5/Extremes.jl.git",
#     devbranch = "dev",
#     versions = ["stable" => "v^", "v#.#"],
#     push_preview = false,
#     target = "build"
#     )
# end

deploydocs(repo = "github.com/jojal5/Extremes.jl.git")


text:
# Function used to create instance of powerSystemData type from MATLAB workspace:
function psDataLoad()
    (Sb, f, t, r, x, b, Y, bustype,
     Gp, Gq, Dp, Dq, Rp, Rq,
     Pmax, Pmin, Qmax, Qmin,
     Plim, Vg, Vceiling, Vfloor,
    busIdx, N, Nr, Ng) = readRTS96Data()

    # Allow each generator to participate equally in droop response.
    # Note: this only applies to analysis types with droop response!
    k = Float64[]
    for i = 1:length(Gp)
        if Gp[i] != 0
            push!(k, 1/length(find(Gp)))
        else
            push!(k,0)
        end
    end

    return powerSystemData(Sb,f,t,r,x,b,Y,bustype,Gp,Gq,Dp,Dq,Rp,Rq,Pmax,Pmin,Qmax,Qmin,Plim,Vg,Vceiling,Vfloor,busIdx,N,Nr,Ng,k)
end

function readRTS96Data();
    # Load RTS96 data
    # uses MAT
    path="../old/caseRTS96.mat" # Assumes current dir is nbs
    caseRTS96 = matread(path); # import MATLAB workspace

    # Connect relevant MATLAB variables to Julia
    bus_i = caseRTS96["bus_i"][:,1]; # vector of unique bus indices (73)
    bus = caseRTS96["bus"][:,1]; # vector of generator bus indices (99)
    Sb = caseRTS96["Sb"]; # base MVA
    bustype = caseRTS96["type"][:,1]

    Gp_long = caseRTS96["Pg"]; # conventional active power output, divide by Sb later
    Gq_long = caseRTS96["Qg"]; # conventional reactive power output, divide by Sb later

    Rp = caseRTS96["Wind"]./Sb; Rp = Rp[:,1]; # renewable active generation, has zero where no farm exists
    Rq = caseRTS96["Wind"].*(0.484/Sb); Rq = Rq[:,1];

    Dp = caseRTS96["Pd"]./Sb + Rp; # I previously subtracted wind from active demand
    Dp = Dp[:,1];
    Dq = caseRTS96["Qd"]./Sb; # reactive power demand vector
    Dq = Dq[:,1];

    Pmax = caseRTS96["Pmax"]./Sb; Pmax = Pmax[:,1]; # conventional gen. max active output
    Pmin = caseRTS96["Pmin"]./Sb; Pmin = Pmin[:,1]; # conventional gen. min active output

    Qmax = caseRTS96["Qmax"]./Sb; Qmax = Qmax[:,1]; # conventional gen. max reactive output
    Qmin = caseRTS96["Qmin"]./Sb; Qmin = Qmin[:,1]; # conventional gen. min reactive output

    Vm = caseRTS96["Vm"][:,1]; # should be vector of ones
    Va = caseRTS96["Va"][:,1]; # should be vector of zeros

    Plim = caseRTS96["rateA"]./Sb; Plim = Plim[:,1];
    Vceiling = caseRTS96["Vmax"][:,1]; # one voltage mag. ceiling constraint per node
    Vfloor = caseRTS96["Vmin"][:,1]; # one voltage mag. floor constraint per node

    Vceiling[bustype .== 2] = 1; # ceiling is 1 for voltage-controlled nodes
    Vfloor[bustype .== 2] = 1; # floor is 1 for voltage-controlled nodes

    Vg = caseRTS96["Vg"][:,1]; # voltage magnitudes from Jenny's solved ACOPF

    f = int(caseRTS96["fbus"][:,1]); # "from bus" ...
    t = int(caseRTS96["tbus"][:,1]); # ... "to bus"
    r = caseRTS96["r"][:,1]; # resistance, pu
    x = caseRTS96["x"][:,1]; # reactance, pu
    b = caseRTS96["b"][:,1];  # susceptance, pu

    # map bus numbers to 1:73
    for i = 1:length(f)
        if f[i] < 201
            f[i] -= 100
        elseif f[i] < 301
            f[i] -= 176;
        else
            f[i] -= 252;
        end
    end
    for i = 1:length(t)
        if t[i] < 201
            t[i] -= 100;
        elseif t[i] < 301
            t[i] -= 176;
        else
            t[i] -= 252;
        end
    end

    busidx = Int64[1];
    for i = 2:length(bus)
        if bus[i] < 201
            push!(busidx, bus[i]-100);
        elseif bus[i] < 301
            push!(busidx, bus[i]-176);
        else
            push!(busidx, bus[i]-252);
        end
    end

    # area 1: buses 1-24
    # area 2: buses 25-48
    # area 3: buses 49-73

    Gp = zeros(length(bus_i))
    Gq = zeros(length(bus_i))

    for i in unique(busidx)
        Gp[i] = sum(Gp_long.*[busidx .== i])/Sb;
        Gq[i] = sum(Gq_long.*[busidx .== i])/Sb;
    end
    # Now Gp and Gq reflect active and reactive generation at buses 1:73 consecutively.

    N = length(bus_i);
    Nr = length(find(Rp));
    Ng = length(find(Gp));

    Y = createY(f,t,r,x,b,true)
    # spy(Y)
    # title(L"Spy plot of $Y_{bus}$")

    return  Sb, f, t, r, x, b, Y, bustype,
            Gp, Gq, Dp, Dq, Rp, Rq,
            Pmax, Pmin, Qmax, Qmin,
            Plim, Vg, Vceiling, Vfloor,
            bus_i, N, Nr, Ng
end


text:
export Day6, Day6_1, Day6_2
abstract type Day6 <: Day end
abstract type Day6_1 <: Day6 end
abstract type Day6_2 <: Day6 end

solve(::Type{T}, lines::Vector) where T <: Day6 = begin
    lines |> split_on_empty_line |> lines -> day6_cond(T, lines)
end

day6_cond(::Type{Day6_1}, groups) = (groups .|> x-> join(x,"") |> Set |> length) |> sum
day6_cond(::Type{Day6_2}, groups) = (groups .|> group -> map(Set,group) |> x->intersect(x...) |> length) |> sum


text:
# This file was generated, do not modify it. # hide
X = MLJ.table(X)
y = categorical(y);


text:
"""
    GlobalFFTParams{T, N, inplace}

Specifies the global parameters for an N-dimensional distributed transform.
These include the element type `T` and global data sizes of input and output
data, as well as the transform types to be performed along each dimension.

---

    GlobalFFTParams(size_global, transforms, [real_type=Float64])

Define parameters for N-dimensional transform.

`transforms` must be a tuple of length `N` specifying the transforms to be
applied along each dimension. Each element must be a subtype of
[`Transforms.AbstractTransform`](@ref). For all the possible transforms, see
[`Transform types`](@ref Transforms).

The element type must be a real type accepted by FFTW, i.e. either `Float32` or
`Float64`.

Note that the transforms are applied one dimension at a time, with the leftmost
dimension first for forward transforms.

# Example

To perform a 3D FFT of real data, first a real-to-complex FFT must be applied
along the first dimension, followed by two complex-to-complex FFTs along the
other dimensions:

```jldoctest
julia> size_global = (64, 32, 128);  # size of real input data

julia> transforms = (Transforms.RFFT(), Transforms.FFT(), Transforms.FFT());

julia> fft_params = PencilFFTs.GlobalFFTParams(size_global, transforms)
Transforms: (RFFT, FFT, FFT)
Input type: Float64
Global dimensions: (64, 32, 128) -> (33, 32, 128)
```

"""
struct GlobalFFTParams{T, N, inplace, F <: AbstractTransformList{N}}
    # Transforms to be applied along each dimension.
    transforms :: F

    size_global_in  :: Dims{N}
    size_global_out :: Dims{N}

    function GlobalFFTParams(size_global::Dims{N},
                             transforms::AbstractTransformList{N},
                             ::Type{T}=Float64,
                            ) where {N, T <: FFTReal}
        F = typeof(transforms)
        size_global_out = length_output.(transforms, size_global)
        inplace = is_inplace(transforms...)
        if inplace === nothing
            throw(ArgumentError(
                "cannot combine in-place and out-of-place transforms: $(transforms)"))
        end
        new{T, N, inplace, F}(transforms, size_global, size_global_out)
    end
end

Base.ndims(::Type{<:GlobalFFTParams{T,N}}) where {T,N} = N
Base.ndims(g::GlobalFFTParams) = ndims(typeof(g))
Transforms.is_inplace(g::GlobalFFTParams{T,N,I}) where {T,N,I} = I

function Base.show(io::IO, g::GlobalFFTParams)
    print(io, "Transforms: ", g.transforms)
    print(io, "\nInput type: ", input_data_type(g))
    print(io, "\nGlobal dimensions: ",
          g.size_global_in, " -> ", g.size_global_out)
    nothing
end

# Determine input data type for multidimensional transform.
input_data_type(g::GlobalFFTParams{T}) where T =
    _input_data_type(T, g.transforms...) :: DataType

function _input_data_type(::Type{T}, transform::AbstractTransform,
                          next_transforms::Vararg{AbstractTransform}) where {T}
    Tin = eltype_input(transform, T) :: DataType
    if Tin === Nothing
        # This is the case if `transform` can take both real and complex data
        # (for instance for `NoTransform` or `R2R`).
        # We check the next transform type.
        return _input_data_type(T, next_transforms...)
    end
    Tin
end

# If all calls to `eltype_input` return Nothing, then we return the given real
# type. This will be the case for combinations of real-to-real transforms.
_input_data_type(::Type{T}) where {T} = T


text:
"""
ap_test_cross(re, factors; cov_method = "time iid", lag = 0, test_method = "Shanken")

CROSS-SECTIONAL ASSET PRICING TESTS

Asset pricing models predict a restriction on the error terms of the second-stage cross-sectional regressions.
The asset pricing model is valid if the error terms of cross-sectional regressions are jointly zero-mean.

Two-Pass Regression Estimation Procedure
- First-stage: time-series regressions of excess returns on factors to obtain estimates of betas for each assets
- Second-stage: cross-sectional regressions (without intercept) of time-series average returns on estiamted betas to obtain risk premium for each assets

INPUTS
`re': T x N matrix of excess returns, where T is the number of periods for each test assets and N is the number of test assets
`factors': T x K matrix of factors. Factors must be excess returns for time-series tests
`cov_method': method for calculating the variance-covariance matrix of the error terms:
            - time iid: error terms are iid across time but cross-sectionally correlated
            - iid: error terms are iiid both across time and assets
            - NW: Newey-West standard errors with m lag
            - HH: Hansen-Hodrick standard errors with m lags
`lag': number of lags used in the Newey-West or Hansen-Hodrick standard errors
`test_method': method of the time-series asset pricing tests:
            - OLS: OLS
            - GLS: GLS
            - Shanken: OLS with Shanken Correction (1992)
            - GLS Shanken: GLS with Shanken Correction (1992)
            - GMM: GMM exploiting the moment conditions from the first-stage time-series regressions and second-stage cross-sectional regressions. The model is overidentified
            - GLS GMM: GLS GMM exploiting the moment conditions from the first-stage time-series regressions and second-stage cross-sectional regressions. The model is overidentified

OUTPUTS
`Ere' : N x 1 vector of expected excess returns for N test assets
`lambda': K x 1 vector of risk premiums for K factors
`alpha': N x 1 vector of alphas for N test assets

`test_statistics': Wald-statistics of the asset pricing tests
`pvalue': pvalue of the test statistics. The test statistics has a Chi-square distribution.

"""

function ap_test_cross(re, factors, cov_method = "time iid", lag = 0, test_method = "Shanken")
    (T, K) = size(factors);
    N = size(re, 2);
    cov_f = cov(factors); # coviance matrix for the factors K x K

    X = [ones(T,1) factors];
    (~, beta, resid, cov_resid) = OLSFn(re, X, cov_method, lag); # First-stage time-series regressions

    Ere = mean(re, dims = 1)'; # N x 1 vector of sample average excess returns for each test assets

    # Methods of cross-sectional tests
    if test_method == "OLS" # OLS
        lambda = (beta' * beta) \ (beta' * Ere); # OLS second stage regressions - Risk premiums K x 1
        alpha = Ere- beta * lambda; # Alaphas N x 1
        Q = I - beta * ((beta' * beta) \ beta');
        cov_alpha = 1.0 ./ T * Q * cov_resid * Q;
        test_statistics = alpha' * (cov_alpha \ alpha);
        pvalue = 1.0 .- cdf.(Chisq(N - K), test_statistics);

    elseif test_method == "GLS" # GLS
        lambda = (beta' * (cov_resid \ beta)) \ beta' * (cov_resid \ Ere); # GLS second stage regressions - Risk premiums K x 1
        alpha = Ere- beta * lambda; # Alaphas N x 1
        test_statistics = T * alpha' * (cov_resid \ alpha);
        pvalue = 1.0 .- cdf.(Chisq(N - K), test_statistics);

    elseif test_method == "Shanken" # OLS Shanken Correction (1992)
        lambda = beta' * beta \ (beta' * Ere); # OLS second stage regressions - Risk premiums K x 1
        alpha = Ere- beta * lambda; # Alaphas N x 1
        Q = I - beta * ((beta' * beta) \ beta');
        cov_alpha = 1.0 ./ T * Q * cov_resid * Q .* (1.0 .+ lambda' * (cov_f \ lambda));
        test_statistics = alpha' * (cov_alpha \ alpha);
        pvalue = 1.0 .- cdf.(Chisq(N - K), test_statistics);

    elseif test_method == "GLS Shanken" # GLS Shanken Correction (1992)
        lambda = (beta' * (cov_resid \ beta)) \ beta' * (cov_resid \ Ere); # GLS second stage regressions - Risk premiums K x 1
        alpha = Ere - beta * lambda; # Alaphas N x 1
        cov_alpha = (1.0 ./ T) * (cov_resid - beta * (beta' * (cov_resid \ beta) \ beta')) .* (1.0 .+ lambda' * (cov_f \ lambda));
        test_statistics = T * (1.0 .+ lambda' * (cov_f \ lambda)) * alpha' * (cov_resid \ alpha);
        pvalue = 1.0 .- cdf.(Chisq(N - K), test_statistics);

    elseif test_method == "GMM" || test_method == "GLS GMM" # OLS GMM and GLS GMM
        lambda = beta' * beta \ (beta' * Ere); # OLS second stage regressions - Risk premiums K x 1
        alpha = Ere- beta * lambda; # Alaphas N x 1
        h = [ones(T, 1) factors];
        f = kron(ones(1, K+1), resid) .* kron(h, ones(1, N)); # Time-series moment conditions generalized to k factors
        f = [f re .- kron(ones(T, 1), (beta * lambda)')]; # Adding cross-sectional moment conditions gernerlized to k factors - p.225 in Cochrane (2000)
        d = - kron(h' * h / T, I(N)); # Jacobian of the time-seires GMM moment conditions
        d = [d zeros((1 + K) * N, K); kron([0 lambda'],I(N)) -beta]; # Adding the Jacobian of the cross-sectional GMM moment conditions

        # Weigthing matrix
        if test_method == "GMM"
            a = [I zeros((K + 1) * N, N); zeros(K, (K + 1) * N) beta']; # weighting matrix for GMM
        else
            a = [I zeros((K + 1) * N, N); zeros(K, (K + 1) * N) beta' * (cov_resid \ I)]; # weighting matrix for GLS GMM
        end

        invad = (a * d) \ I;
        premom = I(size(d * invad * a)[1]) - d * invad * a; #  for the sampling distributionof the moments - p.190 eqn (147) in Cochrane (2000)
        Ef2 = CovMoment(f, cov_method, lag); # spectral density matrix
        varb2 = (1.0 ./ T) * invad * a * Ef2 * a' * invad; # variance of the coefficients - p.190 eqn (146) in Cochrane (2000)
        varmom = (1.0 ./ T) * premom * Ef2 * premom'; # asymptotic covariance matrix of the moments - p.190 eqn (147) in Cochrane (2000)

        test_statistics = alpha' * pinv(varmom[N * (K + 1) + 1:end, N * (K + 1) + 1:end]) * alpha; # test the moment conditions that the the error terms(alphas) in the second-stage cross-sectional regressions are jointly zero
        pvalue = 1.0 .- cdf.(Chisq(N - K), test_statistics);
    end

    return(Ere, lambda, alpha, test_statistics, pvalue)
end


text:
func(x,y) = x*y
dfunc(x,y) = ForwardDiff.derivative(x->func(x,y),x)


text:
@testset "generate.jl" begin
    manager = ExperimentManager(db="choleskytest")
    init_manager(manager)
    @testset "gendecomposition" begin
        instance_name, oid = gendecomposition(manager, "./case9.m",
                                              Dict(:degree_max => Dict(:degree => 10)),
                                              Dict(:degree_max => Dict(:degree => 10)),
                                              3)
        @test isdecomposition(manager.decompositions, instance_name, oid)
        decomposition = getdecomposition(manager.decompositions, instance_name, oid)
        @test haskey(decomposition["options_src"], "degree_max")
        @test haskey(decomposition["options_src"]["degree_max"], "degree")
        @test decomposition["options_src"]["degree_max"]["degree"] == 10
        @test haskey(decomposition["options_dst"], "degree_max")
        @test haskey(decomposition["options_dst"]["degree_max"], "degree")
        @test decomposition["options_dst"]["degree_max"]["degree"] == 10
        @test decomposition["nb_added_edges"] == 3
    end
    @testset "gendecomposition_all" begin
        instances_names, oids = gendecomposition_all(manager,
                                                     Dict(:degree_min => Dict(:degree => 0)),
                                                     Dict(:degree_min => Dict(:degree => 0)),
                                                     2)
        for (instance_name, oid) in zip(instances_names, oids)
            @test isdecomposition(manager.decompositions, instance_name, oid)
            decomposition = getdecomposition(manager.decompositions, instance_name, oid)
            @test haskey(decomposition["options_src"], "degree_min")
            @test haskey(decomposition["options_src"]["degree_min"], "degree")
            @test decomposition["options_src"]["degree_min"]["degree"] == 0
            @test haskey(decomposition["options_dst"], "degree_min")
            @test haskey(decomposition["options_dst"]["degree_min"], "degree")
            @test decomposition["options_dst"]["degree_min"]["degree"] == 0
            @test decomposition["nb_added_edges"] == 2
        end
    end
    @testset "gendecomposition_all_json" begin
        instances_names, oids = gendecomposition_all_json(manager, "config.json")
        for (instance_name, oid) in zip(instances_names, oids)
            @test isdecomposition(manager.decompositions, instance_name, oid)
            decomposition = getdecomposition(manager.decompositions, instance_name, oid)
            @test haskey(decomposition["options_src"], "size_biggest_clique_max")
            @test haskey(decomposition["options_src"]["size_biggest_clique_max"], "size")
            @test decomposition["options_src"]["size_biggest_clique_max"]["size"] == 9
            @test haskey(decomposition["options_dst"], "size_biggest_clique_max")
            @test haskey(decomposition["options_dst"]["size_biggest_clique_max"], "size")
            @test decomposition["options_dst"]["size_biggest_clique_max"]["size"] == 9
            @test decomposition["nb_added_edges"] == 4
        end
    end
    @testset "gendecomposition_all_jsons" begin
        ids = gendecomposition_all_jsons(manager, "config_mult.json")
        for (instance_name, oid) in zip(ids["gen1"]...)
            @test isdecomposition(manager.decompositions, instance_name, oid)
            decomposition = getdecomposition(manager.decompositions, instance_name, oid)
            @test haskey(decomposition["options_src"], "degree")
            @test haskey(decomposition["options_src"]["degree"], "degrees")
            @test decomposition["options_src"]["degree"]["degrees"] == [1, 2, 3, 4, 5, 6, 7, 8, 9]
            @test haskey(decomposition["options_dst"], "degree")
            @test haskey(decomposition["options_dst"]["degree"], "degrees")
            @test decomposition["options_dst"]["degree"]["degrees"] == [1, 2, 3, 4, 5, 6, 7, 8, 9]
        end
        for (instance_name, oid) in zip(ids["gen2"]...)
            @test isdecomposition(manager.decompositions, instance_name, oid)
            decomposition = getdecomposition(manager.decompositions, instance_name, oid)
            @test haskey(decomposition["options_src"], "size_smallest_clique_max")
            @test haskey(decomposition["options_src"]["size_smallest_clique_max"], "size")
            @test decomposition["options_src"]["size_smallest_clique_max"]["size"] == 9
            @test haskey(decomposition["options_dst"], "size_smallest_clique_max")
            @test haskey(decomposition["options_dst"]["size_smallest_clique_max"], "size")
            @test decomposition["options_dst"]["size_smallest_clique_max"]["size"] == 9
            @test decomposition["nb_added_edges"] == 4
        end
    end
    @testset "genoriginal" begin
        instance_name, oid = genoriginal(manager, "case9.m")
        @test isdecomposition(manager.decompositions, instance_name, oid)
        decomposition = getdecomposition(manager.decompositions, instance_name, oid)
        @test isempty(decomposition["options_src"])
        @test isempty(decomposition["options_dst"])
        @test decomposition["nb_added_edges"] == 0
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using CSV
using MetaGraphs
using GraphDataFrameBridge
using LightGraphs

path = "examples/knowledge_graph/data/kg_vertex_types.csv"
edgepath = "examples/knowledge_graph/data/kg_edge_types.csv"
vf = CSV.read(path)
ef = CSV.read(edgepath)
g = GraphDataFrameBridge.metagraph_from_dataframe(MetaDiGraph, ef, :src_type, :dst_type)
set_indexing_prop!(g, :name)
println(g)
for v in vertices(g)
    name = props(g, v)[:name]
    set_prop!(g, v, :label, name)
end

for row in eachrow(ef)
    @show row
    @show v, u = g[row.src_type, :name], g[row.dst_type, :name]
    @show set_prop!(g, v,u, :label, row.edge_relation)
end

open("doc/src/img/olog.dot", "w") do fp
    MetaGraphs.savedot(fp, g)
end


text:
using SpecialFunctions

# https://github.com/JuliaMath/SpecialFunctions.jl/blob/d18ff04178dd37a60cc716a60ab3caf1a6903f43/src/gamma.jl#L20
@cufunc function SpecialFunctions.digamma(x::T) where {T<:Real}
    ψ = zero(x)
    if x ≤ 0 # reflection formula
        ψ -= π / tan(π * x)
        x = 1 - x
    end

    if x < 7
        # shift using recurrence formula
        ν = one(x)
        n = 7 - floor(x)
        while ν <= n - 1
            ψ -= inv(x + ν)
            ν += one(x)
        end
        ψ -= inv(x)
        x += n
    end
    t = inv(x)
    ψ += log(x) - t / 2
    t *= t # 1/z^2
    # the coefficients here are Float64(bernoulli[2:9] .// (2*(1:8)))
    ψ -= (
        t *
        @evalpoly(
            t,
            T(0.08333333333333333), T(-0.008333333333333333),
            T(0.003968253968253968), T(-0.004166666666666667),
            T(0.007575757575757576), T(-0.021092796092796094),
            T(0.08333333333333333), T(-0.4432598039215686)
        )
    )
    return ψ
end

function _cutrigamma(x::T) where {T<:Real}
    ψ = zero(x)
    if x < 8
        # shift using recurrence formula
        n = 8 - CUDA.floor(x)
        ψ += inv(x)^2
        ν = one(x)
        while ν <= n - 1
            ψ += inv(x + ν)^2
            ν += one(x)
        end
        x += n
    end
    t = inv(x)
    w = t * t # 1/z^2
    ψ += t + w / 2
    # the coefficients here are Float64(bernoulli[2:9])
    p = @evalpoly(
        w,
        T(0.16666666666666666), T(-0.03333333333333333),
        T(0.023809523809523808), T(-0.03333333333333333),
        T(0.07575757575757576), T(-0.2531135531135531),
        T(1.1666666666666667), T(-7.092156862745098)
    )
    ψ += t * w * p
    return ψ
end

@cufunc function SpecialFunctions.trigamma(x)
    if x ≤ 0 # reflection formula
        return (π / CUDA.sin(π * x))^2 - _cutrigamma(1 - x)
    else
        return _cutrigamma(x)
    end
end

@cufunc SpecialFunctions.loggamma(x) = CUDA.lgamma(x)

@cufunc_register SpecialFunctions.logbeta
@cufunc_function SpecialFunctions.logbeta(a, b) = (
      SpecialFunctions.loggamma(a)
    + SpecialFunctions.loggamma(b)
    - SpecialFunctions.loggamma(a + b)
)

@cufunc SpecialFunctions.erf(x) = CUDA.erf(x)
@cufunc SpecialFunctions.erfc(x) = CUDA.erfc(x)
@cufunc SpecialFunctions.erfcx(x) = CUDA.erfcx(x)
@cufunc SpecialFunctions.erfinv(x) = CUDA.erfinv(x)
@cufunc SpecialFunctions.erfcinv(x) = CUDA.erfcinv(x)

#####################
### Distributions ###
#####################
using StatsFuns
import StatsFuns: gammalogpdf, zval, xval, normlogpdf, normcdf, log2π, invsqrt2

@cufunc StatsFuns.gammalogpdf(k::Real, θ::Real, x::Number) = -loggamma(k) - k * log(θ) + (k - 1) * log(x) - x / θ

@cufunc function StatsFuns.xval(μ::Real, σ::Real, z::Number)
    if isinf(z) && iszero(σ)
        μ + one(σ) * z
    else
        μ + σ * z
    end
end
@cufunc StatsFuns.zval(μ::Real, σ::Real, x::Number) = (x - μ) / σ

# logpdf
@cufunc_register StatsFuns.normlogpdf
@cufunc_function StatsFuns.normlogpdf(z::Number) = -(abs2(z) + log2π)/2
@cufunc_function function StatsFuns.normlogpdf(μ::Real, σ::Real, x::Number)
    z = zval(μ, σ, x)
    normlogpdf(z) - log(σ)
end

# cdf
@cufunc_register StatsFuns.normcdf
@cufunc_function StatsFuns.normcdf(z::Number) = erfc(-z * invsqrt2)/2
@cufunc_function function StatsFuns.normcdf(μ::Real, σ::Real, x::Number)
    z = zval(μ, σ, x)
    normcdf(z)
end

# logcdf
@cufunc_register StatsFuns.normlogcdf
@cufunc_function function StatsFuns.normlogcdf(z::Number)
    invsqrt2 = inv(sqrt(oftype(z, 2)))
    if z < -1.0
        log(SpecialFunctions.erfcx(-z * invsqrt2)/2) - abs2(z)/2
    else
        SpecialFunctions.log1p(-SpecialFunctions.erfc(z * invsqrt2)/2)
    end
end
@cufunc_function function StatsFuns.normlogcdf(μ::Real, σ::Real, x::Number)
    z = zval(μ, σ, x)
    normlogcdf(z)
end

# truncated normlogpdf
@cufuncf function truncatednormlogpdf(μ, σ, x, lb, ub)
    logtp = StatsFuns.normlogcdf(μ, σ, ub) - StatsFuns.normlogcdf(μ, σ, lb)
    # TODO: deal with outside of boundary
    StatsFuns.normlogpdf(μ, σ, x) - logtp
    # TODO: seems like there's something messed up with the way we return `Inf`
    # if lb <= x <= ub
    #     StatsFuns.normlogpdf(μ, σ, x) - logtp
    # else
    #     TF = float(eltype(x))
    #     -TF(Inf)
    # end
end


text:
# Helper functions for pvt properties.

#%% Gas

"""
`LeeGasViscosity(specificGravity, psiAbs, tempF, Z)`

Gas viscosity (μ_g) in centipoise.

Takes gas specific gravity, psia, °F, Z (deviation factor).

Lee et al 1966 method.
"""
function LeeGasViscosity(specificGravity, psiAbs, tempF, Z)

  tempR = tempF + 459.67
  molecularWeight = 28.967 * specificGravity
  density = (psiAbs * molecularWeight * 0.00149406) / (Z*tempR)

  K = ((0.00094 + 2.0* 10^-6.0 *molecularWeight) * tempR^1.5 ) / (200.0 + 19.0*molecularWeight + tempR)
  X = 3.5 + (986.0/tempR) + 0.01*molecularWeight
  Y = 2.4 - 0.2*X

  return K * exp(X * density^Y ) #in centipoise
end

"""
`HankinsonWithWichertPseudoCriticalTemp(specificGravity, molFracCO2, molFracH2S)`

Pseudo-critical temperature, adjusted pseudo-critical temperature in °R, and Wichert correction factor

Takes gas specific gravity, mol fraction of CO₂, mol fraction of H₂S.

Hankin-Thomas-Phillips method, with Wichert and Aziz correction for sour components.
"""
function HankinsonWithWichertPseudoCriticalTemp(specificGravity, molFracCO2, molFracH2S)
  #Hankinson-Thomas-Phillips pseudo-parameter:
  tempPseudoCritical = 170.5 + (307.3*specificGravity)

  #Wichert and Aziz pseudo-param correction for sour gas components:
  A = molFracCO2 + molFracH2S
  B = molFracH2S
  correctionFactor = 120.0*(A^0.9 - A^1.6) + 15.0*(B^0.5 - B^4.0)

  tempPseudoCriticalMod = tempPseudoCritical - correctionFactor

  return tempPseudoCritical, tempPseudoCriticalMod, correctionFactor
end

"""
`HankinsonWithWichertPseudoCriticalPressure(specificGravity, molFracCO2, molFracH2S)`

Pseudo-critical pressure in psia.

Takes gas specific gravity, mol fraction of CO₂, mol fraction of H₂S.

Hankin-Thomas-Phillips method, with Wichert and Aziz correction for sour components.
"""
function HankinsonWithWichertPseudoCriticalPressure(specificGravity, molFracCO2, molFracH2S)
  #Hankinson-Thomas-Phillips pseudo-parameter:
  pressurePseudoCritical = 709.6 - (58.7*specificGravity) #I think the 58.7 is incorrectly identified as 56.7 in Takacs text

  #pseudo-temperatures:
  tempPseudoCritical, tempPseudoCriticalMod, correctionFactor = HankinsonWithWichertPseudoCriticalTemp(specificGravity, molFracCO2, molFracH2S)
  B = molFracH2S

  return pressurePseudoCritical * tempPseudoCriticalMod / (tempPseudoCritical + B*(1-B)*correctionFactor)
end

"""
`PapayZFactor(pressurePseudoCritical, tempPseudoCriticalRankine, psiAbs, tempF)`

Natural gas compressibility deviation factor (Z).

Take pseudocritical pressure (psia), pseudocritical temperature (°R), pressure (psia), temperature (°F).

Papay 1968 method.
"""
function PapayZFactor(pressurePseudoCritical, tempPseudoCriticalRankine, psiAbs, tempF)

  pressurePseudoCriticalReduced = psiAbs/pressurePseudoCritical

  tempPseudoCriticalReduced = (tempF + 459.67)/tempPseudoCriticalRankine

  return 1 - (3.52*pressurePseudoCriticalReduced)/(10^ (0.9813*tempPseudoCriticalReduced)) +
  (0.274*pressurePseudoCriticalReduced*pressurePseudoCriticalReduced)/(10^(0.8157*tempPseudoCriticalReduced))
end

"""
`KareemEtAlZFactor(pressurePseudoCritical, tempPseudoCriticalRankine, psiAbs, tempF)`

Natural gas compressibility deviation factor (Z).

Take pseudocritical pressure (psia), pseudocritical temperature (°R), pressure (psia), temperature (°F).

Can use gauge pressures so long as unit basis matches.

Direct correlation continuous over 0.2 ≤ P_pr ≤ 15.

Kareem, Iwalewa, Al-Marhoun, 2016.
"""
function KareemEtAlZFactor(pressurePseudoCritical, tempPseudoCriticalRankine, psiAbs, tempF)

  P_pr = psiAbs/pressurePseudoCritical

  T_pr = (tempF + 459.67)/tempPseudoCriticalRankine

  if !(0.2 ≤ P_pr ≤ 15) | !(1.15 ≤ T_pr ≤ 3)
    @info "Using Kareem et al Z-factor correlation with values outside of 0.2 ≤ P_pr ≤ 15, 1.15 ≤ T_pr ≤ 3."
  end

  a1 = 0.317842
  a2 = 0.382216
  a3 = -7.768354
  a4 = 14.290531
  a5 = 0.000002
  a6 = -0.004693
  a7 = 0.096254
  a8 = 0.16672
  a9 = 0.96691
  a10 = 0.063069
  a11 = -1.966847
  a12 = 21.0581
  a13 = -27.0246
  a14 = 16.23
  a15 = 207.783
  a16 = -488.161
  a17 = 176.29
  a18 = 1.88453
  a19 = 3.05921

  t = 1 / T_pr
  A = a1 * t * exp(a2 * (1-t)^2) * P_pr
  B = a3 * t + a4 * t^2 + a5 * t^6 * P_pr^6
  C = a9 + a8 * t * P_pr + a7 * t^2 * P_pr^2 + a6 * t^3 * P_pr^3
  D = a10 * t * exp(a11 * (1-t)^2)
  E = a12 * t + a13 * t^2 + a14 * t^3
  F = a15 * t + a16 * t^2 + a17 * t^3
  G = a18 + a19 * t

  y = D * P_pr / ((1 + A^2) / C - (A^2 * B) / C^3)
  z = D * P_pr * (1 + y + y^2 - y^3) / (D * P_pr + E * y^2 - F * y^G) / (1 - y)^3

  return z
end

"""
`KareemEtAlZFactor_simplified(pressurePseudoCritical, tempPseudoCriticalRankine, psiAbs, tempF)`

Natural gas compressibility deviation factor (Z).

Take pseudocritical pressure (psia), pseudocritical temperature (°R), pressure (psia), temperature (°F).

Linearized form from Kareem, Iwalewa, Al-Marhoun, 2016.
"""
function KareemEtAlZFactor_simplified(pressurePseudoCritical, tempPseudoCriticalRankine, psiAbs, tempF)

  P_pr = psiAbs/pressurePseudoCritical

  T_pr = (tempF + 459.67)/tempPseudoCriticalRankine

  if !(0.2 ≤ P_pr ≤ 15) | !(1.15 ≤ T_pr ≤ 3)
    @info "Using Kareem et al Z-factor correlation with values outside of 0.2 ≤ P_pr ≤ 15, 1.15 ≤ T_pr ≤ 3."
  end

  a1 = 0.317842
  a2 = 0.382216
  a3 = -7.768354
  a4 = 14.290531
  a5 = 0.000002
  a6 = -0.004693
  a7 = 0.096254
  a8 = 0.16672
  a9 = 0.96691

  t = 1 / T_pr
  A = a1 * t * exp(a2 * (1-t)^2) * P_pr
  B = a3 * t + a4 * t^2 + a5 * t^6 * P_pr^6
  C = a9 + a8 * t * P_pr + a7 * t^2 * P_pr^2 + a6 * t^3 * P_pr^3

  z = (1+ A^2)/C - A^2 * B / C^3

  return z
end

"""
`gasVolumeFactor(pressureAbs, Z, tempF)`

Corrected gas volume factor (B_g).

Takes absolute pressure (psia), Z-factor, temp (°F).
"""
function gasVolumeFactor(pressureAbs, Z, tempF)

  return 0.0283 * (Z*(tempF+459.67)/pressureAbs)
end

"""
`gasDensity_insitu(specificGravityGas, Z_factor, abspressure, tempF)`

In-situ gas density in lb/ft³ (ρ_g).

Takes gas s.g., Z-factor, absolute pressure (psia), temperature (°F).
"""
function gasDensity_insitu(specificGravityGas, Z_factor, abspressure, tempF)

  tempR = tempF + 459.67

  return 2.7 * specificGravityGas * abspressure / (Z_factor * tempR)
end

#%% Oil
#TODO: add Vasquez-Beggs for solution GOR
#TODO: add Hanafy et al correlations

"""
`StandingBubblePoint(APIoil, sg_gas, R_b, tempF)`

Bubble point pressure in psia.

Takes oil gravity (°API), gas specific gravity, total solution GOR at pressures above bubble point (R_b, scf/bbl), temp (°F).
"""
function StandingBubblePoint(APIoil, sg_gas, R_b, tempF)
    y = 0.00091 * tempF - 0.0125 * APIoil
    return 18.2 * ((R_b/sg_gas)^0.83 * 10^y - 1.4)
end

"""
`StandingSolutionGOR(APIoil, specificGravityGas, psiAbs, tempF)`

Solution GOR (Rₛ) in scf/bbl.

Takes oil gravity (°API), gas specific gravity, pressure (psia), temp (°F), total solution GOR (R_b, scf/bbl), and bubblepoint value (psia).

Standing method.
"""
function StandingSolutionGOR(APIoil, specificGravityGas, psiAbs, tempF, R_b, bubblepoint::Real)

    if psiAbs >= bubblepoint
        return R_b
    else
        y = 0.00091*tempF - 0.0125*APIoil

        return specificGravityGas * (psiAbs / (18 * 10^y) )^1.205 #scf/bbl
    end
end

"""
`StandingSolutionGOR(APIoil, specificGravityGas, psiAbs, tempF)`

Solution GOR (Rₛ) in scf/bbl.

Takes oil gravity (°API), gas specific gravity, pressure (psia), temp (°F), total solution GOR (R_b, scf/bbl), and bubblepoint function.

Standing method.
"""
function StandingSolutionGOR(APIoil, specificGravityGas, psiAbs, tempF, R_b, bubblepoint::Function)

    if psiAbs >= bubblepoint(APIoil, specificGravityGas, R_b, tempF)
        return R_b
    else
        y = 0.00091*tempF - 0.0125*APIoil

        return specificGravityGas * (psiAbs / (18 * 10^y) )^1.205 #scf/bbl
    end
end

"""
`StandingOilVolumeFactor(APIoil, specificGravityGas, solutionGOR, psiAbs, tempF)`

Oil volume factor (Bₒ).

Takes oil gravity (°API), gas specific gravity, solution GOR (scf/bbl), absolute pressure (psia), temp (°F).

Standing method.
"""
function StandingOilVolumeFactor(APIoil, specificGravityGas, solutionGOR, psiAbs, tempF)

  fFactor = solutionGOR * sqrt(specificGravityGas/(141.5/(APIoil + 131.5))) + 1.25*tempF

  return 0.972 + (0.000147) * fFactor^1.175
end

"""
`oilDensity_insitu(APIoil,  specificGravityGas,  solutionGOR,  oilVolumeFactor)`

Oil density (ρₒ) in mass-lbs per ft³.

Takes oil gravity (°API), gas specific gravity, solution GOR (R_s, scf/bbl), oil volume factor.
"""
function oilDensity_insitu(APIoil,  specificGravityGas,  solutionGOR,  oilVolumeFactor)

  return (141.5/(APIoil + 131.5)*350.4 + 0.0764*specificGravityGas*solutionGOR)/(5.61 * oilVolumeFactor) #mass-lbs per ft³
end

"""
`BeggsAndRobinsonDeadOilViscosity(APIoil, tempF)`

Dead oil viscosity (μ_oD) in centipoise.

Takes oil gravity (°API), temp (°F).

Use with caution at 100-150° F: viscosity can be significantly overstated.

Beggs and Robinson method.
"""
function BeggsAndRobinsonDeadOilViscosity(APIoil, tempF)

  if 100 <= tempF <= 155
    @info "Warning: using Beggs and Robinson for dead oil viscosity at $(tempF)° F--consider using another correlation for 100-150° F."
  end

  y = 10^(3.0324 - 0.02023*APIoil)
  x = y*tempF^-1.163

  return 10^x - 1
end

"""
`GlasoDeadOilViscosity(APIoil, tempF)`

Dead oil viscosity (μ_oD) in centipoise.

Takes oil gravity (°API), temp (°F).

Glaso method.
"""
function GlasoDeadOilViscosity(APIoil, tempF)

  return ((3.141* 10^10) / tempF^3.444) * log10(APIoil)^(10.313*log10(tempF) - 36.447)
end

"""
`ChewAndConnallySaturatedOilViscosity(deadOilViscosity,  solutionGOR)`

Saturated oil viscosity (μₒ) in centipoise.

Takes dead oil viscosity (cp), solution GOR (scf/bbl).

Chew and Connally method to correct from dead to live oil viscosity.
"""
function ChewAndConnallySaturatedOilViscosity(deadOilViscosity,  solutionGOR)

   A = 0.2 + 0.8 / 10.0^(0.00081*solutionGOR)
   b = 0.43 + 0.57 / 10.0^(0.00072*solutionGOR)

  return A * deadOilViscosity^b
end

#%% Water

"""
`waterDensity_stb(waterGravity)`

Water density in lb per ft³.

Takes water specific gravity.
"""
function waterDensity_stb(waterGravity)
  return waterGravity * 62.4 #lb per ft^3
end

"""
`waterDensity_insitu(waterGravity, B_w)`

Water density in lb per ft³.

Takes specific gravity, B_w.
"""
function waterDensity_insitu(waterGravity, B_w)
  return waterGravity * 62.4 / B_w #lb per ft^3
end

"""
`GouldWaterVolumeFactor(pressureAbs,  tempF)`

Water volume factor (B_w).

Takes absolute pressure (psia), temp (°F).

Gould method.
"""
function GouldWaterVolumeFactor(pressureAbs,  tempF)

  return 1.0 + 1.21 * 10^-4 * (tempF-60)  + 10^-6 * (tempF-60)^2 - 3.33 * pressureAbs * 10^-6
end

const assumedWaterViscosity = 1.0 #centipoise

#%% Interfacial tension

"""
`gas_oil_interfacialtension(APIoil, pressureAbsolute, tempF)`

Gas-oil interfactial tension in dynes/cm.

Takes oil API, absolute pressure (psia), temp (°F).

Possibly Baker Swerdloff method; same method utilized by [Fekete](https://www.ihsenergy.ca/support/documentation_ca/Harmony/content/html_files/reference_material/calculations_and_correlations/pressure_loss_calculations.htm)
"""
function gas_oil_interfacialtension(APIoil, pressureAbsolute, tempF)

  if tempF <= 68.0
    σ_dead = 39 - 0.2571 * APIoil
  elseif tempF >= 100.0
    σ_dead = 37.5 - 0.2571 * APIoil
  else
    σ_68 = 39 - 0.2571 * APIoil
    σ_100 = 37.5 - 0.2571 * APIoil
    σ_dead = σ_68 + (tempF - 68.0) * (σ_100 - σ_68) / (100.0 - 68.0) # interpolate
  end

  C = 1.0 - 0.024 * pressureAbsolute^0.45
  return C * σ_dead
end

"""
`gas_water_interfacialtension(pressureAbsolute, tempF)`

Gas-water interfactial tension in dynes/cm.

Takes absolute pressure (psia), temp (°F).

Possibly Baker Swerdloff method; same method utilized by [Fekete](https://www.ihsenergy.ca/support/documentation_ca/Harmony/content/html_files/reference_material/calculations_and_correlations/pressure_loss_calculations.htm)
"""
function gas_water_interfacialtension(pressureAbsolute, tempF)

  if tempF <= 74
    return 75 - 1.018 * pressureAbsolute^0.349
  elseif tempF >= 280
    return 53 - 0.1048 * pressureAbsolute^0.637
  else
    σ_74 = 75 - 1.018 * pressureAbsolute^0.349
    σ_280 = 53 - 0.1048 * pressureAbsolute^0.637
    return σ_74 + (tempF - 74.0) * (σ_280 - σ_74) / (280.0 - 74.0) #interpolate
  end
end