text:
"""
    sysinv = inv(sys; atol = 0, atol1 = atol, atol2 = atol, rtol, checkinv = true)

Compute for a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`, a descriptor realization of its inverse system
`sysinv = (Ai-λEi,Bi,Ci,Di)`, such that the transfer function matrix `Ginv(λ)` of `sysinv` is the inverse of `G(λ)` (i.e., `G(λ)*Ginv(λ) = I`).
The realization of `sysinv` is determined using inversion-free formulas and the invertibility condition is checked, unless `checkinv = false`.

The keyword arguments `atol1`, `atol2` and `rtol` specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `n*ϵ`, where `n` is the order of the square matrices `A` and `E`, and  `ϵ` is the working machine epsilon.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.
"""
function inv(sys::DescriptorStateSpace{T}; checkinv::Bool = true,
             atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
             rtol::Real = (size(sys.A,1)+1)*eps(real(float(one(T))))*iszero(min(atol1,atol2))) where T
   p, m = sys.ny, sys.nu
   n = sys.nx
   m == p  || error("The system must have the same number of inputs and outputs")
   Ai = [sys.A sys.B; sys.C sys.D]
   Bi = [zeros(T,n,m); -Matrix{T}(I,m,m)]
   Ci = [zeros(T,m,n) Matrix{T}(I,m,m)]
   Di = zeros(T,m,m)
   Ei = (sys.E == I) ? blockdiag(Matrix{T}(I,n,n),zeros(T,m,m)) :
                       blockdiag(sys.E,zeros(T,m,m))
   checkinv && (MatrixPencils.isregular(Ai, Ei, atol1 = atol1, atol2 = atol2, rtol = rtol ) ||
                error("The system is not invertible"))
   return DescriptorStateSpace{T}(Ai, Ei, Bi, Ci, Di, sys.Ts)
end
"""
    sysldiv = ldiv(sys1, sys2; atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ, checkinv = true)
    sysldiv = sys1 \\ sys2

Compute for the descriptor systems `sys1 = (A1-λE1,B1,C1,D1)` with the transfer function matrix `G1(λ)` and
`sys2 = (A2-λE2,B2,C2,D2)` with the transfer function matrix `G2(λ)`, a descriptor realization
`sysldiv = (Ai-λEi,Bi,Ci,Di)` of `sysldiv = inv(sys1)*sys2`,
whose transfer-function matrix `Gli(λ)` represents the result of the left division `Gli(λ) = inv(G1(λ))*G2(λ)`.
The realization of `sysldiv` is determined using inversion-free formulas and the invertibility condition for `sys1` is checked,
unless `checkinv = false`.

The keyword arguments `atol1`, `atol2` and `rtol` specify, respectively,
the absolute tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`,
the absolute tolerance for the nonzero elements of `E1` and `E2`,
and the relative tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`, `E1` and `E2`.
The default relative tolerance is `n*ϵ`, where `n` is the maximum of orders of the square matrices `A1` and `A2`, and  `ϵ` is the working machine epsilon.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.
"""
function ldiv(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2};
              atol::Real = zero(real(T1)), atol1::Real = atol, atol2::Real = atol,
              rtol::Real = (max(size(sys1.A,1),size(sys2.A,1))*eps(real(float(one(T1)))))*iszero(min(atol1,atol2))) where {T1,T2}
   p1, m1 = sys1.ny, sys1.nu
   m1 == p1  || error("The system SYS1 must be square")
   p2, m2 = sys2.ny, sys2.nu
   p1 == p2  || error("The systems SYS1 and SYS2 must have the same number of outputs")
   sys1.Ts == sys2.Ts ||  error("The systems SYS1 and SYS2 must have same sampling time")
   T = promote_type(T1,T2)
   n1 = sys1.nx
   n2 = sys2.nx
   if n1 > 0
      A, E, B, C, D = dssdata(T,[sys1 sys2])
      Ai = [A B[:,1:m1]; C D[:,1:m1]]
      Ei = E == I ? [I zeros(T,n1+n2,m1); zeros(T,m1,n1+n2+m1)]  : blockdiag(E,zeros(T,m1,m1))
      MatrixPencils.isregular(Ai, Ei, atol1 = atol1, atol2 = atol2, rtol = rtol) ||
                 error("The system SYS1 is not invertible")
      Bi = [B[:,m1+1:m1+m2]; D[:,m1+1:m1+m2]]
      Ci = [zeros(T,p1,n1+n2) -I]
      Di = zeros(T,p1,m2)
   else
      D1 = copy_oftype(sys1.D,T)
      LUD = lu(D1)
      (norm(D1,Inf) <= atol1 || rcond(LUD.U) <= 10*m1*eps(real(float(one(T))))) &&
               error("The system SYS2 is not invertible")
      Ai, Ei, Bi, Ci, Di = dssdata(T,sys2)
      ldiv!(LUD,Ci); ldiv!(LUD,Di)
   end

   return DescriptorStateSpace{T}(Ai, Ei, Bi, Ci, Di, sys2.Ts)
end
function (\)(sys1::DescriptorStateSpace, sys2::DescriptorStateSpace; kwargs...)
   ldiv(sys1,sys2; kwargs...)
end
"""
    sysrdiv = rdiv(sys1, sys2; atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ, checkinv = true)
    sysrdiv = sys1 / sys2

Compute for the descriptor systems `sys1 = (A1-λE1,B1,C1,D1)` with the transfer function matrix `G1(λ)` and
`sys2 = (A2-λE2,B2,C2,D2)` with the transfer function matrix `G2(λ)`, a descriptor realization
`sysrdiv = (Ai-λEi,Bi,Ci,Di)` of `sysrdiv = sys1*inv(sys2)`,
whose transfer-function matrix `Gri(λ)` represents the result of the right division `Gri(λ) = G1(λ)*inv(G2(λ))`.
The realization of `sysrdiv` is determined using inversion-free formulas and the invertibility condition for `sys2` is checked,
unless `checkinv = false`.

The keyword arguments `atol1`, `atol2` and `rtol` specify, respectively,
the absolute tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`,
the absolute tolerance for the nonzero elements of `E1` and `E2`,
and the relative tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`, `E1` and `E2`.
The default relative tolerance is `n*ϵ`, where `n` is the maximum of orders of the square matrices `A1` and `A2`,
and  `ϵ` is the working machine epsilon.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.
"""
function rdiv(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2};
              atol::Real = zero(real(T1)), atol1::Real = atol, atol2::Real = atol,
              rtol::Real = max(size(sys1.A,1),size(sys2.A,1))*eps(real(float(one(T1))))*iszero(min(atol1,atol2))) where {T1,T2}
   p2, m2 = sys2.ny, sys2.nu
   m2 == p2  || error("The system SYS2 must be square")
   p1, m1 = sys1.ny, sys1.nu
   m1 == m2  || error("The systems SYS1 and SYS2 must have the same number of inputs")
   sys1.Ts == sys2.Ts ||  error("The systems SYS1 and SYS2 must have same sampling time")
   T = promote_type(T1,T2)
   n1 = sys1.nx
   n2 = sys2.nx
   if n2 > 0
      A, E, B, C, D = dssdata(T,[sys2; sys1])
      Ai = [A B; C[1:p2,:] D[1:p2,:]]
      Ei = E == I ? [I zeros(T,n1+n2,p2); zeros(T,p2,n1+n2+p2)] : blockdiag(E,zeros(T,p2,p2))
      MatrixPencils.isregular(Ai, Ei, atol1 = atol1, atol2 = atol2, rtol = rtol) ||
                 error("The system SYS2 is not invertible")
      Ci = [C[p2+1:p1+p2,:] D[p2+1:p1+p2,:]]
      Bi = [zeros(T,n1+n2,m1); -I]
      Di = zeros(T,p1,m1)
   else
      D2 = copy_oftype(sys2.D,T)
      LUD = lu(D2)
      (norm(D2,Inf) <= atol1 || rcond(LUD.U) <= 10*m1*eps(real(float(one(T))))) &&
                  error("The system SYS2 is not invertible")
      Ai, Ei, Bi, Ci, Di = dssdata(T,sys1)
      rdiv!(Bi,LUD); rdiv!(Di,LUD)
   end

   return DescriptorStateSpace{T}(Ai, Ei, Bi, Ci, Di, sys1.Ts)
end
function (/)(sys1::DescriptorStateSpace, sys2::DescriptorStateSpace; kwargs...)
    rdiv(sys1,sys2; kwargs...)
end
"""
    sysdual = gdual(sys, rev = false)
    sysdual = transpose(sys, rev = false)

Compute for a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`,
the descriptor system realization of its dual system
`sysdual = (Ad-λEd,Bd,Cd,Dd)`, where `Ad = transpose(A)`, `Ed = transpose(E)`, `Bd = transpose(C)`,
`Cd = transpose(B)` and `Dd = transpose(D)`,
such that the transfer function matrix `Gdual(λ)` of `sysdual` is the transpose of `G(λ)`
(i.e., `Gdual(λ) = transpose(G(λ))`).

If `rev = true`, the tranposition is combined with the reverse permutation of the state variables, such that
`sysdual = (P*Ad*P-λP*Ed*P,P*Bd,Cd*P,Dd)`, where `P` is the permutation matrix with ones down the second diagonal.
"""
function gdual(SYS::DescriptorStateSpace{T};kwargs...) where T
    transpose(SYS;kwargs...)
end
function transpose(SYS::DescriptorStateSpace{T}; rev = false) where T
   if rev
       return DescriptorStateSpace{T}(reverse(reverse(transpose(SYS.A),dims=1),dims=2),
                                      SYS.E == I ? I : reverse(reverse(transpose(SYS.E),dims=1),dims=2),
                                      reverse(transpose(SYS.C),dims=1), reverse(transpose(SYS.B),dims=2), copy(transpose(SYS.D)), SYS.Ts)
   else
       return DescriptorStateSpace{T}(copy(transpose(SYS.A)), SYS.E == I ? I : copy(transpose(SYS.E)), copy(transpose(SYS.C)), copy(transpose(SYS.B)),
                                      copy(transpose(SYS.D)), SYS.Ts)
   end
end

"""
    sysconj = ctranspose(sys)
    sysconj = sys'

Compute the conjugate transpose (or adjoint) of a descriptor system (see [`adjoint`](@ref)).
"""
function ctranspose(sys::DescriptorStateSpace{T}) where T
   if sys.Ts == 0
      # continuous-time case
      return DescriptorStateSpace{T}(copy(-sys.A'), sys.E == I ? I : copy(sys.E'), copy(-sys.C'), copy(sys.B'), copy(sys.D'), sys.Ts)
   else
      # discrete-time case
      (nx, ny, nu) = (sys.nx, sys.ny, sys.nu)
      n = nx+ny
      return DescriptorStateSpace{T}(sys.E == I ? Matrix{T}(I,n,n) : blockdiag(sys.E',Matrix{T}(I,ny,ny)),
                                     [sys.A' sys.C'; zeros(T,ny,n)],
                                     [zeros(T,nx,ny) ; -I],
                                     [sys.B' zeros(T,nu,ny)], copy(sys.D'), sys.Ts)
   end
end
"""
    sysconj = adjoint(sys)
    sysconj = sys'

Compute for a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`,
the descriptor system realization of its adjoint (also called _conjugate transpose_) system
`sysconj = (Ac-λEc,Bc,Cc,Dc)`, such that the transfer function matrix `Gconj(λ)` of `sysconj`
is the appropriate conjugate transpose of `G(λ)`, as follows:
for a continuous-time system with `λ = s`, `Gconj(s) := transpose(G(-s))`, while
for a discrete-time system with `λ = z`, `Gconj(z) := transpose(G(1/z))`.
"""
function adjoint(sys::DescriptorStateSpace)
    ctranspose(sys)
end


text:
#--- set the seed or generate a random one
# const seed = Base.Random.RANDOM_SEED
const seed = Uint32[2657077506, 531413685, 532641760, 1678774672] # big bump
srand(seed)

# ---- set directory name for the simulation
const scriptname = "scriptTrack"
const savepath = joinpath(pwd(), "simulations", "$(scriptname)_$(seed[1])") #<--change the directory name here

# ---- parameters of the simulation run
const maxiter = 5_000  # sets the maximum number of gibbs iterations for each worker
const percentNyqForC = 0.5   # used for T l_max
const numofparsForP  = 1500  # used for P l_max
const hrfactor = 2.0
const pixel_size_arcmin = 2.0
const n = 2.0^9
const beamFWHM = 0.0
const nugget_at_each_pixel = (4.0)^2
begin
	local deltx =  pixel_size_arcmin * π / (180 * 60) # rads
	local period = deltx * n # side length in rads
	local deltk =  2 * π / period
	local nyq = (2 * π) / (2 * deltx)
	const maskupP  = √(deltk^2 * numofparsForP / π)  #l_max for for phi
	const maskupC  = min(9000.0, percentNyqForC * (2 * π) / (2 * pixel_size_arcmin * π / (180*60))) #l_max for for phi
	println("muK_per_arcmin = $(sqrt(nugget_at_each_pixel * (pixel_size_arcmin^2)))") # muK per arcmin
	println("maskupP = $maskupP") # muK per arcmin
	println("maskupC = $maskupC") # muK per arcmin
	println("$(seed[1])") # muK per arcmin
end
const scale_grad =  2.0e-3
const scale_hmc  =  2.0e-3

# ------------ load modules and functions
push!(LOAD_PATH, pwd()*"/src")
using Interp
FFTW.set_num_threads(15)
BLAS.blas_set_num_threads(15)
require("cmb.jl")
require("fft.jl")
require("funcs.jl") # use reload after editing funcs.jl
require("Interp.jl")

# --------- generate cmb spectrum class for high res and low res
parlr = setpar(
	pixel_size_arcmin, n, beamFWHM, nugget_at_each_pixel,
	maskupC, maskupP
)
parhr = setpar(
	pixel_size_arcmin./hrfactor, hrfactor*n, beamFWHM, nugget_at_each_pixel,
	maskupC, maskupP
)

# -------- Simulate data: ytx, maskvarx, phix, tildetx
ytk_nomask, tildetk, phix, tx_hr = simulate_start(parlr);
phik = fft2(phix, parlr)
tmpdo = maximum(parlr.grd.x)*0.3
tmpup = maximum(parlr.grd.x)*0.4
maskboolx = tmpdo .<= parlr.grd.x .<= tmpup
maskvarx  = parlr.nugget_at_each_pixel .* ones(size(parlr.grd.x))
maskvarx[maskboolx] = Inf
ytx = ifft2r(ytk_nomask, parlr)
ytx[maskboolx] = 0.0
ytk = fft2(ytx, parlr)

# -----------  set the save directory and save stuff
isdir(savepath) && run(`rm -rf $savepath`)
run(`mkdir $savepath`)
run(`cp src/funcs.jl $savepath/funcs_save.jl`)
run(`cp scripts/$scriptname.jl $savepath/$(scriptname)_save.jl`)
writecsv("$savepath/seed.csv", seed)
writecsv("$savepath/x.csv", parlr.grd.x)
writecsv("$savepath/y.csv", parlr.grd.y)
writecsv("$savepath/k1.csv", parlr.grd.k1)
writecsv("$savepath/k2.csv", parlr.grd.k2)
writecsv("$savepath/ytx.csv", ytx)
writecsv("$savepath/phix.csv", phix)
writecsv("$savepath/tildetx_lr.csv", ifft2r(tildetk, parlr))
writecsv("$savepath/tx_lr.csv", tx_hr[1:3:end, 1:3:end] )
writecsv("$savepath/maskvarx.csv", maskvarx)

# ------------------ initalized and run the gibbs
function gibbsloop(its, parhr, parlr, ytx, maskvarx, maskupC, savepath, scale_grad, scale_hmc, hrfactor, record_run)
	acceptclk   = [1] #initialize acceptance record
	tx_hr_curr  = zero(parhr.grd.x)
	ttx         = zero(parhr.grd.x)
	p1hr, p2hr  = zero(parhr.grd.x), zero(parhr.grd.x)
	phik_curr   = zero(fft2(ytx, parlr))
	phix_curr   = zero(ytx)
	tx_lr_curr  = zero(ytx)
	tildetx_hr_curr = zero(parhr.grd.x)

	for bglp = 1:its
		tic()
		# ----- update tildetx_hr_curr
		if bglp % 100 == 1
			p1hr[:], p2hr[:] = gibbspass_t!(tx_hr_curr, ttx, phik_curr, ytx,
				maskvarx, parlr, parhr, repstretch(2.0, 1.2maskupC, 75, 1000)
			)
		end
		p1hr[:], p2hr[:] = gibbspass_t!(tx_hr_curr, ttx, phik_curr, ytx,
			maskvarx, parlr, parhr, fill(Inf, 400)
		)

		tildetx_hr_curr[:] = spline_interp2(
			parhr.grd.x, parhr.grd.y, tx_hr_curr,
			parhr.grd.x + p1hr, parhr.grd.y + p2hr
		)

		# ----- update phik_curr
		push!(acceptclk, hmc!(phik_curr, tildetx_hr_curr, parlr, parhr, scale_hmc))

		# ----- now save
		phix_curr[:]  = ifft2r(phik_curr, parlr)
		tx_lr_curr[:] =  tx_hr_curr[1:int(hrfactor):end,1:int(hrfactor):end]
		record_run = hcat(record_run,
			[tx_lr_curr[400,300],
			 tx_lr_curr[50,50],
			 phix_curr[400,300],
			 phix_curr[50,50],
			 phik_curr[10,10],
			 phik_curr[4,7]])
		if bglp % 50 == 1
			writecsv("$savepath/trace.csv", record_run)
			writecsv("$savepath/phix_curr_$bglp.csv", phix_curr)
			writecsv("$savepath/acceptclk.csv", acceptclk)
		end # if
		toc()
	end # for
end # function

tx_lr = tx_hr[1:3:end,1:3:end]
record_run = [tx_lr[400,300],
			 tx_lr[50,50],
			 phix[400,300],
			 phix[50,50],
			 phik[10, 10],
			 phik[4, 7]]
gibbsloop(maxiter, parhr, parlr, ytx, maskvarx, maskupC, savepath, scale_grad, scale_hmc, hrfactor, record_run)

ext: jl

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
const tol = 10*eps()

"""
A struct detailing an ito integral. It contains a UnivariateFunction detailing the integrand as well as a symbol detailing an id of the integral's processes.

Usual (and most general) contructor is:

    ItoIntegral(brownian_id_::Symbol, f_::UnivariateFunction)

Convenience constructor for `ItoIntegral`s where the integrand is a flat function is:

    ItoIntegral(brownian_id_::Symbol, variance_::Real)
"""
struct ItoIntegral
    brownian_id_::Symbol
    f_::UnivariateFunction
    function ItoIntegral(brownian_id_::Symbol, variance_::Real)
        return new(brownian_id_, PE_Function(variance_, 0.0, 0.0, 0))
    end
    function ItoIntegral(brownian_id_::Symbol, f_::UnivariateFunction)
        return new(brownian_id_, f_)
    end
end

"""
    variance(ito::ItoIntegral, from::Real, to::Real)
    variance(ito::ItoIntegral, base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime})
Get the variance of an `ItoIntegral` from one point of time to another.
### Inputs
* `ito` - The `ItoIntegral` you want the variance for.
* `from` - The time at which the integration starts.
* `to` - The time at which the integration ends.
### Outputs
* A scalar
"""
function variance(ito::ItoIntegral, from::Real, to::Real)
    if from + eps() > to
        return 0.0;
    end
    return evaluate_integral(ito.f_^2, from, to)
end
function variance(ito::ItoIntegral, base::Union{Date,DateTime},
                           from::Union{Date,DateTime}, to::Union{Date,DateTime})
    from_fl = years_between(from, base)
    to_fl   = years_between(to, base)
    return variance(ito, from_fl, to_fl)
end

"""
    volatility(ito::ItoIntegral, on::Union{Date,DateTime})
Get the volatility of an `ItoIntegral` on a certain date.
### Inputs
* `ito` - The `ItoIntegral` you want the volatility for.
* `on` - What instant do you want the volatility for.
### Outputs
* A scalar
"""
function volatility(ito::ItoIntegral, on::Union{Date,DateTime})
    return ito.f_(on)
end
function volatility(ito::ItoIntegral, on::Real)
    return ito.f_(on)
end

"""
    covariance(ito1::ItoIntegral,ito2::ItoIntegral, from::Real, to::Real, gaussian_correlation::Real)
    covariance(ito1::ItoIntegral,ito2::ItoIntegral, base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime}, gaussian_correlation::Real)
Get the covariance of two `ItoIntegral`s over a certain period given the underlying Brownian processes have a correlation of `gaussian_correlation`.
### Inputs
* `ito1` - The first `ItoIntegral`
* `ito2` - The second `ItoIntegral`
* `from` - The start of the period
* `to` - The end of the period
* `gaussian_correlation` - The correlation between the brownians for each of the two itos. This should be in the range [-1,1].
* `on` - What instant do you want the volatility for.
### Outputs
* A scalar
"""
function covariance(ito1::ItoIntegral,ito2::ItoIntegral, from::Real, to::Real, gaussian_correlation::Real)
    if from + eps() >= to
        return 0.0;
    end
    return gaussian_correlation * evaluate_integral(ito1.f_ * ito2.f_, from, to)
end
function covariance(ito1::ItoIntegral,ito2::ItoIntegral, base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime}, gaussian_correlation::Real)
    from_fl = years_between(from, base)
    to_fl   = years_between(to, base)
    return covariance(ito1, ito2, from_fl, to_fl, gaussian_correlation)
end

"""
    brownians_in_use(itos::Array{ItoIntegral,1}, brownians::Array{Symbol,1})
Determine which Browninan processes are used in an array of `ItoIntegral`s.
### Inputs
* itos - A dict containing each of the ito integrals
* brownians - All possible brownians
### Outputs
* A `Vector` of what brownians are in use in itos
* A `Vector` with the indices of these brownians.
"""
function brownians_in_use(itos::Dict{Symbol,ItoIntegral}, brownians::Array{Symbol,1})
    all_brownians_in_use = unique(map(x -> x.brownian_id_ , values(itos)))
    indices_in_use   = unique(findall(map(x -> x in all_brownians_in_use , brownians)))
    reduced_brownian_list = brownians[indices_in_use]
    return all_brownians_in_use, indices_in_use, reduced_brownian_list
end

"""
Creates an ItoSet. This contains :
* A correlation matrix of brownian motions.
* A vector giving the axis labels for this correlation matrix.
* A dict of `ItoInteral`s. Here the keys should be ids for the ito integrals and the values should be `ItoIntegral`s.
Determine which Brownian processes are used in an array of `ItoIntegral`s.
"""
struct ItoSet{T<:Real}
    brownian_correlation_matrix_::Hermitian{T}
    brownian_ids_::Array{Symbol,1}
    ito_integrals_::Dict{Symbol,ItoIntegral}
    function ItoSet(brownian_corr_matrix::Hermitian{T}, brownian_ids::Array{Symbol,1},
                    ito_integrals::Dict{Symbol,ItoIntegral}) where T<:Real
        if size(brownian_ids)[1] != size(brownian_corr_matrix)[1]
            error("The shape of brownian_ids_ must match the number of rows/columns of brownian_correlation_matrix_")
        end
        all_brownians_in_use, used_brownian_indices, brown_ids = brownians_in_use(ito_integrals, brownian_ids)
        if length(setdiff(all_brownians_in_use, brownian_ids)) > 0
            error(string("In creating an ItoSet there are some brownian motions referenced by ito integrals for which there are no corresponding entries in the correlation matrix for brownian motions. Thus an ItoSet cannot be built. These include ", setdiff(all_brownians_in_use, brownian_ids)))
        end
        brownian_corr_matrix_subset = Hermitian(brownian_corr_matrix[used_brownian_indices,used_brownian_indices])
        return new{T}(brownian_corr_matrix_subset, brown_ids, ito_integrals)
    end
end

"""
    correlation(ito::ItoSet, index1::Integer, index2::Integer)
    correlation(ito::ItoSet, brownian_id1::Symbol, brownian_id2::Symbol)
Get correlation between brownian motions in an `ItoSet`.
### Inputs
* `ito` - An `ItoSet` that you want the correlation for two itos within.
* `index1` or `brownian_id1` - The index/key for the first ito integral.
* `index2` or `brownian_id2` - The index/key for the second ito integral.
### Returns
* A scalar
"""
function correlation(ito::ItoSet, index1::Integer, index2::Integer)
    return ito.brownian_correlation_matrix_[index1, index2]
end
function correlation(ito::ItoSet, brownian_id1::Symbol, brownian_id2::Symbol)
    index1 = findall(brownian_id1 .== ito.brownian_ids_)[1]
    index2 = findall(brownian_id2 .== ito.brownian_ids_)[1]
    return correlation(ito, index1, index2)
end

"""
    volatility(ito::ItoSet, ito_integral_id::Symbol, on::Union{Date,DateTime})
Get volatility of an `ito_integral` on a date.
### Inputs
* `ito` - An `ItoSet` that you want the volatility for.
* `ito_integral_id` - The key of the ito dict that you are interested in
* `on` The time or instant you want the volatility for.
### Returns
* A scalar
"""
function volatility(ito::ItoSet, ito_integral_id::Symbol, on::Union{Date,DateTime})
    ito_integral = ito.ito_integrals_[ito_integral_id]
    return volatility(ito_integral, on)
end

"""
    make_covariance_matrix(ito_set_::ItoSet, from::Real, to::Real)
Make a covariance matrix given an ItoSet and a period of time. This returns a
Hermitian covariance matrix as well as a vector of symbols representing the axis
labelling on this `Hermitian`.

### Inputs
* `ito_set_` - An `ItoSet` you want to make a covariance matrix from.
* `from` - The (numeric) time from which the covariance span starts.
* `to` - The (numeric) time at which the covariance span ends.
### Returns
* A `Hermitian` covariance matrix.
* A `Vector{Symbol}` of labels for the covariance matrix.
"""
function make_covariance_matrix(ito_set_::ItoSet{T}, from::Real, to::Real) where T<:Real
    number_of_itos = length(ito_set_.ito_integrals_)
    ito_ids = collect(keys(ito_set_.ito_integrals_))
    cov = Array{T,2}(undef, number_of_itos,number_of_itos)
    for r in 1:number_of_itos
        rito = ito_set_.ito_integrals_[ito_ids[r]]
        for c in r:number_of_itos
            cito = ito_set_.ito_integrals_[ito_ids[c]]
            cr_correlation = correlation(ito_set_, rito.brownian_id_, cito.brownian_id_)
            cov[r,c] = covariance(rito, cito, from, to, cr_correlation)
        end
    end
    return Hermitian(cov), ito_ids
end

"""
    StochasticIntegralsCovariance
This is an abstract type which represents structs that represent covariances over spans of time.
The concrete instances of this type should support extracting correlations, covariances, volatilities
and random number generation.
"""
abstract type StochasticIntegralsCovariance end

"""
    ForwardCovariance

Creates an `ForwardCovariance` struct. This contains :

* An `Itoset`.
* The time from which the covariance period starts.
* The time to which the covariance period extends to.
And in the constructor the following items are generated and stored in the object:
* A covariance matrix
* Labels for the covariance matrix.
* The cholesky decomposition of the covariance matrix.
* The inverse of the covariance matrix.
* The determinant of the covariance matrix.

The constructors are:

    ForwardCovariance(ito_set_::ItoSet, from_::Real, to_::Real;
                      calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
    ForwardCovariance(ito_set_::ItoSet, from::Union{Date,DateTime}, to::Union{Date,DateTime})
    ForwardCovariance(old_ForwardCovariance::ForwardCovariance, from::Real, to::Real)
    ForwardCovariance(old_ForwardCovariance::ForwardCovariance, from::Union{Date,DateTime}, to::Union{Date,DateTime})
"""
struct ForwardCovariance <:StochasticIntegralsCovariance
    ito_set_::ItoSet
    from_::Real
    to_::Real
    covariance_::Hermitian
    covariance_labels_::Array{Symbol,1}
    chol_::Union{Missing,LowerTriangular}
    inverse_::Union{Missing,Hermitian}
    determinant_::Union{Missing,Real}
    function ForwardCovariance(ito_set_::ItoSet, from_::Real, to_::Real;
             calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
        covariance_, covariance_labels_ = make_covariance_matrix(ito_set_, from_, to_)
        chol_ = missing
        inverse_ = missing
        determinant_ = missing
        if calculate_chol chol_ = LowerTriangular(cholesky(covariance_).L) end
        if calculate_inverse inverse_           = Hermitian(inv(covariance_)) end
        if calculate_determinant determinant_       = det(covariance_) end
        return new(ito_set_, from_, to_, covariance_, covariance_labels_, chol_, inverse_, determinant_)
    end
    function ForwardCovariance(ito_set_::ItoSet, from::Union{Date,DateTime}, to::Union{Date,DateTime};
             calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
        from_ = years_from_global_base(from)
        to_   = years_from_global_base(to)
        return ForwardCovariance(ito_set_, from_, to_; calculate_chol = calculate_chol,
                  calculate_inverse = calculate_inverse, calculate_determinant = calculate_determinant)
    end
    function ForwardCovariance(old_ForwardCovariance::ForwardCovariance, from::Real, to::Real)
        return ForwardCovariance(old_ForwardCovariance.ito_set_, from, to)
    end
    function ForwardCovariance(old_ForwardCovariance::ForwardCovariance, from::Union{Date,DateTime}, to::Union{Date,DateTime})
        from_ = years_from_global_base(from)
        to_   = years_from_global_base(to)
        return ForwardCovariance(old_ForwardCovariance.ito_set_, from_, to_)
    end
end

"""

    SimpleCovariance

Creates an `SimpleCovariance` struct. This is a simplified version of `ForwardCovariance`
and is designed for `ItoIntegral`s that are constant (so correlations do not need to be recalculated).
This computational saving is the only advantage - `ForwardCovariance` is more general.
It contains the same elements as a `ForwardCovariance`:
* An Itoset
* The time from which the covariance period starts.
* The time to which the covariance period extends to.
And in the constructor the following items are generated and stored in the object:
* A covariance matrix
* Labels for the covariance matrix.
* The cholesky decomposition of the covariance matrix.
* The inverse of the covariance matrix.
* The determinant of the covariance matrix.

The constructors are:

    SimpleCovariance(ito_set_::ItoSet, from_::Real, to_::Real;
         calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
"""
mutable struct SimpleCovariance <: StochasticIntegralsCovariance
    ito_set_::ItoSet
    from_::Real
    to_::Real
    covariance_::Hermitian
    covariance_labels_::Array{Symbol,1}
    chol_::Union{Missing,LowerTriangular}
    inverse_::Union{Missing,Hermitian}
    determinant_::Union{Missing,Real}
    function SimpleCovariance(ito_set_::ItoSet, from_::Real, to_::Real;
             calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
        fc = ForwardCovariance(ito_set_, from_, to_; calculate_chol = calculate_chol, calculate_inverse = calculate_inverse, calculate_determinant = calculate_determinant)
        return new(fc.ito_set_, fc.from_, fc.to_, fc.covariance_, fc.covariance_labels_, fc.chol_, fc.inverse_, fc.determinant_)
    end
end

"""
    update!(sc::SimpleCovariance, from::Real, to::Real)

This takes a `SimpleCovariance` and updates it for a new  span in time.
The new span in time is between from and to. For `SimpleCovariance` this is done by
just adjusting the covariances for the new time span (with corresponding adjustments)
to the cholesky, inverse, etc.

The corresponding technique for a `ForwardCovariance` (which is also a `StochasticIntegralsCovariance`)
is to feed it into a new `ForwardCovariance` constructor which will recalculate for the new span.
### Inputs
* `sc` - The `SimpleCovariance` struct.
* `from` - The time from which you want the covariance for.
* `to` - The time to which you want the covariance for.
### Returns
Nothing. It juts updates the sc struct you pass in as an input.
"""
function update!(sc::SimpleCovariance, from::Real, to::Real)
    old_duration = sc.to_ - sc.from_
    new_duration = to - from
    relative_duration = new_duration/old_duration
    sc.covariance_  = sc.covariance_ * (relative_duration)
    sc.chol_ = ismissing(sc.chol_) ? missing : LowerTriangular(sc.chol_ .* sqrt(relative_duration))
    sc.inverse_ = ismissing(sc.inverse_) ? missing : Hermitian(sc.inverse_ ./  (relative_duration))
    sc.determinant_ = ismissing(sc.determinant_) ? missing : sc.determinant_ *  ((relative_duration)^length(sc.covariance_labels_))
    sc.from_ = from
    sc.to_ = to
end

"""
    volatility(covar::ForwardCovariance, index::Integer, on::Union{Date,DateTime})
    volatility(covar::ForwardCovariance, id::Symbol, on::Union{Date,DateTime})
Get the volatility of an `ForwardCovariance` on a date.
### Inputs
* `covar` - An `ForwardCovariance` that you want the volatility for.
* `index` - The key of the ito dict that you are interested in
* `on` The time or instant you want the volatility for.
### Returns
* A scalar
"""
function volatility(covar::ForwardCovariance, index::Integer, on::Union{Date,DateTime})
    return volatility(covar.ito_set_, index, on)
end
function volatility(covar::ForwardCovariance, id::Symbol, on::Union{Date,DateTime})
    return volatility(covar.ito_set_, id, on)
end

"""
    variance(covar::ForwardCovariance, id::Symbol)
    variance(covar::ForwardCovariance, index::Integer)
Get the variance of an `ForwardCovariance` over a period.
### Inputs
* `covar` - An `ForwardCovariance` that you want the variance for.
* `id` or `index` - The key/index of the ito dict that you are interested in
### Returns
* A scalar
"""
function variance(covar::ForwardCovariance, id::Symbol)
        index = findall(id .== covar.covariance_labels_)[1]
        return variance(covar, index)
end
function variance(covar::ForwardCovariance, index::Integer)
    return covar.covariance_[index,index]
end

"""
    covariance(covar::ForwardCovariance, index_1::Integer, index_2::Integer)
    covariance(covar::ForwardCovariance, id1::Symbol, id2::Symbol)
Get the covariance of two ito integrals in a `ForwardCovariance` over a period.
### Inputs
* `covar` - An `ForwardCovariance` that you want the covariance for.
* `index_1` or `id1` - The key/index of the first ito that you are interested in
* `index_2` or `id2` - The key/index of the second ito that you are interested in
### Returns
* A scalar
"""
function covariance(covar::ForwardCovariance, index_1::Integer, index_2::Integer)
    return covar.covariance_[index_1,index_2]
end
function covariance(covar::ForwardCovariance, id1::Symbol, id2::Symbol)
    index_1 = findall(id1 .== covar.covariance_labels_)[1]
    index_2 = findall(id2 .== covar.covariance_labels_)[1]
    return covariance(covar, index_1, index_2)
end

"""
    correlation(covar::ForwardCovariance, index_1::Integer, index_2::Integer)
    correlation(covar::ForwardCovariance, id1::Symbol, id2::Symbol)
Get the correlation of two `ItoIntegral`s over a period.
### Inputs
* `covar` - An `ForwardCovariance` that you want the correlation for.
* `index_1` or `id1` - The key/index of the first ito that you are interested in
* `index_2` or `id2` - The key/index of the second ito that you are interested in
### Returns
* A scalar
"""
function correlation(covar::ForwardCovariance, index_1::Integer, index_2::Integer)
    cova = covariance(covar, index_1, index_2)
    var1  = variance(covar, index_1)
    var2  = variance(covar, index_2)
    return cova/sqrt(var1 * var2)
end
function correlation(covar::ForwardCovariance, id1::Symbol, id2::Symbol)
    index_1 = findall(id1 .== covar.covariance_labels_)[1]
    index_2 = findall(id2 .== covar.covariance_labels_)[1]
    return correlation(covar, index_1, index_2)
end

## Random draws
"""
    get_draws(covar::ForwardCovariance; uniform_draw::Array{T,1} = rand(length(covar.covariance_labels_))) where T<:Real

### Inputs
* `covar` - An `ForwardCovariance` or `SimpleCovariance` struct that you want to draw from.
* `uniform_draw`- The draw vector (from the uniform distribution)
### Returns
* A Dict of draws.
"""
function get_draws(covar::Union{ForwardCovariance,SimpleCovariance}; uniform_draw::Array{T,1} = rand(length(covar.covariance_labels_))) where T<:Real
    number_of_itos = length(covar.covariance_labels_)
    normal_draw = quantile.(Ref(Normal()), uniform_draw)
    scaled_draw = covar.chol_ * normal_draw
    first_set_of_draws = Dict{Symbol,Real}(covar.covariance_labels_ .=> scaled_draw)
    return first_set_of_draws
end
"""
    get_draws(covar::ForwardCovariance, num::Integer; twister::MersenneTwister = MersenneTwister(1234), antithetic_variates = false)
get pseudorandom draws from a `ForwardCovariance` struct. Other schemes (like quasirandom) can be done by inserting quasirandom
numbers in as the uniform_draw.
If the `antithetic_variates` control is set to true then every second set of draws will be antithetic to the previous.

### Inputs
* `covar` - An `ForwardCovariance` or `SimpleCovariance` struct that you want to draw from.
* `num`- The number of draws you want
* `twister`- The number of draws you want
### Returns
* A `Vector` of `Dict`s of draws.
"""
function get_draws(covar::Union{ForwardCovariance,SimpleCovariance}, num::Integer; number_generator::NumberGenerator = Mersenne(MersenneTwister(1234), length(covar.covariance_labels_)), antithetic_variates = false)
    if antithetic_variates
        half_num = convert(Int, round(num/2))
        array_of_dicts = Array{Dict{Symbol,Real}}(undef, half_num*2)
        number_of_itos = length(covar.covariance_labels_)
        for i in 1:half_num
            first_entry = 2*(i-1)+1
            draw = next!(number_generator)
            array_of_dicts[first_entry] = get_draws(covar; uniform_draw = draw)
            array_of_dicts[first_entry + 1] = get_draws(covar; uniform_draw = 1.0 .- draw)
        end
        return array_of_dicts
    else
        array_of_dicts = Array{Dict{Symbol,Real}}(undef, num)
        number_of_itos = length(covar.covariance_labels_)
        for i in 1:num
            draw =  next!(number_generator)
            array_of_dicts[i] = get_draws(covar; uniform_draw = draw)
        end
        return array_of_dicts
    end
end

# This is most likely useful for bug hunting.
"""
    get_zero_draws(covar::ForwardCovariance)
get a draw of zero for all `ItoIntegral`s. This may be handy for bug hunting.
### Inputs
* covar - the `ForwardCovariance` struct that you want zero draws from
### Outputs
* A `Dict` of zero draws
"""
function get_zero_draws(covar::ForwardCovariance)
    return Dict{Symbol,Real}(covar.covariance_labels_ .=> 0.0)
end
# This is most likely useful for bug hunting.
"""
    get_zero_draws(covar::ForwardCovariance, num::Integer)
get an array of zero draws for all `ItoIntegral`s. May be handy for bug hunting.
### Inputs
* `covar` - the `ForwardCovariance` struct that you want zero draws from
* `num` - The number of zero draws you want.
### Outputs
* A `Vector` of `Dict`s of zero draws
"""
function get_zero_draws(covar::ForwardCovariance, num::Integer)
    array_of_dicts = Array{Dict{Symbol,Real}}(undef, num)
    for i in 1:num
        array_of_dicts[i] = get_zero_draws(covar)
    end
    return array_of_dicts
end

"""
    pdf(covar::ForwardCovariance, coordinates::Dict{Symbol,Real})
get the value of the pdf at some coordinates. Note that it is assumed that
the mean of the multivariate gaussian is the zero vector.
### Inputs
* `covar` - the `ForwardCovariance` struct that you want to evaluate the pdf
* `coordinates` - The coordinates you want to examine.
### Outputs
* A scalar
"""
function pdf(covar::ForwardCovariance, coordinates::Dict{Symbol,Real})
    # The pdf is det(2\pi\Sigma)^{-0.5}\exp(-0.5(x - \mu)^\prime \Sigma^{-1} (x - \mu))
    # Where Sigma is covariance matrix, \mu is means (0 in this case) and x is the coordinates.
    rank_of_matrix = length(covar.covariance_labels_)
    x = get.(Ref(coordinates), covar.covariance_labels_, 0)
    one_on_sqrt_of_det_two_pi_covar     = 1/(sqrt(covar.determinant_) * (2*pi)^(rank_of_matrix/2))
    return one_on_sqrt_of_det_two_pi_covar * exp(-0.5 * x' * covar.inverse_ * x)
end

"""
    log_likelihood(covar::ForwardCovariance, coordinates::Dict{Symbol,Real})
get the log likelihood at some coordinates. Note that it is assumed that
the mean of the multivariate gaussian is the zero vector.
### Inputs
* `covar` - the `ForwardCovariance` struct that you want to evaluate the log likelihood.
* `coordinates` - The coordinates you want to examine.
### Outputs
* A scalar
"""
function log_likelihood(covar::ForwardCovariance, coordinates::Dict{Symbol,R}) where R<:Real
    # The pdf is det(2\pi\Sigma)^{-0.5}\exp(-0.5(x - \mu)^\prime \Sigma(x - \mu))
    # Where Sigma is covariance matrix, \mu is means (0 in this case) and x is the coordinates.
    rank_of_matrix = length(covar.covariance_labels_)
    x = get.(Ref(coordinates), covar.covariance_labels_, 0)
    one_on_sqrt_of_det_two_pi_covar     = -0.5*log(covar.determinant_)  +  (rank_of_matrix/2)*log(2*pi)
    return one_on_sqrt_of_det_two_pi_covar + (-0.5 * x' * covar.inverse_ * x)
end


text:
### Testing quadrature related function
using Jacobi
using Polynomials
@static if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

Nmax = 200
a = rand(1:6, Nmax)  # Coefficient for a ninth degree polynomial

N = 10
Q = 5
p = Poly(a[1:N])
p_int = polyint(p)
X = polyval(p_int, 1.0) - polyval(p_int, -1.0)
q = Quadrature(GJ, Q, 0.0, 0.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# I'm not getting the precision I need for BigFloats. Should investigate further
#X = polyval(p1_int, BigFloat(1.0)) - polyval(p1_int, BigInt(-1.0))
#q2 = Quadrature(GJ, 5, BigFloat(0), BigFloat(0), BigFloat)
#y = polyval(p1, qzeros(q2))
#Xg= sum(qweights(q2) .* y)
#@test abs(Xg-X) 0.0 200*eps(Xg)

# Weights a=1, b=2
p2 = p * Poly([1, -1]) * Poly([1, 1])^2
p2_int = polyint(p2)
X = polyval(p2_int, 1.0) - polyval(p2_int, -1.0)
q = Quadrature(GJ, Q, 1.0, 2.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# Test Gauss-Radau:
# GRJM
N = 9
Q = 5
p = Poly(a[1:N])
p_int = polyint(p)
X = polyval(p_int, 1.0) - polyval(p_int, -1.0)
q = Quadrature(GRJM, Q, 0.0, 0.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# GRJP
q = Quadrature(GRJP, Q, 0.0, 0.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# Weights a=1, b=2 GRJM
p2 = p * Poly([1, -1]) * Poly([1, 1])^2
p2_int = polyint(p2)
X = polyval(p2_int, 1.0) - polyval(p2_int, -1.0)
q = Quadrature(GRJM, Q, 1.0, 2.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# Weights a=1, b=2 GRJP
p2 = p * Poly([1, -1]) * Poly([1, 1])^2
p2_int = polyint(p2)
X = polyval(p2_int, 1.0) - polyval(p2_int, -1.0)
q = Quadrature(GRJP, Q, 1.0, 2.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# Test Gauss-Lobatto
N = 8
Q = 5
p = Poly(a[1:N])
p_int = polyint(p)
X = polyval(p_int, 1.0) - polyval(p_int, -1.0)
q = Quadrature(GLJ, Q, 0.0, 0.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# Weights a=1, b=2 GLJ
p2 = p * Poly([1, -1]) * Poly([1, 1])^2
p2_int = polyint(p2)
X = polyval(p2_int, 1.0) - polyval(p2_int, -1.0)
q = Quadrature(GRJM, Q, 1.0, 2.0, Float64)
y = polyval(p, qzeros(q))
Xg = sum(qweights(q) .* y)
@test Xg ≈ X

# Testing derivatives:
# GJ
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GJ, Q, 0, 0)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# GRJM
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GRJM, Q, 0, 0)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# GRJP
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GRJP, Q, 0, 0)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# GLJ
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GLJ, Q, 0, 0)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# Testing with weights a and b
a1 = 0.5
b1 = 0.9

# GJ
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GJ, Q, a1, b1)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# GRJM
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GRJM, Q, a1, b1)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# GRJP
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GRJP, Q, a1, b1)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# GLJ
Q = 10
N = 10
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GLJ, Q, a1, b1)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

Q = 200
N = 200
p = Poly(a[1:N])
dp = polyder(p)
q = Quadrature(GLJ, Q, a1, b1)
z = qzeros(q)
y = polyval(p, z)
dy = polyval(dp, z)
dyg = qdiff(q) * y
@test dy ≈ dyg

# Test interpolation:
# GJ

N = 10
Q = 10
p = Poly(a[1:N])
x = -1.0:0.02:1.0
ue = polyval(p, x)
q = Quadrature(GLJ, Q, 1.0, 0.5)
z = qzeros(q)
u = polyval(p, z)
I = interp_mat(x, z)
ui = I * u
@test ui ≈ ue


text:
## This file contains functions and types used for in slico digestioon
## define restriction enzyme [name, cut-sequence, cut-site, overhang]
mutable struct restrictionEnzyme
    name::String # name of the restriction enzyme
    cutSeq::Array{String} # recognition sequence
    cutPlace::Array{Int64} # where it gets cut
    overhang::Array{String} # overhang of the restriction enzyme if exists
    restrictionEnzyme(name, cutSeq, cutPlace, overhang) =
        new(name, cutSeq, cutPlace, overhang)
end

## list of restraiction enzymes
ApeKI = restrictionEnzyme("ApeKI", ["GCAGC", "GCTGC"], [2, 2], ["CTG", "CAG"]); # Elshire et al. 2011
PstI = restrictionEnzyme("PstI", ["CTGCAG"], [2], ["CTGCA"]); # Poland et al. 2012
MspI = restrictionEnzyme("MspI", ["CCGG"], [2], ["CGG"]); # Poland et al. 2012

## nucleotides and its reverse complements
nucl = Dict(
    'R' => rand(['A', 'G']),
    'Y' => rand(['C', 'T']),
    'S' => rand(['G', 'C']),
    'W' => rand(['A', 'T']),
    'K' => rand(['G', 'T']),
    'M' => rand(['A', 'C']),
    'B' => rand(['C', 'G', 'T']),
    'D' => rand(['A', 'G', 'T']),
    'H' => rand(['A', 'C', 'T']),
    'V' => rand(['A', 'C', 'G']),
); # randomly assign IUPAC codes for uncertain ones in the assembly
comp = Dict('A' => 'T', 'C' => 'G', 'G' => 'C', 'T' => 'A'); # complements of DNA nucleotides

## function to split sequence
function splitSequence(seq, cutSeq, startSeq, endSeq)
    cuts = split(uppercase(seq), cutSeq) # split DNA sequence
    if length(cuts) > 1
        cuts[2:length(cuts)] = startSeq .* cuts[2:length(cuts)] # re-attach the strat seq to the fragment after cut-place
        cuts[1:(length(cuts)-1)] = cuts[1:(length(cuts)-1)] .* endSeq # re-attach the end seq to the fragment before cut-place
    else
        cuts
    end
    cuts
end

## function to digest sequence using specified restriction enzyme
function digest(sequence, re::restrictionEnzyme)
    cutSeq = re.cutSeq
    cutPlace = re.cutPlace
    overhang = re.overhang
    startSeq = [SubString(cutSeq[i], cutPlace[i]) for i = 1:length(cutSeq)]
    endSeq =
        [SubString(cutSeq[i], 1, cutPlace[i] - 1) * overhang[i] for i = 1:length(cutSeq)]
    for k = 1:length(cutSeq)
        cuts = sequence
        sequence = []
        for s in cuts
            sequence = [sequence; splitSequence(s, cutSeq[k], startSeq[k], endSeq[k])]
        end
        sequence
    end
    if length(sequence) > 2
        sequence[2:length(sequence)]
    else
        sequence
    end
end

## function to implement double-digestion (or enzyme has more than one recongnition site)
function digestSecond(sequence, re::restrictionEnzyme)
    cutSeq = re.cutSeq
    cutPlace = re.cutPlace
    overhang = re.overhang
    startSeq = [SubString(cutSeq[i], cutPlace[i]) for i = 1:length(cutSeq)]
    endSeq =
        [SubString(cutSeq[i], 1, cutPlace[i] - 1) .* overhang[i] for i = 1:length(cutSeq)]
    goodSeq = []
    for k = 1:length(sequence)
        s = sequence[k]
        sp = splitSequence(uppercase(s), cutSeq[1], startSeq[1], endSeq[1])
        if length(sp) > 1
            if k == 1
                s2 = startSeq .* sp[length(sp)]
                goodSeq = [goodSeq; s2]
            elseif k == length(sequence)
                s1 = sp[1] .* endSeq
                goodSeq = [goodSeq; s1]
            else
                s1 = sp[1] .* endSeq
                s2 = startSeq .* sp[length(sp)]
                goodSeq = [goodSeq; s1; s2]
            end
        else
            goodSeq = [goodSeq; sp]
        end
    end
    goodSeq
end

## function to implement virtual digestion
function digestGenome(
    genofile::String,
    re::Array{restrictionEnzyme,1},
    useChr::Array{Int64},
    useChrLen::Array{Float64},
    lowerThresh::Int64,
    upperThresh::Int64,
    plotOutput::Bool,
    writeOutput::Bool,
)
    ## module 1: read genome and cut it into fragments
    genome = readdlm(GZip.open(genofile), '\t') # load genome
    index = LinearIndices(genome) # index genome
    header = index[findall(x -> x[1] == '>', genome)] # identify headers within each FASTA file
    global chrLen = Array{Float64}(undef, 0) # empty array to store estimated chromosome length
    seqFrag = Array{String}(undef, 0) # empty array to store the DNA sequence
    chrFrag = Array{Int64}(undef, 0) # empty array to store chromosome
    idFrag = Array{Int64}(undef, 0) # empty array to store the ID
    posFrag = Array{Int64}(undef, 0) # empty array to store the position
    lenFrag = Array{Int64}(undef, 0) # empty array to store the position
    ## (optional) subsetting genome
    if length(useChr) == 1
        global numChr = useChr[1] # useChr is considered to be # of chromosomes if an integer is supplied
        useChr = collect(1:useChr[1]) # select first n chromosomes
        starts = (header.+1)[useChr] # subset genome
        ends = [header[2:end] .- 1; size(genome)[1]][useChr]
        println("INFO: Using chromsome $(useChr)!")
    elseif length(useChr) == length(header)
        global numChr = length(useChr) # useChr is considered to be a set of chromosomes if a vector is supplied
        starts = (header.+1)[useChr] # select chromosomes (specified by [useChr])
        ends = [header[2:end] .- 1; size(genome)[1]][useChr] # subset genome
        println("INFO: Using chromsome $(useChr) in the simulation!")
    else
        error("Wrong chromosome number/set specified, please check your parameters!")
    end
    for c = 1:numChr
        seq = [join(genome[starts[c]:ends[c]])] # read genomic sequence
        ## virtual digestion
        if size(re, 1) == 1
            frag = digest(seq, re[1]) # single-digestion
            reName = re[1].name
        else
            raw = digest(seq, re[1])
            frag = digestSecond(raw, re[2]) # double-digestion if more than one enzymes is supplied
            reName = [re[1].name re[2].name]
        end
        numFrag = length(frag) # number of fragments
        println(
            "CHROMOSOME $c: $numFrag GBS fragments are generated through virtual digestion using $(reName).",
        )
        len = length.(frag) # length of fragments
        pos = [0; cumsum(len)[1:(end-1)]] .+ 1 # starting position of each fragment
        chrLenEst = pos[end] / 1e6 # estimate chromosome length
        println("CHROMOSOME $c: Estimated chromosome length equals to $chrLenEst Mb.")
        ## (optional) subsetting chromosome
        if useChrLen == [] || size(useChrLen, 1) < numChr
            keep = [1:length(frag)...] # keep all fragments if subsetting is not required
            chrLen = [chrLen; chrLenEst] # store chromosome lengths
            println(
                "CHROMOSOME $c: Subsetting chromosome is not required. Keeping all $numFrag available GBS fragments.",
            )
        elseif size(useChrLen, 1) == numChr || useChrLen[c] <= estLen
            keep = findall(x -> x < useChrLen[c] * 1e6, pos) # keep only fragments found within specified length (varied at each chromosome)
            chrLen = useChrLen # store user-defined chromosome lengths
            println(
                "CHROMOSOME $c (up to $(useChrLen[c]) cM): $(length(keep)) out of $numFrag available GBS fragments are kept.",
            )
        elseif size(useChrLen, 1) == 1 # keep only fragments found within specified length (uniform across chromosomes)
            keep = findall(x -> x < useChrLen[1] * 1e6, pos)
            chrLen = [chrLen; useChrLen[1]]
            println(
                "CHROMOSOME $c (up to $(useChrLen[c]) cM): $(length(keep)) out of $numFrag available GBS fragments are kept.",
            )
        else
            error(
                "CHROMOSOME $c: Wrong chromosome length(s) specified, please check your parameters!",
            )
        end
        idFrag = [idFrag; [1:length(keep)...]] # identifier of each fragment
        chrFrag = [chrFrag; fill(c, length(keep))] # origin of chromosome per fragment
        posFrag = [posFrag; pos[keep]] # starting position per fragment
        lenFrag = [lenFrag; len[keep]] # length of each fragment
        seqFrag = [seqFrag; frag[keep]] # genomic sequence of each fragment
    end
    ## module 2: fragment size-selection
    fragRaw = [idFrag chrFrag posFrag lenFrag seqFrag] # dataset stores all fragments
    numFragTotal = length(seqFrag) # total number of GBS fragments
    selected = findall(x -> lowerThresh <= x <= upperThresh, lenFrag) # fragment size-selection
    numFragSelected = length(selected) # number of fragments selected in size-selection
    println(
        "INFO: $numFragSelected out of $numFragTotal GBS fragments are selected after size-selection with lower and upper thresholds equals to [$lowerThresh, $upperThresh].",
    )
    fragGBS = fragRaw[selected, :] # dataset conatins GBS fragments (after size-selection)
    global GBSFrag = (
        frag = fragGBS[:, 5],
        len = fragGBS[:, 4],
        pos = fragGBS[:, 3],
        chr = fragGBS[:, 2],
        id = fragGBS[:, 1],
    )
    ## module 3: estimate coverage of selected GBS fragments (per Mb)
    GBSCoverChr = Array{Int64}(undef, 0)
    GBSFragWinStatrs = Array{Int64}(undef, 0)
    GBSFragWinEnds = Array{Int64}(undef, 0)
    GBSFragCover = Array{Float64}(undef, 0)
    for c = 1:numChr
        numWin = Int(round(chrLen[c]))
        winStarts = Int.([0:(numWin-1)...] .* 1e6)
        winEnds = Int.([1:numWin...] .* 1e6)
        GBSFragChr = findall(x -> x == c, GBSFrag.chr)
        GBSFragPos = GBSFrag.pos[GBSFragChr]
        GBSFragLen = GBSFrag.len[GBSFragChr]
        GBSFragWin = [
            findall(x -> x <= winEnds[w] && x >= winStarts[w], GBSFragPos) for w = 1:numWin
        ]
        GBSFragCoverWin = [sum(GBSFragLen[GBSFragWin[w]]) / 1e6 for w = 1:numWin]
        GBSFragCover = [GBSFragCover; GBSFragCoverWin]
        GBSFragWinEnds = [GBSFragWinEnds; winEnds]
        GBSFragWinStatrs = [GBSFragWinStatrs; winStarts]
        GBSCoverChr = [GBSCoverChr; fill(c, numWin)]
    end
    GBSCover = [GBSCoverChr GBSFragWinStatrs GBSFragWinEnds GBSFragCover]
    println(
        "INFO: Expected sequencing coverage based on $numFragSelected selected GBS fragments is approximately $(round(mean(GBSFragCover)*100;digits=2))%.",
    )
    ## module 4: writing out and plotting fragment-szie distribution and GBS coverage if required
    if plotOutput == true
        histogram(
            fragRaw[:, 4],
            normalize = :probability,
            title = "Size Distribution of Raw GBS Fragements",
            label = "",
            fmt = :png,
        )
        savefig("RawFragSize")
        histogram(
            fragGBS[:, 4],
            normalize = :probability,
            title = "Size Distribution of GBS Fragements (after Size-Selection)",
            label = "",
            fmt = :png,
        )
        savefig("GBSFragSize")
        histogram(
            GBSFragCover,
            normalize = :probability,
            title = "Coverage of GBS Fragments per Mb",
            label = "",
            fmt = :png,
        )
        savefig("GBSCoverage")
    end
    if writeOutput == true
        writedlm("RawFrag.txt", fragRaw)
        writedlm("GBSFrag.txt", fragGBS)
        writedlm("GBSCoverage.txt", GBSCover)
    end
end


text:
# ------------------------------------------------------------------------------------------
# Surrogate functions
# ------------------------------------------------------------------------------------------
Base.broadcastable(c::Surrogate) = Ref(c)

struct Hinge{T<:Real} <: Surrogate
    ϑ::T
end

value(l::Hinge, s) = max(0, 1 + l.ϑ*s)

struct Quadratic{T<:Real} <: Surrogate
    ϑ::T
end

value(l::Quadratic, s) = max(0, 1 + l.ϑ*s)^2

# ------------------------------------------------------------------------------------------
# Model state
# ------------------------------------------------------------------------------------------
mutable struct State{T<:Real}
    s::Vector{T}
    αβ::Vector{T}
    δ::T
    αsum::T
    βsort::Vector{T}

    State(T) = new{T}()
end

function find_βmax(model::Model, βk)
    βsort = model.state.βsort
    return βsort[1] != βk ? βsort[1] : βsort[2]
end

function update_βmax!(model::Model, K::KernelMatrix, k, β, βnew)
    k > K.nα || return

    βsort = model.state.βsort
    deleteat!(βsort, searchsortedfirst(βsort, β; rev = true))
    insert!(βsort, searchsortedfirst(βsort, βnew; rev = true), βnew)
    return
end

function update_αsum!(model::Model, K::KernelMatrix, k, l, Δ)
    if k <= K.nα && l > K.nα
        model.state.αsum += Δ
    end
    return
end

# ------------------------------------------------------------------------------------------
# Update rule
# ------------------------------------------------------------------------------------------
struct UpdateRule{T<:Real}
    L::T
    Δ::T
    num::T
    den::T
    Δlb::T
    Δub::T
    k::Int
    l::Int
    δ::T
end

UpdateRule(k, l; T = Float32) = UpdateRule{T}(-Inf, 0, 0, 0, 0, 0, k, l, 0)

function UpdateRule(num::T, den::Real, lb::Real, ub::Real, k::Int, l::Int) where T
    Δ = compute_Δ(num, den, lb, ub)
    L = - den*Δ^2/2 - num*Δ

    return UpdateRule{T}(L, Δ, num, den, lb, ub, k, l, 0)
end

function UpdateRule(
    model::AbstractPatMat,
    K::KernelMatrix,
    num::T,
    den::Real,
    lb::Real,
    ub::Real,
    k::Int,
    l::Int,
    δ::Real,
) where T

    Δ = compute_Δ(num, den, lb, ub)
    L = - den*Δ^2/2 - num*Δ - (δ - model.state.δ)*K.nβ*model.τ

    return UpdateRule{T}(L, Δ, num, den, lb, ub, k, l, δ)
end

compute_Δ(num, den, lb, ub) = min(max(lb, - num/den), ub)

function select(r1::UpdateRule{T}, r2::UpdateRule{T}) where T
    r1.Δ == 0 || isinf(r1.L) || isnan(r1.L) && return r2
    r2.Δ == 0 || isinf(r2.L) || isnan(r2.L) && return r1

    return if r1.L == r2.L
        rand(Bool) ? r1 : r2
    elseif r1.L > r2.L
        r1
    else
        r2
    end
end

function update!(model::Model, K::KernelMatrix; k  = rand(1:K.n))
    best = UpdateRule(k, k; T = eltype(K))
    K[k,1]
    for l in 1:K.n
        l == k && continue
        update = if k <= K.nα && l <= K.nα
            rule_αα(model, K, k, l)
        elseif k <= K.nα && l > K.nα
            rule_αβ(model, K, k, l)
        elseif k > K.nα && l <= K.nα
            rule_αβ(model, K, l, k)
        else
            rule_ββ(model, K, k, l)
        end
        best = select(best, update)
    end
    update!(model, K, best)
    return k == best.l ? (best.Δ, best.l, best.k) : (best.Δ, best.k, best.l)
end

function update!(model::Model, K::KernelMatrix, update::UpdateRule)
    iszero(update.Δ) && return
    @unpack k, l, Δ = update

    y = (k <= K.nα && l > K.nα) ? -1 : 1
    if !isa(model, TopPush)
        αβ = model.state.αβ
        update_βmax!(model, K, k, αβ[k], αβ[k] + Δ)
        update_βmax!(model, K, l, αβ[l], αβ[l] - y*Δ)
        update_αsum!(model, K, k, l, Δ)
    end

    model.state.s .+= Δ .* (K[k, :] - y*K[l, :])
    model.state.αβ[k] += Δ
    model.state.αβ[l] -= y*Δ
    if isa(model, AbstractPatMat)
        model.state.δ = update.δ
    end
    return
end


text:
function blink(t::D3Tree)
    w = Window()
    str = stringmime(MIME("text/html"), t)
    # println(str)
    body!(w, str)
    return w
end

function inchrome(t::D3Tree)
    fname = joinpath(mktempdir(), "tree.html")
    open(fname, "w") do f
        show(f, MIME("text/html"), t)
    end
    if Sys.iswindows()
        run(`cmd /C start chrome "$fname"`)
    elseif Sys.isapple()
        run(`open -a "Google Chrome" $fname`)
    else
        run(`google-chrome $fname`)
    end
end

"""
    inbrowser(t::D3Tree, browsername::String)

Open an html document with the D3Tree in a browser with a platform-specific command.
"""
function inbrowser(t::D3Tree, browsername::String)
    if Sys.iswindows()
        inbrowser(t, `cmd /C start $browsername`)
    elseif Sys.isapple()
        inbrowser(t, `open -a $browsername`)
    else
        inbrowser(t, `$browsername`)
    end
end

"""
    inbrowser(t::D3Tree, command::Cmd)

Open an html document with the D3Tree in a program launched with the specified command.
"""
function inbrowser(t::D3Tree, command::Cmd)
    fname = joinpath(mktempdir(), "tree.html")
    open(fname, "w") do f
        show(f, MIME("text/html"), t)
    end
    run(`$command $fname`)
end


text:
"""
$(SIGNATURES)

Compute `log(sum(exp, X))` in a numerically stable way that avoids intermediate over- and
underflow.

`X` should be an iterator of real numbers. The result is computed using a single pass over
the data.

# References

[Sebastian Nowozin: Streaming Log-sum-exp Computation](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)
"""
logsumexp(X) = _logsumexp_onepass(X)

"""
$(SIGNATURES)

Compute `log.(sum(exp.(X); dims=dims))` in a numerically stable way that avoids
intermediate over- and underflow.

The result is computed using a single pass over the data.

# References

[Sebastian Nowozin: Streaming Log-sum-exp Computation](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)
"""
logsumexp(X::AbstractArray{<:Real}; dims=:) = _logsumexp(X, dims)

_logsumexp(X::AbstractArray{<:Real}, ::Colon) = _logsumexp_onepass(X)
function _logsumexp(X::AbstractArray{<:Real}, dims)
    # Do not use log(zero(eltype(X))) directly to avoid issues with ForwardDiff (#82)
    FT = float(eltype(X))
    xmax_r = reduce(_logsumexp_onepass_op, X; dims=dims, init=(FT(-Inf), zero(FT)))
    return @. first(xmax_r) + log1p(last(xmax_r))
end

function _logsumexp_onepass(X)
    # fallback for empty collections
    isempty(X) && return log(sum(X))
    return _logsumexp_onepass_result(_logsumexp_onepass_reduce(X, Base.IteratorEltype(X)))
end

# function barrier for reductions with single element and without initial element
_logsumexp_onepass_result(x) = float(x)
_logsumexp_onepass_result((xmax, r)::Tuple) = xmax + log1p(r)

# iterables with known element type
function _logsumexp_onepass_reduce(X, ::Base.HasEltype)
    # do not perform type computations if element type is abstract
    T = eltype(X)
    isconcretetype(T) || return _logsumexp_onepass_reduce(X, Base.EltypeUnknown())

    FT = float(T)
    return reduce(_logsumexp_onepass_op, X; init=(FT(-Inf), zero(FT)))
end

# iterables without known element type
_logsumexp_onepass_reduce(X, ::Base.EltypeUnknown) = reduce(_logsumexp_onepass_op, X)

## Reductions for one-pass algorithm: avoid expensive multiplications if numbers are reduced

# reduce two numbers
function _logsumexp_onepass_op(x1, x2)
    a = x1 == x2 ? zero(x1 - x2) : -abs(x1 - x2)
    xmax = x1 > x2 ? oftype(a, x1) : oftype(a, x2)
    r = exp(a)
    return xmax, r
end

# reduce a number and a partial sum
function _logsumexp_onepass_op(x, (xmax, r)::Tuple)
    a = x == xmax ? zero(x - xmax) : -abs(x - xmax)
    if x > xmax
        _xmax = oftype(a, x)
        _r = (r + one(r)) * exp(a)
    else
        _xmax = oftype(a, xmax)
        _r = r + exp(a)
    end
    return _xmax, _r
end
_logsumexp_onepass_op(xmax_r::Tuple, x) = _logsumexp_onepass_op(x, xmax_r)

# reduce two partial sums
function _logsumexp_onepass_op((xmax1, r1)::Tuple, (xmax2, r2)::Tuple)
    a = xmax1 == xmax2 ? zero(xmax1 - xmax2) : -abs(xmax1 - xmax2)
    if xmax1 > xmax2
        xmax = oftype(a, xmax1)
        r = r1 + (r2 + one(r2)) * exp(a)
    else
        xmax = oftype(a, xmax2)
        r = r2 + (r1 + one(r1)) * exp(a)
    end
    return xmax, r
end


text:
using HDF5
using JSON
using Flux

"""
Loads the weights from the weights.h5 file.
"""
function weights(file="weights.h5")
    f = h5open(file, "r")
#    for ele in keys(read(f))
#        if contains(ele ,"/")
#            temp = split(ele, "/")
#            read(f)[temp[1]][temp[1]] = read(f)[temp[1]][temp[2]][temp[1]][temp[2]]
#        end
#    end
    weight = read(f)
    for ele in keys(weight)
        if !isempty(weight[ele])
            if haskey(weight[ele], ele*"_1")
                for ele2 in keys(weight[ele][ele*"_1"])
                    weight[ele][ele*"_1"][ele2] = convert(Array{Float64, N} where N, weight[ele][ele*"_1"][ele2])
                end
            elseif haskey(weight[ele], ele)
                for ele2 in keys(weight[ele][ele])
                    weight[ele][ele][ele2] = convert(Array{Float64, N} where N, weight[ele][ele][ele2])
                end
            elseif ele == "conv1"
                #weight[ele]["conv"][ele]["conv"] = convert(Array{Float64, N} where N, weight[ele]["conv"][ele]["conv"]["kernel:0"])
                weight[ele*"/conv"] = Dict{Any, Any}()
                weight[ele*"/conv"][ele*"/conv"] = Dict{Any, Any}()
                weight[ele*"/conv"][ele*"/conv"]["kernel:0"] = convert(Array{Float64, N} where N, weight[ele]["conv"][ele]["conv"]["kernel:0"])
            end
        end
    end
    return weight
end

"""
Load the structure of the model from the .JSON
file.
"""
function load_structure(file="structure.json")
    res = JSON.parse(String(read(open(file, "r"))))
    return res["config"]
end

"""
Check if model uses sequential/functional API.
"""
function check_modeltype(file)
    res = JSON.parse(String(read(open(file, "r"))))
    return res["class_name"]
end

struct new_type
    layer_type::Symbol
    input_nodes::Any
    fields::Any
end

"""
Get structure information from class_name
"""
function layer_type(a)
    if (a["class_name"] == "Conv2D" || a["class_name"] == "Conv3D")
        return :Conv
    elseif (a["class_name"] == "MaxPooling2D" || a["class_name"] == "MaxPooling3D")
        return :MaxPool
    else
        return Symbol(a["class_name"])
    end
    #    if (a["class_name"] == "Dropout")
    #    return :Dropout
    #elseif (a["class_name"] == "Flatten")
    #    return :Flatten
    #elseif (a["class_name"] == "Dense")
    #    return :Dense
    #elseif (a["class_name"] == "Activation")
    #    return :Activation
    #elseif (a["class_name"] == "Reshape")
    #    return :Reshape
    #elseif (a["class_name"] == "BatchNormalization")
    #    return :BatchNormalization
    #elseif (a["class_name"] == "InputLayer")
    #    return :InputLayer
    #elseif (a["class_name"] == "ZeroPadding2D")
    #    return :ZeroPadding2D
    #elseif (a["class_name"] == "Add")
    #    return :Add
    #elseif (a["class_name"] == "AveragePooling2D")
    #    return :AveragePooling2D
    #elseif (a["class_name"] == "Concatenate")
    #    return :Concatenate
    #elseif (a["class_name"] == "GlobalAveragePooling2D")
    #    return :GlobalAveragePooling2D
    #elseif (a["class_name"] == "Embedding")
    #    return :Embedding
    #elseif (a["class_name"] == "LSTM")
    #    return :LSTM
    #elseif (a["class_name"] == "MaxPooling1D")
    #    return :MaxPooling1D
    #elseif (a["class_name"] == "AveragePooling1D")
    #    return :AveragePooling1D
    #elseif (a["class_name"] == "Conv1D")
    #    return :Conv1D
    #end
end

"""
Extract necessary fields only from the type of layer.
"""
function fields(a)
    if layer_type(a) == :Conv
        return ["name", "strides", "activation", "kernel_size", "padding"]
    elseif layer_type(a) == :MaxPool
        return ["name", "strides", "padding", "pool_size"]
    elseif layer_type(a) == :Dropout
        return ["name", "rate"]
    elseif layer_type(a) == :Flatten
        return ["name"]
    elseif layer_type(a) == :Dense
        return ["name", "activation"]
    elseif layer_type(a) == :relu
        return ["name", "activation"]
    elseif layer_type(a) == :softmax
        return ["name", "activation"]
    elseif layer_type(a) == :Reshape
        return ["name", "target_shape"]
    elseif layer_type(a) == :BatchNormalization
        return ["name", "momentum", "epsilon"]
    elseif layer_type(a) == :InputLayer
        return ["name"]
    elseif layer_type(a) == :ZeroPadding2D
        return ["name", "padding"]
    elseif layer_type(a) == :Add
        return ["name"]
    elseif layer_type(a) == :AveragePooling2D
        return ["name","strides","padding","pool_size"]
    elseif layer_type(a) == :Concatenate
        return ["name","axis"]
    elseif layer_type(a) == :GlobalAveragePooling2D
        return ["name"]
    elseif layer_type(a) == :Embedding
        return ["name", "output_dim", "input_dim", "input_length"]
    elseif layer_type(a) == :LSTM
        return ["name", "activation", "units", "return_sequences"]
    elseif layer_type(a) == :Activation
        return ["name", "activation"]
    elseif layer_type(a) == :MaxPooling1D
        return ["name","strides", "pool_size", "padding"]
    elseif layer_type(a) == :AveragePooling1D
        return ["name","strides", "pool_size", "padding"]
    elseif layer_type(a) == :Conv1D
        return ["name", "activation", "strides", "dilation_rate", "kernel_size", "filters", "padding"]
    end
end

"""
Extract data layerwise from the structure["layers"]
array.
"""
function load_layers(a::Array{Any, 1})
    res = Array{Any, 1}()
    for ele in a
        d = Dict{Any, Any}()
        for ele2=1:length(fields(ele))
            d[fields(ele)[ele2]] = ele["config"][fields(ele)[ele2]]
        end
        if !haskey(ele, "inbound_nodes")
            ele["inbound_nodes"] = nothing
        end
        inputs = ele["inbound_nodes"]
    push!(res, new_type(layer_type(ele), inputs, d))
    end
    return res
end


text:
module SemiLagMPI

import Base.Threads: @sync, @spawn, nthreads, threadid

include("mpiinterface.jl")

include("util.jl")
include("fftbig.jl")
include("mesh.jl")
include("interpolation.jl")
include("lagrange.jl")
include("spline.jl")
include("bspline.jl")
include("bsplinelu.jl")
include("bsplinefft.jl")
include("advection.jl")
include("util_poisson.jl")
include("poisson.jl")
include("rotation.jl")
include("translation.jl")

export UniformMesh, Advection1d, Advection1dData, AbstractInterpolation, AbstractInterpolation2d
export Lagrange, B_SplineLU, B_SplineFFT, interpolate!
export compute_charge!,
    compute_elfield!, compute_elfield, compute_ee, compute_ke, advection!
export dotprod, getpoissonvar, getrotationvar, gettranslationvar
export TimeOptimization, NoTimeOpt, SimpleThreadsOpt, SplitThreadsOpt, MPIOpt
export get_type, sizeall, getdata

end


text:
function check_vote(name)
	if voted.get(name)
		println("kick them out!")
	else
		voted[name] = true
		println("let them vote!")
	end
end

voted = Dict()
check_vote("tom")
check_vote("mike")
check_vote("mike")


text:
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/build_literate/Simulation_Pulse.ipynb)

# # RF-Pulse Simulation
# The following code replicates the RF-pulse simulation of Fig. 3 and plots the ``z^s``-magnetization at the end of respective pulse.

# For these simulations we need the following packages:
using MRIgeneralizedBloch
using QuadGK
using DifferentialEquations
using SpecialFunctions
using Plots
plotlyjs(bg = RGBA(31/255,36/255,36/255,1.0), ticks=:native); #hide #!nb

# and we simulate an isolated semi-solid spin pool with the following parameters:
R₁ = 1 # 1/s
T₂ˢ = 10e-6; # s

# Here, we simulate π-pulses with the following parameters:
α = π
Tʳᶠ = exp.(range(log(2e-7), log(1e-1), length=100)) # s
ω₁ = α ./ Tʳᶠ # rad/s
ω₀ = 0; # rad/s
# Replace first line with `α = π/4` or `α = π/2` to simulate the other two rows of Fig. 3.

# ## Lorentzian Lineshape
# In this script, we simulate the three lineshapes separately, starting with the Lorentzian lineshape for which the Bloch model provides a ground truth.

# ### Bloch Model
# We can formulate the [Bloch model](http://dx.doi.org/10.1103/PhysRev.70.460) as
# ```math
# \partial_t \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} = \begin{pmatrix}
# -R_2 & -ω_0 & ω_1 & 0 \\
# ω_0 & -R_2 & 0 & 0 \\
# -ω_1 & 0 & -R_1 & R_1 \\
# 0 & 0 & 0 & 0
# \end{pmatrix} \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} ,
# ```
# where the matrix is the Hamiltonian of the Bloch model. For a constant ``ω_0`` and ``ω_1``, we can evaluate the Bloch model by taking the  matrix exponential of its Hamiltonian:
H(ω₁, ω₀, R₂, R₁) = [-R₂ -ω₀  ω₁  0;
                      ω₀ -R₂   0  0;
                     -ω₁   0 -R₁ R₁;
                       0   0   0  0]

z_Bloch = similar(Tʳᶠ)
for i = 1:length(Tʳᶠ)
    (_, _, z_Bloch[i], _)  = exp(H(ω₁[i], ω₀, 1 / T₂ˢ, R₁) * Tʳᶠ[i]) * [0; 0; 1; 1]
end

# ### Graham's Spectral Model
# [Graham's spectral model](http://doi.org/10.1002/jmri.1880070520) is derived by integrating over the lineshape multiplied by the spectral response function of the RF-pulse. This results in the RF-induced saturation rate `Rʳᶠ` that is used in an exponential model:

Rʳᶠ = @. ω₁^2 * T₂ˢ * ((exp(-Tʳᶠ / T₂ˢ) -1) * T₂ˢ + Tʳᶠ) / Tʳᶠ
z_Graham_spec_Lorentzian = @. (Rʳᶠ * exp(-Tʳᶠ * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ);

# ### Graham's Single Frequency Approximation
# In the [single frequency approximation](http://doi.org/10.1002/jmri.1880070520), Graham assumes that the RF-pulse has only a single frequency, which reduces `Rʳᶠ` to

g_Lorentzian(ω₀) = T₂ˢ / π ./ (1 .+ (T₂ˢ .* ω₀).^2)
Rʳᶠ = @. π * ω₁^2 * g_Lorentzian(ω₀)
z_Graham_SF_approx_Lorentzian = @. (Rʳᶠ * exp(-Tʳᶠ * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ);

# where `g_Lorentzian(ω₀)` denotes the Lorentzian lineshape.

# ### Sled's Model
# [Sled's model](http://dx.doi.org/10.1006/jmre.2000.2059) is given by the ordinary differential equation (ODE)
# ```math
# \partial_t z(t) = \left(-\pi \int_0^t G(t-τ) \omega_1(τ)^2 dτ \right) z(t)  + R_1 (1-z),
# ```
# where ``G(t-τ)`` is the Green's function. The Hamiltonian of this ODE is implemented in [`apply_hamiltonian_sled!`](@ref) and the ODE can be solved with the [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/) package:

z₀ = [1] # initial z-magnetization
z_Sled_Lorentzian = similar(Tʳᶠ)
for i = 1:length(Tʳᶠ)
    param = (ω₁[i], 1, ω₀, R₁, T₂ˢ, greens_lorentzian)
    prob = ODEProblem(apply_hamiltonian_sled!, z₀, (0, Tʳᶠ[i]), param)
    z_Sled_Lorentzian[i] = solve(prob)[end][1]
end

# ### Generalized Bloch Model
# The generalized Bloch model is given by the integro-differential equation (IDE)
# ```math
# \partial_t z(t) = - ω_1(t) \int_0^t G(t,τ) ω_1(τ) z(τ) dτ + R_1 (1 - z(t)) ,
# ```
# or by
# ```math
# \partial_t z(t) = - ω_y(t) \int_0^t G(t,τ) ω_y(τ) z(τ) dτ - ω_x(t) \int_0^t G(t,τ) ω_x(τ) z(τ) dτ + R_1 (1 - z(t)) ,
# ```
# for off-resonant RF-pulses with ``ω_1 = ω_x + i ω_y``. The Hamiltonian of the IDE is implemented in [`apply_hamiltonian_gbloch!`](@ref) and we can solve this IDE with the [delay-differential equation (DDE)](https://diffeq.sciml.ai/stable/tutorials/dde_example/) solver of the [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/) package:

z_fun(p, t) = [1.0]; # initialize history function (will be populated with an interpolation by the DDE solver)

z_gBloch_Lorentzian = similar(Tʳᶠ)
for i = 1:length(Tʳᶠ)
    param = (ω₁[i], 1, ω₀, R₁, T₂ˢ, greens_lorentzian)
    prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, z_fun, (0, Tʳᶠ[i]), param)
    z_gBloch_Lorentzian[i] = solve(prob)[end][1]
end

# Now we have solved all five models and can plot the solutions for comparison:

p = plot(xaxis=:log, legend=:bottomright, xlabel="Tʳᶠ [s]", ylabel="zˢ(Tʳᶠ)")
plot!(p, Tʳᶠ, z_gBloch_Lorentzian, label="generalized Bloch model")
plot!(p, Tʳᶠ, Tʳᶠ .* 0 .+ cos(α), label="cos(α)")
plot!(p, Tʳᶠ, z_Sled_Lorentzian, label="Sled's model")
plot!(p, Tʳᶠ, z_Graham_spec_Lorentzian, label="Graham's spectral model")
plot!(p, Tʳᶠ, z_Graham_SF_approx_Lorentzian, label="Graham's single frequency approximation")
plot!(p, Tʳᶠ, z_Bloch, label="Bloch model")
#md Main.HTMLPlot(p) #hide

# ## Gaussian Lineshape
# We can repeat these simulations (with the exception of the Bloch model) for the Gaussian lineshape:

Rʳᶠ = @. ω₁^2 * T₂ˢ * (2 * T₂ˢ * (exp(-(Tʳᶠ/T₂ˢ)^2/2)-1) + sqrt(2π) * Tʳᶠ * erf(Tʳᶠ/T₂ˢ/sqrt(2))) / (2 * Tʳᶠ)
z_Graham_spec_Gaussian = @. (Rʳᶠ * exp(-Tʳᶠ * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ)

g_Gaussian(ω₀) = T₂ˢ / sqrt(2π) * exp(-(T₂ˢ * ω₀)^2 / 2) # lineshape
Rʳᶠ = @. π * ω₁^2 * g_Gaussian(ω₀)
z_Graham_SF_approx_Gaussian = @. (Rʳᶠ * exp(-Tʳᶠ * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ)

z_Sled_Gaussian = similar(Tʳᶠ)
for i = 1:length(Tʳᶠ)
    param = (ω₁[i], 1, ω₀, R₁, T₂ˢ, greens_gaussian)
    prob = ODEProblem(apply_hamiltonian_sled!, z₀, (0, Tʳᶠ[i]), param)
    z_Sled_Gaussian[i] = solve(prob)[end][1]
end

z_gBloch_Gaussian = similar(Tʳᶠ)
for i = 1:length(Tʳᶠ)
    param = (ω₁[i], 1, ω₀, R₁, T₂ˢ, greens_gaussian)
    prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, z_fun, (0, Tʳᶠ[i]), param)
    z_gBloch_Gaussian[i] = solve(prob)[end][1]
end

p = plot(xaxis=:log, legend=:bottomright, xlabel="Tʳᶠ [s]", ylabel="zˢ(Tʳᶠ)")
plot!(p, Tʳᶠ, z_gBloch_Gaussian, label="generalized Bloch model")
plot!(p, Tʳᶠ, Tʳᶠ .* 0 .+ cos(α), label="cos(α)")
plot!(p, Tʳᶠ, z_Sled_Gaussian, label="Sled's model")
plot!(p, Tʳᶠ, z_Graham_spec_Gaussian, label="Graham's spectral model")
plot!(p, Tʳᶠ, z_Graham_SF_approx_Gaussian, label="Graham's single frequency approximation")
#md Main.HTMLPlot(p) #hide

# ## Super-Lorentzian Lineshape
# Further, we can repeat these simulations for the [super-Lorentzian lineshape](http://dx.doi.org/10.1002/mrm.1910330404) with the exception of Graham's single frequency approximation, as the super-Lorentzian lineshape diverges at ``ω_0 → 0``.

G_superLorentzian = interpolate_greens_function(greens_superlorentzian, 0, maximum(Tʳᶠ)/T₂ˢ)

f_PSD(τ) = quadgk(ct -> (4 / τ / abs(1 - 3 * ct^2) * (exp(- τ^2 / 8 * (1 - 3 * ct^2)^2) - 1) + sqrt(2π) * erf(τ / 2 / sqrt(2) * abs(1 - 3 * ct^2))) / abs(1 - 3 * ct^2), 0.0, 1.0)[1]
Rʳᶠ = @. f_PSD(Tʳᶠ / T₂ˢ) * ω₁^2 * T₂ˢ
z_Graham_spec_superLorentzian = @. (Rʳᶠ * exp(-Tʳᶠ * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ)

z_Sled_superLorentzian = similar(Tʳᶠ)
for i = 1:length(Tʳᶠ)
    param = (ω₁[i], 1, ω₀, R₁, T₂ˢ, G_superLorentzian)
    prob = ODEProblem(apply_hamiltonian_sled!, z₀, (0, Tʳᶠ[i]), param)
    z_Sled_superLorentzian[i] = solve(prob)[end][1]
end

z_gBloch_superLorentzian = similar(Tʳᶠ)
for i = 1:length(Tʳᶠ)
    param = (ω₁[i], 1, ω₀, R₁, T₂ˢ, G_superLorentzian)
    prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, z_fun, (0, Tʳᶠ[i]), param)
    z_gBloch_superLorentzian[i] = solve(prob)[end][1]
end

p = plot(xaxis=:log, legend=:bottomright, xlabel="Tʳᶠ [s]", ylabel="zˢ(Tʳᶠ)")
plot!(p, Tʳᶠ, z_gBloch_superLorentzian, label="generalized Bloch model")
plot!(p, Tʳᶠ, Tʳᶠ .* 0 .+ cos(α), label="cos(α)")
plot!(p, Tʳᶠ, z_Sled_superLorentzian, label="Sled's model")
plot!(p, Tʳᶠ, z_Graham_spec_superLorentzian, label="Graham's spectral model")
#md Main.HTMLPlot(p) #hide

# This simulation reveals the most pronounced deviations of the generalized Bloch model from established models due to the slower decay of the super-Lorentzian Green's function.

# ### Error Analysis
# Assuming a super-Lorentzian lineshape, we quantify the deviations of Sled's model from the generalized Bloch model:

Tʳᶠᵢ = 1e-3 # s
ω₁ᵢ = α / Tʳᶠᵢ # rad/s
param = (ω₁ᵢ, 1, ω₀, R₁, T₂ˢ, G_superLorentzian)

prob = ODEProblem(apply_hamiltonian_sled!, z₀, (0, Tʳᶠᵢ), param)
z_Sled_superLorentzian_i = solve(prob)[end][1]

prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, z_fun, (0, Tʳᶠᵢ), param)
z_gBloch_superLorentzian_i = solve(prob)[end][1]

z_Sled_superLorentzian_i - z_gBloch_superLorentzian_i

# For ``T_{\text{RF}} = 1``ms, the deviations are small compared to the thermal equilibrium magnetization ``z^s_0 = 1``, but with ``T_{\text{RF}} = 0.1``ms, this deviation becomes sizable:

Tʳᶠᵢ = 1e-4 # s
ω₁ᵢ = α / Tʳᶠᵢ # rad/s
param = (ω₁ᵢ, 1, ω₀, R₁, T₂ˢ, G_superLorentzian)

prob = ODEProblem(apply_hamiltonian_sled!, z₀, (0, Tʳᶠᵢ), param)
z_Sled_superLorentzian_i = solve(prob)[end][1]

prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, z_fun, (0, Tʳᶠᵢ), param)
z_gBloch_superLorentzian_i = solve(prob)[end][1]

z_Sled_superLorentzian_i - z_gBloch_superLorentzian_i

#src ###############################################################
#src # export data for plotting
#src ###############################################################
using Printf #src
io = open(expanduser(string("~/Documents/Paper/2021_MT_IDE/Figures/Pulse_Response_FA", round(Int, α / pi * 180), "deg.txt")), "w") #src
write(io, "TRF_s z_Bloch z_gBloch_Lorentzian z_gBloch_Gaussian z_gBloch_superLorentzian z_Sled_Lorentzian z_Sled_Gaussian z_Sled_superLorentzian z_Graham_spec_Lorentzian z_Graham_spec_Gaussian z_Graham_spec_superLorentzian z_Graham_SF_approx_Lorentzian z_Graham_SF_approx_Gaussian \n") #src
for i = 1:length(Tʳᶠ) #src
    write(io, @sprintf("%1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e %1.3e \n", #src
    Tʳᶠ[i],  #src
    z_Bloch[i],  #src
    z_gBloch_Lorentzian[i],  #src
    z_gBloch_Gaussian[i],  #src
    z_gBloch_superLorentzian[i],  #src
    z_Sled_Lorentzian[i],  #src
    z_Sled_Gaussian[i],  #src
    z_Sled_superLorentzian[i],  #src
    z_Graham_spec_Lorentzian[i],  #src
    z_Graham_spec_Gaussian[i],  #src
    z_Graham_spec_superLorentzian[i],  #src
    z_Graham_SF_approx_Lorentzian[i],  #src
    z_Graham_SF_approx_Gaussian[i] #src
    )) #src
end #src
close(io) #src


text:
using Random, Statistics
Random.seed!(0)

trueVar, trueStd = 1/12, sqrt(1/12)

function estVar(n)
    sample = rand(n)
    sum((sample .- 0.5).^2)/n
end

N = 10^7
for n in 5:5:30
    biasVar = mean([estVar(n) for _ in 1:N]) - trueVar
    biasStd = mean([sqrt(estVar(n)) for _ in 1:N]) - trueStd
    println("n = ",n, " Var bias: ", round(biasVar, digits=5),
           	    "\t Std bias: ", round(biasStd, digits=5))
end

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Tmatrix
using Zygote

rx, rz = 1e-6, 1e-6
n_max = 1
k1_r = 1e7;
k1_i = 0.0;
k2_r = 1.5e7;
k2_i = 1e3;

k1 = Complex(k1_r, k1_i)
k2 = Complex(k2_r, k2_i)

n_θ_points = 10
θ_array = collect(LinRange(1e-6, pi, n_θ_points))
r_array = rx * ones(size(θ_array))
ϕ_array = zeros(size(θ_array))
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

# testing out T-matrix calculation and automatic differentiation using <<T_matrix_SeparateRealImag_arbitrary_mesh>>
@time T_sph = calculate_Tmatrix_for_spheroid_SeparateRealImag(
    rx,
    rz,
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i;
    n_θ_points = n_θ_points,
    rotationally_symmetric = true,
    symmetric_about_plane_perpendicular_z = false,
);
println("T-matrix using spheroid <<calculate_Tmatrix_for_spheroid_SeparateRealImag>>")
display(T_sph)

@time T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
    BigFloat_precision,
);
println("T-matrix using an arbitrary mesh <<T_matrix_SeparateRealImag_arbitrary_mesh>>")
display(T)

@time ∂T = Zygote.jacobian(
    Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh,
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
    BigFloat_precision,
)
println("Jacobian for <<T_matrix_SeparateRealImag_arbitrary_mesh>> calcuated by Zygote:")
display(∂T)

# testing gradient descent optimization for <<sin>>
function target_to_be_minimized(t)
    return sin(t)
end

function ∂target_to_be_minimized(t)
    return Zygote.gradient(target_to_be_minimized, t)[1]
end

"""
learning_rate = 0.01
loss_here = 1e6
t_0 = 0.3
∂loss = 0
while loss_here > 0.01
    t_0 -= learning_rate * ∂loss
    loss_here = target_to_be_minimized(t_0)
    ∂loss = ∂target_to_be_minimized(t_0)
    println("loss_here = loss_here, ∂loss = ∂loss")
end
"""

# testing gradient descent optimization for <<T_matrix_SeparateRealImag_arbitrary_mesh>>
# TODO: use an optimization package instead, so that it has better stop criteria and tricks for faster convergence
rx, rz = 1e-6, 1e-6
n_max = 2
k1_r = 1e7;
k1_i = 0.0;
k2_r = 1.5e7;
k2_i = 1e3;

n_θ_points = 20
θ_array = collect(LinRange(1e-6, pi, n_θ_points))
r_array = rx * ones(size(θ_array))
ϕ_array = zeros(size(θ_array))
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

function objective_function_normal_mesh_density(r_array)
    # println("I am <<objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")
    ϕ_array = zeros(size(θ_array))
    T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max,
        k1_r,
        k1_i,
        k2_r,
        k2_i,
        (collect(r_array)),
        θ_array,
        ϕ_array,
        rotationally_symmetric,
        symmetric_about_plane_perpendicular_z,
        BigFloat_precision,
    )
    # return abs(T[1]) # if we need to minimize the first element of T-matrix

    # return sum(abs.(T)) # if we need to minimize the sum of all real and imag

    # return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, Complex(k1_r, k1_i)) # if we need to minimize scattering cross section

    return 1 * Tmatrix.get_OrentationAv_emissivity_from_Tmatrix(
        T,
        Complex(k1_r, k1_i),
        Tmatrix.calculate_surface_area_of_axisymmetric_particle(r_array, θ_array),
    )
end

function objective_function_double_mesh_density(r_array_input)
    r_θ_array = double_mesh_density(r_array_input, θ_array)
    r_array_double_density = r_θ_array[:, 1]
    θ_array_double_density = r_θ_array[:, 2]
    ϕ_array_double_density = zeros(size(θ_array_double_density))

    # println("I am <<objective_function>>, and the type of r_array_double_density is $(typeof(r_array_double_density))")
    println("r_array_input = $r_array_input")
    println("r_array_double_density = $r_array_double_density")

    T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max,
        k1_r,
        k1_i,
        k2_r,
        k2_i,
        (collect(r_array_double_density)),
        θ_array_double_density,
        ϕ_array_double_density,
        rotationally_symmetric,
        symmetric_about_plane_perpendicular_z,
        BigFloat_precision,
    )
    # return abs(T[1]) # if we need to minimize the first element of T-matrix

    # return sum(abs.(T)) # if we need to minimize the sum of all real and imag

    # return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, Complex(k1_r, k1_i)) # if we need to minimize scattering cross section

    return 1 * Tmatrix.get_OrentationAv_emissivity_from_Tmatrix(
        T,
        Complex(k1_r, k1_i),
        Tmatrix.calculate_surface_area_of_axisymmetric_particle(
            r_array_double_density,
            θ_array_double_density,
        ),
    )
end

# objective_function = objective_function_normal_mesh_density
objective_function = objective_function_double_mesh_density

function ∂objective_function(r_array)
    # println("I am <<∂objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")
    return Zygote.gradient(objective_function, (collect(r_array)))[1]
end

using Plots
n_max = 2
n_θ_points = 20
θ_array = collect(LinRange(1e-6, pi, n_θ_points))
r_array_initial, _ = Tmatrix.ellipsoid(0.1e-7, 0.12e-7, θ_array)
# r_array_initial = 0.1e-7 * ones(size(θ_array))
r_array = r_array_initial
ϕ_array = zeros(size(θ_array))
xyz_initial = vcat(
    Tmatrix.convert_coordinates_Sph2Cart.(
        r_array_initial,
        θ_array,
        zeros(size(r_array_initial)),
    )...,
)
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

learning_rate = 0.5e-18
# learning_rate = 1e-2
global loss_here = -1e6
global ∂loss = 0
println("Starting optimization for T-matrix")
global n_iteration = 0
loss_array = []
loss_old = 9999999999999
while (n_iteration < 50) # (loss_here < 0)
    n_iteration += 1
    r_array = r_array .- learning_rate .* ∂loss

    loss_here = objective_function(r_array)
    ∂loss = ∂objective_function(r_array)
    append!(loss_array, loss_here)

    println()
    println("iteration #$n_iteration: loss_here = $loss_here, ∂loss = $∂loss")
    println("r_array = $r_array")

    xyz = vcat(
        Tmatrix.convert_coordinates_Sph2Cart.(r_array, θ_array, zeros(size(r_array)))...,
    )
    p1 = plot(
        xyz_initial[:, 1],
        xyz_initial[:, 3],
        aspect_ratio = :equal,
        label = "initial particle",
    )
    p1 = plot!(
        xyz[:, 1],
        xyz[:, 3],
        aspect_ratio = :equal,
        label = "particle after $n_iteration iterations",
        title = "scattering cross section = $loss_here m",
    )
    p2 = plot(
        1:length(loss_array),
        loss_array,
        xlabel = "iteration #",
        ylabel = "scattering cross section (m)",
    )
    fig = plot(p1, p2, layout = (1, 2), size = (1200, 800))
    savefig(
        fig,
        "cache/iteration_particle_plots/maximizing_emissivity_07_12_2021__/particle_geom_iteration_$n_iteration.png",
    )
    # if (loss_here - loss_old) < 1e-6; break; end
    loss_old = loss_here
end


text:
include("../../../util/src/util_module.jl")

module gridded
using ..util

include("general.jl")
include("curve.jl")
include("curve_2d.jl")

end # curve

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
export inverse_hasmethod
"""
    inverse_hasmethod(f::Function)
Returns the types that a function hasmethod for.
```julia
inverse_hasmethod(string)
```
"""
function inverse_hasmethod(f::Function)
    allmethods = methods(f).ms
    out = Type[]
    for method in allmethods
        t = typereturn(method.sig)
        if !(typeof(t) in [NoArg])
            push!(out, t)
        end
    end
    return out
end

struct NoArg end

function typereturn(sig)
    sigparams = sig.parameters
    sigparams_length = length(sigparams)
    if sigparams_length == 1
        return NoArg()
    else
        return sigparams[end]
    end
end


text:
using DataVisualization: log_scaler, logistic_scaler, standardizer
using DataVisualization: min_max_scaler, max_abs_scaler
using DataVisualization: binner, quantile_binner, uniform_binner

@testset "Scaler test" begin
    isf = logistic_scaler([1.0])
    @test 15 |> isf |> inv(isf) ≈ 15
    isf = log_scaler([1.0])
    @test 15 |> isf |> inv(isf) ≈ 15
    isf = standardizer([1, 10])
    @test 15 |> isf |> inv(isf) ≈ 15
    isf = min_max_scaler([1, 10])
    @test 15 |> isf |> inv(isf) ≈ 15
    isf = max_abs_scaler([1, 10])
    @test 15 |> isf |> inv(isf) ≈ 15
end

@testset "Binning test" begin
    column = [0, 11, 23, 34, 50, 75, 98]
    bins = 0:20:100
    result = [1, 1, 2, 2, 3, 4, 5]
    @test binner(column, bins) ≈ result
    result = [1, 1, 2, 2, 3, 4, 6]
    @test uniform_binner(column, 6) ≈ result
    result = [1, 1, 2, 3, 4, 5, 6]
    @test quantile_binner(column, 6) ≈ result
end


text:
using Printf

function getlgc(r::Integer, a::Integer, c::Integer, m::Integer, sh::Integer)
    state = r
    return function lgcrand()
        state = mod(a * state + c, m)
        return state >> sh
    end
end

seed, nrep = 0, 10
bsdrand = getlgc(seed, 1103515245, 12345, 2 ^ 31, 0)

println("The first $nrep results for a BSD rand seeded with $seed:")
for _ in 1:nrep
    @printf("%14d\n", bsdrand())
end

msrand = getlgc(seed, 214013, 2531011, 2 ^ 31, 16)

println("\nThe first $nrep results for a M\$ rand seeded with $seed:")
for _ in 1:nrep
    @printf("%14d\n", msrand())
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using MyGA
using Test
using Combinatorics
using DataFrames

correct_bag_data_type = ([1 2 3],[1,2,3],[1;2;3],(1,2,3),
                         [1. 2. 3.],[1.,2.,3.],[1.;2.;3.],(1.,2.,3.))
correct_bag_test = combinations(correct_bag_data_type,2)
correct_bag_test_output = [[1.;2.;3.] [1.;2.;3.]]

non_correct_bag_data_type = ([1 2 3],[1 2 3 4],(1,2,3,4,5))
non_correct_bag_test = combinations(non_correct_bag_data_type,2)

@testset "Turtle.jl" begin
    @testset "create a correct bag" begin
        for datain in correct_bag_test
            @test begin
                mybag = Bags(datain[1],datain[2])
                mybag.vertexes == correct_bag_test_output
            end
        end
    end

    @testset "create a non correct bag" begin
        for datain in non_correct_bag_test
            @test_throws DimensionMismatch Bags(datain[1],datain[2])
        end
    end

end

# desition_tree_data_letters = DataFrame(letter = ["a","b", "a"], number = [0,-1,101], class = ["good","bad","good"])
# desition_tree_data_letters_features = (:letter,:number)
# desition_tree_data_letters_categorycol = :class
#
# desition_tree_data_box1 = DataFrame(x = [0,1,0,1], y = [0,0,1,1], out = [0,0,1,1])
# desition_tree_data_box2 = DataFrame(x = [0,-1,1,0,0], y = [0,0,0,-1,1], out = [0,1,1,1,1])
# desition_tree_data_box_features = (:x,:y)
# desition_tree_data_box_categorycol = :out
#
# desition_tree_data_weather = DataFrame( outlook = ["sunny","sunny","overcast","rain","rain","rain","overcast","sunny","sunny","rain","sunny","overcast","overcast","rain"],
#                                         temperature = ["hot","hot","hot","mild","cool","cool","cool","mild","cool","mild","mild","mild","hot","mild"],
#                                         humidity = ["high","high","high","high","normal","normal","normal","high","normal","normal","normal","high","normal","high"],
#                                         windy = [0,1,0,0,0,1,1,0,0,0,1,1,0,1],
#                                         class = [0,0,1,1,1,0,1,0,1,1,1,1,1,0])
# desition_tree_data_weather_features = (:outlook,:temperature,:humidity,:windy)
# desition_tree_data_weather_categorycol = :class
#
# desition_tree_data_in = (   (desition_tree_data_letters,desition_tree_data_letters_features,desition_tree_data_letters_categorycol),
#                             (desition_tree_data_box1,desition_tree_data_box_features,desition_tree_data_box_categorycol),
#                             (desition_tree_data_box2,desition_tree_data_box_features,desition_tree_data_box_categorycol),
#                             (desition_tree_data_weather,desition_tree_data_weather_features,desition_tree_data_weather_categorycol))
#
# desition_tree_data_res =(   Dict("letter" => Dict("b" => "bad","a" => "good")),
#                             Dict("y" => Dict("1" => "1","0" => "0")),
#                             Dict("y" => Dict("1" => "1","0" => "0")),
#                             Dict("outlook" => Dict(
#                                                     "sunny" => Dict("humidity" => Dict(
#                                                                                "normal" => 1,
#                                                                                "high" => 0)),
#                                                     "rain" => Dict("windy" => Dict(
#                                                                                "1" => 0,
#                                                                                "0" => 1)))
#                                                     "overcast" => 1)) # TODO - need to define the last 2 database result
# @testset "DesitionTree.jl" begin
#     @testset "create desition tree" begin
#         for (data,res) in zip(desition_tree_data_in,desition_tree_data_res)
#             @test create_tree(data[1], data[2], data[3]) = res
#     end
# end


text:
using Test
using brainflow
using JSON

@testset "BoardShim constructor" begin
    # for most boards the master_id == the regular id
    params = brainflow.BrainFlowInputParams()
    board_shim = brainflow.BoardShim(brainflow.CYTON_BOARD, params)
    @test board_shim.master_board_id == Integer(brainflow.CYTON_BOARD)

    # for some special boards you need to provide the master_id in the other_info as string...
    params = brainflow.BrainFlowInputParams(other_info = "1")
    board_shim = brainflow.BoardShim(brainflow.STREAMING_BOARD, params)
    @test board_shim.master_board_id == 1

    # test the error if the other_info cannot be parsed as Int32
    params = brainflow.BrainFlowInputParams(other_info = "")
    @test_throws brainflow.BrainFlowError brainflow.BoardShim(brainflow.STREAMING_BOARD, params)
end

@testset "prepare_session error" begin
    # use a board that is not connected
    params = brainflow.BrainFlowInputParams()
    board_shim = brainflow.BoardShim(brainflow.CYTON_BOARD, params)
    @test_throws brainflow.BrainFlowError("Error in prepare_session 13", 13) brainflow.prepare_session(board_shim)
end

@testset "generated functions" begin
    @test brainflow.get_eeg_channels isa Function
end

@testset "data filtering" begin
    @test Int32(brainflow.CONSTANT) == 1
    @test Int32(brainflow.LINEAR) == 2

    @test Int32(brainflow.MEAN) == 0
    @test Int32(brainflow.MEDIAN) == 1
    @test Int32(brainflow.EACH) == 2

    @test Int32(brainflow.BUTTERWORTH) == 0
    @test Int32(brainflow.CHEBYSHEV_TYPE_1) == 1
    @test Int32(brainflow.BESSEL) == 2

    @test Int32(brainflow.NO_WINDOW) == 0
    @test Int32(brainflow.HANNING) == 1
    @test Int32(brainflow.HAMMING) == 2
    @test Int32(brainflow.BLACKMAN_HARRIS) == 3
end

@testset "model params" begin
    params = BrainFlowModelParams(metric = "concentration")
    @test params.metric == brainflow.Concentration()

    params = BrainFlowModelParams(metric = "Relaxation", classifier = "KNN")
    @test params.metric == brainflow.Relaxation()
    @test params.classifier == brainflow.Knn()

    json = JSON.json(BrainFlowModelParams())
    d = JSON.parse(json)
    @test d["classifier"] == 0
    @test d["metric"] == 0
    @test d["file"] == ""
    @test d["other_info"] == ""
end

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
using MarketData
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("TimeSeries readwrite parses csv file correctly") do

    context("1d values array works") do
        @fact typeof(cl.values) --> Array{Float64,1}
    end

    context("2d values array works") do
        @fact typeof(ohlc.values) --> Array{Float64,2}
    end

    context("timestamp parses to correct type") do
        @fact typeof(cl.timestamp)        --> Vector{Date}
        @fact typeof(datetime1.timestamp) --> Vector{DateTime}
    end

    context("meta field is correctly constructed") do
        @fact AAPL.meta  --> "AAPL"
        @fact BA.meta    --> "BA"
        @fact CAT.meta   --> "CAT"
        @fact DELL.meta  --> "DELL"
        @fact EBAY.meta  --> "EBAY"
        @fact F.meta     --> "F"
        @fact GE.meta    --> "GE"
        @fact mdata.meta --> "Apple"
    end
end

facts("const objects have expected length") do

    context("test objects") do
        @fact length(datetime1) --> 5
        @fact length(op)        --> 500
        @fact length(cl)        --> 500
        @fact length(ohlc)      --> 500
        @fact length(ohlcv)     --> 500
        @fact length(mdata)     --> 500
    end

    context("historical data sets") do
        @fact length(AAPL) --> 8336
        @fact length(BA)   --> 13090
        @fact length(CAT)  --> 13090
        @fact length(DELL) --> 6353
        @fact length(EBAY) --> 3842
        @fact length(F)    --> 10491
        @fact length(GE)   --> 13090
    end
end