text:
include("BorderedArrays.jl")

####################################################################################################
"""
	dt = DotTheta( (x,y) -> dot(x,y) / length(x) )

This parametric type allows to define a new dot product from the one saved in `dt::dot`. More precisely:

	dt(u1, u2, p1::T, p2::T, theta::T) where {T <: Real}

computes, the weigthed dot product ``\\langle (u_1,p_1), (u_2,p_2)\\rangle_\\theta = \\theta \\Re \\langle u_1,u_2\\rangle  +(1-\\theta)p_1p_2`` where ``u_i\\in\\mathbb R^N``. The ``\\Re`` factor is put to ensure a real valued result despite possible complex valued arguments.

	normtheta(u, p::T, theta::T)

Compute, the norm associated to weigthed dot product ``\\langle (u_1,p_1), (u_2,p_2)\\rangle_\\theta``.

!!! info "Info"
    This is used in the pseudo-arclength constraint with the dot product ``\\frac{1}{N} \\langle u_1,u_2\\rangle,\\quad u_i\\in\\mathbb R^N``
"""
struct DotTheta{Tdot}
	dot::Tdot		# defaults to (x,y) -> dot(x,y) / length(x)
end

DotTheta() = DotTheta( (x,y) -> dot(x,y) / length(x))

# Implementation of the dot product associated to DotTheta
function (dt::DotTheta)(u1, u2, p1::T, p2::T, theta::T) where {T <: Real}
	# we restrict the type of the parameters because for complex problems, we still want the parameter to be real
	return real(dt.dot(u1, u2) * theta + p1 * p2 * (one(T) - theta))
end

# Implementation of the norm associated to DotTheta
function (dt::DotTheta)(u, p::T, theta::T) where T
	# we restrict the type of the parameters because for complex problems, we still want the parameter to be real
	return sqrt(dt(u, u, p, p, theta))
end

(dt::DotTheta)(a::BorderedArray{vec, T}, b::BorderedArray{vec, T}, theta::T) where {vec, T} = dt(a.u, b.u, a.p, b.p, theta)

(dt::DotTheta)(a::BorderedArray{vec, T}, theta::T) where {vec, T} = dt(a.u, a.p, theta)

####################################################################################################
# equation of the arc length constraint
function arcLengthEq(dt::DotTheta, u, p, du, dp, theta, ds)
	return dt(u, du, p, dp, theta) - ds
end

####################################################################################################
abstract type AbstractTangentPredictor end
abstract type AbstractSecantPredictor <: AbstractTangentPredictor end

"""
	Natural predictor / corrector
"""
struct NaturalPred <: AbstractSecantPredictor end

"""
	Secant tangent predictor
"""
struct SecantPred <: AbstractSecantPredictor end

"""
	Bordered Tangent predictor
"""
struct BorderedPred <: AbstractTangentPredictor end

function getPredictor!(z_pred::M, z_old::M, tau::M, ds, algo::Talgo) where {T, vectype, M <: BorderedArray{vectype, T}, Talgo <: AbstractTangentPredictor}
	# we perform z_pred = z_old + ds * tau
	copyto!(z_pred, z_old)
	axpy!(ds, tau, z_pred)
end

# generic corrector based on Bordered formulation
function corrector(Fhandle, Jhandle, z_old::M, tau_old::M, z_pred::M,
			ds, theta, contparams, dottheta::DotTheta,
			algo::Talgo, linearalgo = MatrixFreeLBS();
			normC = norm,
			callback = (x, f, J, res, iteration; kwargs...) -> true, kwargs...) where {T, vectype, M<:BorderedArray{vectype, T}, Talgo <: AbstractTangentPredictor}
	return newtonPALC(Fhandle, Jhandle,
			z_old, tau_old, z_pred,
			ds, theta,
			contparams, dottheta; linearbdalgo = linearalgo, normN = normC, callback = callback, kwargs...)
end

# corrector based on natural formulation
function corrector(Fhandle, Jhandle, z_old::M, tau_old::M, z_pred::M,
			ds, theta, contparams, dottheta::DotTheta,
			algo::NaturalPred, linearalgo = MatrixFreeLBS();
			normC = norm,
			callback = (x, f, J, res, iteration; kwargs...) -> true, kwargs...) where {T, vectype, M <: BorderedArray{vectype, T}}
	res = newton(u -> Fhandle(u, z_pred.p),
				 u -> Jhandle(u, z_pred.p),
				 z_pred.u, contparams.newtonOptions;
				 normN = normC, callback = callback, kwargs...)
	return BorderedArray(res[1], z_pred.p), res[2], res[3], res[4]
end

# tangent computation using Natural / Secant predictor
# tau_new is the tangent prediction
function getTangent!(tau_new::M, z_new::M, z_old::M, tau_old::M, F, J,
	ds, theta, contparams, normtheta::DotTheta,
	algo::Talgo, verbosity, linearalgo) where {T, vectype, M <: BorderedArray{vectype, T}, Talgo <: AbstractSecantPredictor}
	(verbosity > 0) && println("--> predictor = ", algo)
	# secant predictor: tau = z_new - z_old; tau *= sign(ds) / normtheta(tau)
	copyto!(tau_new, z_new)
	minus!(tau_new, z_old)
	if algo isa SecantPred
		α = sign(ds) / normtheta(tau_new, theta)
	else
		α = sign(ds) / abs(tau_new.p)
	end
	rmul!(tau_new, α)
end

# tangent computation using Bordered system
# tau_new is the tangent prediction
function getTangent!(tau_new::M, z_new::M, z_old::M, tau_old::M, F, J,
	ds, theta, contparams, dottheta::DotTheta,
	algo::BorderedPred, verbosity, linearbdalgo) where {T, vectype, M <: BorderedArray{vectype, T}}
	(verbosity > 0) && println("--> predictor = Bordered")
	# tangent predictor
	ϵ = contparams.finDiffEps
	# dFdl = (F(z_old.u, z_old.p + ϵ) - F(z_old.u, z_old.p)) / ϵ
	dFdl = F(z_old.u, z_old.p + ϵ)
	minus!(dFdl, F(z_old.u, z_old.p))
	rmul!(dFdl, 1/ϵ)

	# tau = getTangent(J(z_old.u, z_old.p), dFdl, tau_old, theta, contparams.newtonOptions.linsolve)
	tau_normed = copyto!(similar(tau_old), tau_old) #copy(tau_old)
	rmul!(tau_normed, theta / length(tau_old.u), 1 - theta)
	# extract tangent as solution of bordered linear system, using zero(z_old.u)
	tauu, taup, flag, it = linearbdalgo( J(z_old.u, z_old.p), dFdl,
			tau_normed, rmul!(similar(z_old.u), false), T(1), theta)

	# the new tangent vector must preserve the direction along the curve
	α = T(1) / dottheta(tauu, tau_old.u, taup, tau_old.p, theta)

	# tau_new = α * tau
	copyto!(tau_new.u, tauu)
	tau_new.p = taup
	rmul!(tau_new, α)
end
################################################################################################
function arcLengthScaling(theta, contparams, tau::M, verbosity) where {M <: BorderedArray}
	# the arclength scaling algorithm is based on Salinger, Andrew G, Nawaf M Bou-Rabee, Elizabeth A Burroughs, Roger P Pawlowski, Richard B Lehoucq, Louis Romero, and Edward D Wilkes. “LOCA 1.0 Library of Continuation Algorithms: Theory and Implementation Manual,” March 1, 2002. https://doi.org/10.2172/800778.

	thetanew = theta
	g = abs(tau.p * theta)
	(verbosity > 0) && print("Theta changes from $(theta) to ")
	if (g > contparams.gMax)
		thetanew = contparams.gGoal / tau.p * sqrt( abs(1.0 - g^2) / abs(1.0 - tau.p^2) )
		if (thetanew < contparams.thetaMin)
		  thetanew = contparams.thetaMin;
		end
	end
	(verbosity > 0) && print("$(thetanew)\n")
	return thetanew
end
################################################################################################
function stepSizeControl(ds, theta, contparams, converged::Bool, it_newton_number::Int64, tau::M, verbosity) where {T, vectype, M<:BorderedArray{vectype, T}}
	if converged == false
		if  abs(ds) <= contparams.dsmin
			(verbosity > 0) && printstyled("*"^80*"\nFailure to converge with given tolerances\n"*"*"^80, color=:red)
			# we stop the continuation
			return ds, theta, true
		end
		dsnew = sign(ds) * max(abs(ds) / 2, contparams.dsmin);
		(verbosity > 0) && printstyled("Halving continuation step, ds=$(dsnew)\n", color=:red)
	else
		# control to have the same number of Newton iterations
		Nmax = contparams.newtonOptions.maxIter
		factor = (Nmax - it_newton_number) / Nmax
		dsnew = ds * (1 + contparams.a * factor^2)
		(verbosity > 0) && @show 1 + contparams.a * factor^2
	end

	# control step to stay between bounds
	if abs(dsnew) < contparams.dsmin
		dsnew = sign(dsnew) * contparams.dsmin
	end

	if abs(dsnew) > contparams.dsmax
		dsnew = sign(dsnew) * contparams.dsmax
	end

	if contparams.doArcLengthScaling
		thetanew = arcLengthScaling(theta, contparams, tau, verbosity)
	else
		thetanew = theta
	end
	@assert abs(dsnew) >= contparams.dsmin "Error with ds value"
	return dsnew, thetanew, false
end
####################################################################################################
"""
This is the classical matrix-free Newton Solver used to solve `F(x, p) = 0` together
with the scalar condition `n(x, p) = (x - x0) * xp + (p - p0) * lp - n0`
"""
function newtonPALC(F, Jh,
						z0::BorderedArray{vectype, T},
						tau0::BorderedArray{vectype, T},
						z_pred::BorderedArray{vectype, T},
						ds::T, theta::T,
						contparams::ContinuationPar{T},
						dottheta::DotTheta;
						linearbdalgo = BorderingBLS(),
						normN = norm,
						callback = (x, f, J, res, iteration, optionsN; kwargs...) ->  true, kwargs...) where {T, vectype}
	# Extract parameters
	newtonOpts = contparams.newtonOptions
	@unpack tol, maxIter, verbose, alpha, almin, linesearch = newtonOpts
	@unpack finDiffEps = contparams

	N = (x, p) -> arcLengthEq(dottheta, minus(x, z0.u), p - z0.p, tau0.u, tau0.p, theta, ds)
	normAC = (resf, resn) -> max(normN(resf), abs(resn))

	# Initialise iterations
	x = copyto!(similar(z_pred.u), z_pred.u) # copy(z_pred.u)
	p = z_pred.p
	x_pred = copyto!(similar(x), x) # copy(x)

	# Initialise residuals
	res_f = F(x, p);  res_n = N(x, p)

	dX   = copyto!(similar(res_f), res_f) # copy(res_f)
	dp   = T(0)
	# dFdp = (F(x, p + finDiffEps) - res_f) / finDiffEps
	dFdp = copyto!(similar(res_f), F(x, p + finDiffEps)) # copy(F(x, p + finDiffEps))
	minus!(dFdp, res_f)	# dFdp = dFdp - res_f
	rmul!(dFdp, T(1) / finDiffEps)

	res     = normAC(res_f, res_n)
	resHist = [res]
	it = 0

	# Displaying results
	verbose && displayIteration(it, 1, res)
	step_ok = true

	# Main loop
	while (res > tol) & (it < maxIter) & step_ok
		# copyto!(dFdp, (F(x, p + epsi) - F(x, p)) / epsi)
		copyto!(dFdp, F(x, p + finDiffEps)); minus!(dFdp, res_f); rmul!(dFdp, T(1) / finDiffEps)

		J = Jh(x, p)
		u, up, flag, liniter = linearbdalgo(J, dFdp, tau0, res_f, res_n, theta)

		if linesearch
			step_ok = false
			while !step_ok & (alpha > almin)
				# x_pred = x - alpha * u
				copyto!(x_pred, x)
				axpy!(-alpha, u, x_pred)

				p_pred = p - alpha * up
				copyto!(res_f, F(x_pred, p_pred))

				res_n  = N(x_pred, p_pred)
				res = normAC(res_f, res_n)

				if res < resHist[end]
					if (res < resHist[end] / 2) & (alpha < 1)
						alpha *= 2
					end
					step_ok = true
					copyto!(x, x_pred)
					p  = p_pred
				else
					alpha /= 2
				end
			end
		else
			minus!(x, u) 	# x .= x .- u
			p = p - up

			copyto!(res_f, F(x, p))

			res_n  = N(x, p)
			res = normAC(res_f, res_n)
		end

		push!(resHist, res)
		it += 1

		callback(x, res_f, J, res, it, contparams; kwargs...) == false && (it = maxIter)
		verbose && displayIteration(it, 1, res, liniter)

	end
	return BorderedArray(x, p), resHist, resHist[end] < tol, it
end
####################################################################################################


text:
"""
Recursive Brim
"""
function recursive_brim!(M::Modular)
  # -2- R and T sub-matrices
  n_rows, n_cols = size(M.A)
  m = sum(M.A)
  R = M.S[[1:n_cols],:]
  T = M.S[[(n_cols+1):end],:]
  B = M.A - kron(sum(M.A, 1), sum(M.A, 2)) ./ m

  # -3- Recursive BRIM
  old_Q = 0.0
  new_Q = 0.000001
  while new_Q > old_Q
    t_tilde = B'*T
    bestPart!(R, t_tilde)
    r_tilde = B*R
    bestPart!(T, r_tilde)
    old_Q = new_Q
    new_Q = maximum(map((x) -> x/m, [sum((B'*T) .* R), sum((B*R) .* T)]))
  end
  M.S[[1:n_cols],:] = R
  M.S[[(n_cols+1):end],:] = T
  no_empty_modules!(M)
  return M
end

"""
Returns the best partition (used within recursive_brim)
"""
function bestPart!(Target::Array{Int64, 2}, Scores::Array{Float64, 2})
  best_scores = maximum(Scores, 2)
  @inbounds for i = 1:size(Target)[1]
    for c = 1:size(Scores)[2]
      if Scores[i,c] == best_scores[i]
        Target[i,:] = zeros([1:size(Target)[2]])
        Target[i,c] = 1
        break
      else
        Target[i,c] = 0
      end
    end
  end
end

ext: jl


text:
function do_env_test(min::VersionNumber, max::VersionNumber)
    if isdefined(Pkg, :dependencies)
        deps = Pkg.dependencies()
        uuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a")
        @test haskey(deps, uuid) && min <= deps[uuid].version <= max
    else
        deps = Pkg.installed()
        @test haskey(deps, "Example") && min <= deps["Example"] <= max
    end
end

proj = Base.active_project()
try
    project!"""
    [deps]
    Example = "7876af07-990d-54b4-ab0e-23690620f79a"

    [compat]
    Example = "0.4"
    """
    do_env_test(v"0.4", v"0.4.99")

    manifest!"""
    # This file is machine-generated - editing it directly is not advised

    [[Example]]
    git-tree-sha1 = "46e44e869b4d90b96bd8ed1fdcf32244fddfb6cc"
    uuid = "7876af07-990d-54b4-ab0e-23690620f79a"
    version = "0.5.3"
    """
    do_env_test(v"0.5.3", v"0.5.3")
finally
    proj === nothing ? Pkg.activate() : Pkg.activate(proj)
end


text:
###################################
### Last Update: 29 Oct 20
###################################
### ''Simulation of Optimal Proposal''
### Code follows method/notation in Dietrich & Newsam (96)
###################################

using Distributions
using LinearAlgebra
using FFTW
using Random

####################################################################
####################################################################

Random.seed!(123)

#### Parameters for Covariance Function
struct f_th
    s::Float64
    l::Float64
end

#### Parameters for Grid
struct f_gr
    p::Int ## Number of Grid Points Along x-axis is (p+1)
    q::Int ## Number of Grid Points Along y-axis is (q+1)
    P::Float64
    Q::Float64 ## 2D-Field is rectangle [0,P]x[0,Q]
end

#### Parameters for Stations
struct f_st
  st_ij::Array{Int64,2} ## the position of the stations on the grid
                        ## i.e station positions determined using 2 co-ordinates
   st_p::Array{Int64,1} ## the no of grid point occupied by each station
                        ## i.e station positions determined using 1 co-ordinate
end

####################################################################
#### Example #######################################################

th = f_th(1.0, 1.0)

gr = f_gr(20, 20, 40.0, 40.0)

###################################################
#### Example ######################################

p = gr.p; q = gr.q; P = gr.P; Q = gr.Q

Dx = P/p; Dy = Q/q ## Grid Steps Along x and y Directions

###################################################
#### Example ######################################

#### In real case Stations will be read from a file containing co-ordinates;
#### here, I delect position of n2 Stations at random from all grid nodes

n2 = 3; n1 = (p+1)*(q+1)

grid_F = fill(0,(n1,2))

#### the loop contructs 2-dim co-ordinates for all grid points
for j in 0:q

  for i in 0:p

    grid_F[j*(p+1)+i+1,:] = [i, j]

  end

end
####

st_p = sample(1:n1,n2,replace=false)

st_ij = grid_F[st_p,:]

st = f_st(st_ij, st_p)

###################################################
#### Example ######################################

#### I simulate some arbitrary values to be used as observations of height
#### I give an example value to the standard deviation Sobs;
#### I only consider a single time instance

Yobs_t = rand(Normal(0,1), n2)
#### the ordering in the observations is in tandem with the ordering
#### of the stations in variable st

Sobs = 0.01

###################################################
###################################################

#### Example Specification of Covariance Function
function R(x::Float64, y::Float64, th)

  Cov = (th.s^2)*exp(-(abs(x)+abs(y))/(2.0*th.l))

  return(Cov)

end

###################################################
###################################################

#### Specification of 'Extented' Covariance Function
#### following Method/Notation in Dietrich & Newsam (96)
function bar_R(x::Float64, y::Float64, th, gr)

  P = gr.P; Q = gr.Q

  bar_Cov = 0.0

  if ( (0<=x) & (x<=P) & (0<=y) & (y<=Q) )

    bar_Cov = R(x, y, th)

  elseif ( (P<=x) & (x<=2*P) & (0<=y) & (y<=Q) )

    bar_Cov = R(2*P-x, y, th)

  elseif ( (0<=x) & (x<=P) & (Q<=y) & (y<=2*Q) )

    bar_Cov = R(x, 2*Q-y, th)

  elseif ( (P<=x) & (x<=2*P) & (Q<=y) & (y<=2*Q) )

    bar_Cov = R(2*P-x, 2*Q-y, th)

  end

  return(bar_Cov)

end

###################################################
###################################################

#### Obtain (4pq)x(1) Vector bar_rho
#### following Method/Notation in Dietrich & Newsam (96)
function f_bar_rho(th, gr)

  p = gr.p; q = gr.q; P = gr.P; Q = gr.Q

  Dx = P/p; Dy = Q/q

  bar_rho = zeros(4*p*q)

  #########

  for j in 0:(2*q-1)

    for i in 0:(2*p-1)

      bar_rho[(j*2*p)+(i+1)] = bar_R(i*Dx, j*Dy, th, gr)

    end

  end

  #########

  return(bar_rho)

end

###################################################
#### Example ######################################

####bar_rho = f_bar_rho(th, gr)

###################################################
###################################################

#### Construct W Operation on (4pq)x1 Vector V
#### W is 2D Discrete Fourier Transform
#### following Method/Notation in Dietrich & Newsam (96);
#### V is (4pq)x1; result is (4pq)x1
function f_W(V, gr)

  p = gr.p; q = gr.q

  WV = zeros(4*p*q)

  #### Re-arrange (4pq)x(1) Vector V into an (2p)x(2q) matrix
  #### by using first 2p-elements of V as 1st column,
  #### next 2p-elements as 2nd column, and so on.
  bar_E = reshape(V, (2*p,2*q))

  #### Use Julia's 2D-FFT tranform
  E = sqrt(4*p*q)*ifft(bar_E)
  #### E is in general a Complex-Value (2p)x(2q) Matrix

  #### Re-arrange (2p)x(2q) matrix E
  #### into an (4pq)x1 matrix, by using 1st column of E,
  #### then 2nd column of E, and so on.
  WV = reshape(E,(4*p*q))

  return(WV)

end

###################################################
#### Example ######################################

####[bar_rho f_W(bar_rho, gr)]

###################################################
###################################################

#### Construct WH Operation on (4pq)x1 Vector V
#### WH is 2D Inverse Discrete Fourier Transform
#### following Method/Notation in Dietrich & Newsam (96);
#### V is (4pq)x1; result is (4pq)x1
function f_WH(V, gr)

  p = gr.p; q = gr.q

  WHV = zeros(4*p*q)

  #### Re-arrange (4pq)x(1) Vector V into an (2p)x(2q) matrix
  #### by using first 2p-elements of V as 1st column,
  #### next 2p-elements as 2nd column, and so on.
  E = reshape(V, (2*p,2*q))

  #### Use Julia's 2D-IFFT tranform
  bar_E = (1.0/sqrt(4*p*q))*fft(E)
  #### bar_E is in general a Complex-Value (2p)x(2q) Matrix

  #### Re-arrange (2p)x(2q) matrix bar_E
  #### into an (4pq)x1 matrix, by using 1st column of bar_E,
  #### then 2nd column of bar_E, and so on.
  WHV = reshape(bar_E,(4*p*q))

  return(WHV)

end

###################################################
#### Example ######################################

####V = bar_rho; WV = f_W(V,gr); WHWV = f_WH(WV,gr);

#### CHECK!! ######################################

####[V WV WHWV] ## Since WH*W = I, for any V, it should be V = WHWV!!

###################################################
###################################################

###################################################
###################################################
####
#### Construct Covariance Matrix on Original Grid
#### following Method/Notation in Dietrich & Newsam (96)
#### This Function is Needed only for Testing!!
####function f_R11(th, gr)
####
####  p = gr.p; q = gr.q; P = gr.P; Q = gr.Q
####
####  Dx = P/p; Dy = Q/q; n1 = (p+1)*(q+1)
####
####  R11 = zeros(n1, n1)
####
####  for m in 0:q
####
####    for l in 0:p
####
####      for j in 0:q
####
####        for i in 0:p
####
####          d1 = m*(p+1)+l+1; d2 = j*(p+1)+i+1
####
####          R11[d1,d2] = R(abs(l*Dx-i*Dx), abs(m*Dy-j*Dy), th)
####
####        end
####
####      end
####
####    end
####
####  end
####
####  return(R11)
####
####end
####
###################################################
#### Example ######################################
####
####R11 = f_R11(th, gr)
####
#### Construct Covariance Matrix on Extended Grid
#### following Method/Notation in Dietrich & Newsam (96)
#### This Function is Needed only for Testing!!
####function f_bar_R11(th, gr)
####
####  p = gr.p; q = gr.q; P = gr.P; Q = gr.Q
####
####  Dx = P/p; Dy = Q/q; bar_n1 = 4*p*q
####
####  bar_R11 = zeros(bar_n1, bar_n1)
####
####  for m in 0:(2*q-1)
####
####    for l in 0:(2*p-1)
####
####      for j in 0:(2*q-1)
####
####        for i in 0:(2*p-1)
####
####          d1 = m*2*p+l+1; d2 = j*2*p+i+1
####
####          bar_R11[d1,d2] = bar_R(abs(l*Dx-i*Dx), abs(m*Dy-j*Dy), th, gr)
####
####        end
####
####      end
####
####    end
####
####  end
####
####  return(bar_R11)
####
####end
####
###################################################
#### Example ######################################
####
####bar_R11 = f_bar_R11(th, gr)
####
###################################################
###################################################

#### Construct Covariance between Stations and Extended Grid
#### following Method/Notation in Dietrich & Newsam (96);
#### result is matrix of dimnesions n2xbar_n1
function f_bar_R21(th, st, gr)

  st_ij = st.st_ij; n2 = size(st_ij)[1];
  ## No of Stations

  p = gr.p; q = gr.q; P = gr.P; Q = gr.Q

  Dx = P/p; Dy = Q/q

  bar_R21 = zeros(n2, 4*p*q)

  for m in 1:n2

    for j in 0:(2*q-1)

      for i in 0:(2*p-1)

        d2 = j*2*p+i+1

        bar_R21[m,d2] = bar_R(abs(st_ij[m,1]-i)*Dx, abs(st_ij[m,2]-j)*Dy, th, gr)

      end

    end

  end

  return(bar_R21)

end

###################################################
#### Example ######################################

####bar_R21 = f_bar_R21(th, st, gr)

###################################################
###################################################

#### Construct Covariance between Observations
#### following Method/Notation in Dietrich & Newsam (96);
#### result is matrix of dimnesions n2xn2
function f_R22(th, st, Sobs, gr)

  st_ij = st.st_ij; n2 = size(st_ij)[1]
  ## No of Stations

  p = gr.p; q = gr.q; P = gr.P; Q = gr.Q

  Dx = P/p; Dy = Q/q

  ############

  R22 = zeros(n2, n2)

  for i in 1:n2

    xi = st_ij[i,1]
    yi = st_ij[i,2]

    for j in 1:n2

      xj = st_ij[j,1]
      yj = st_ij[j,2]

      R22[i,j] = R(abs(xj-xi)*Dx,abs(yj-yi)*Dy,th)

      if (i==j)

        R22[i,j] = R22[i,j] + Sobs^2

      end

    end

  end

  return(R22)

end

###################################################
#### Example ######################################

####R22 = f_R22(th, st, Sobs, gr)

###################################################
###################################################

#### Construct Covariance between Original Grid & Station Positions
#### following Method/Notation in Dietrich & Newsam (96)
#### result is matrix of dimnesions n1xn2
function f_R12(th, st, gr)

  st_ij = st.st_ij; n2 = size(st_ij)[1]
  ## No of Stations

  p = gr.p; q = gr.q; P = gr.P; Q = gr.Q

  Dx = P/p; Dy = Q/q; n1 = (p+1)*(q+1)

  R12 = zeros(n1, n2)

  for m in 1:n2

    for j in 0:q

      for i in 0:p

        d1 = j*(p+1)+i+1

        R12[d1,m] = R(abs(st_ij[m,1]-i)*Dx, abs(st_ij[m,2]-j)*Dy, th)

      end

    end

  end

  return(R12)

end

###################################################
#### Example ######################################

####R12 = f_R12(th, st, gr)

###################################################
###################################################
####
#### Construct K Matrix
#### following Method/Notation in Dietrich & Newsam (96)
####function f_K(th, st, gr)
####
####  st_ij = st.st_ij; n2 = size(st_ij)[1]
####  ## No of Stations
####
####  p = gr.p; q = gr.q
####
####  bar_rho = f_bar_rho(th, gr)
####
####  Lambda = sqrt(4*p*q)*Diagonal(real.(f_W(bar_rho, gr)))
####  #### THEORY CHECKPOINT:
####  #### All values in Lambda must be POSITIVE;
####  #### OTHERWISE, the METHOD FAILS!!!
####
####  bar_R21 = f_bar_R21(th, st, gr); bar_R12 = bar_R21'
####
####  WHbar_R12 = zeros(Complex, 4*p*q, n2)
####  for i in 1:n2
####
####    WHbar_R12[:,i] = f_WH(bar_R12[:,i], gr)
####
####  end
####
####  KH = Lambda^(-1/2)*WHbar_R12
####
####  K = KH' ## this gives Cnjugate Transpose
####
####  return(K)
####
####end
####
###################################################
### Example #######################################

#### K = f_K(th, st, gr)

##### Check!! #####################################

#### sort(eigen(bar_R11).values)
#### THEORY CHECKPOINT:
#### All eigenvalues of bar_R11 must be POSITIVE;
#### OTHERWISE, the METHOD FAILS!!!

#### A = bar_R21*inv(bar_R11)*(bar_R21)'

#### B = real.(K*K')

#### Check: Must Have A=B!!

####################################

#### bar_RR = [bar_R11 (bar_R21)'; bar_R21 R22]

#### THEORY CHECKPOINT:
#### All eigenvalues of bar_RR must be POSITIVE;
#### OTHERWISE, the METHOD FAILS!!!
#### sort(eigen(bar_RR).values)

#### Lambda = sqrt(4*p*q)*Diagonal(real.(f_W(bar_rho, gr)))

#### [sort(eigen(bar_R11).values) sort(eigen(Lambda).values)]
#### Check: The two colums must be the same !!

###################################################
###################################################
#### V is (n2)x1; result is (n1)x1
function f_A1_T(V, st_p, n1)

  n2 = length(V); A1_T_V = zeros(n1)

  for i in 1:n2

    A1_T_V[st_p[i]] = V[i]

  end

  return(A1_T_V)

end

###################################################
### Example #######################################

####V = rand(n2)

####A1_T_V = f_A1_T(V, st_p, n1)

###################################################
###################################################
#### U is (n1)x1; result is (n2)x1
function f_A1(U, st_p)

  n1 = length(U); n2 = length(st_p)

  A1_U = zeros(n2)

  for i in 1:n2

    A1_U[i] = U[st_p[i]]

  end

  return(A1_U)

end

###################################################
### Example #######################################

####U = rand(n1)

####A1_U = f_A1(U, st_p)

###################################################
###################################################
#### V is (n1)x1; result is (n1)x1
function f_Sigma(V, th, gr)

  p = gr.p; q = gr.q; n1 = (p+1)*(q+1); bar_n1 = 4*p*q

  ############

  bar_V = zeros(bar_n1)

  for j in 0:q

    bar_V[(j*2*p+1):(j*2*p+p+1)] = V[(j*(p+1)+1):((j+1)*(p+1))]

  end

  bar_rho = f_bar_rho(th, gr)

  Lambda = sqrt(bar_n1)*Diagonal(real.(f_W(bar_rho, gr)))
  #### THEORY CHECKPOINT:
  #### All values in Lambda must be POSITIVE;
  #### OTHERWISE, the METHOD FAILS!!!

  bar_Sigma_V = f_W(Lambda*f_WH(bar_V, gr),gr)

  ############

  Proj_Index = fill(0, n1)

  for j in 0:q

    Proj_Index[(j*(p+1)+1):((j+1)*(p+1))] = collect((j*2*p+1):(j*2*p+p+1))

  end

  Sigma_V = bar_Sigma_V[Proj_Index]

  return(real.(Sigma_V))

end

###################################################
### Example #######################################

####V = rand(n1); Sigma_V = f_Sigma(V, th, gr)

####maximum(abs.(R11*V-Sigma_V)) #### difference must be 0

###################################################
###################################################

####################################################################
####################################################################
### ALL FUNCTIONS SO FAR ARE STAND-ALONE ONES;
### THE ONE THAT FOLLOW NEED INTEGRATION WITH TUOMAS/MOSE CODE;
### THE POSITION ASSUMED IN THE CODE THAT FOLLOWS IS THAT WE ARE
### AT TIME t-1, AND HAVE N PARTICLES, EACH OF DIMENSION 3x(n1),
### REPRESENTING THE 3-DIMENSIONAL VECTOR FIELD ON A GRID OF n1 POSITIONS,
### AT TIME t-1.
### THE FUNCTIONS THE FOLLOW APPLY A SINGLE PARTICLE FILTER STEP.
### IN PARTICULAR, THE CODE THAT FOLLOWS DOES THE FOLLOWING:
### - IT SIMULATES THE OPTIMAL PROPOSAL FOR THE HEIGHT FIELD AT TIME t.
### - IT *DOES NOT* SIMULATE THE PROPOSAL FOR THE TWO VELOCITY FIELDS.
### - FOR THE LATTER PURPOSE THE CODE MUST BE COMBINED WITH THE CODE OF
###   OF TUOMAS/MOSE
###################################################################
####################################################################

#### Calculation of the mean for the optimal proposal for the
#### Height Field at time t; it is assumed that INTEGRATION with TUOMAS/MOSE code
#### it is assumed that INTEGRATION with TUOMAS/MOSE code will provide the
#### (N)x(n1) field FH_t which is the push-forward of all N particles at time t-1
#### under the PDE dynamics; FH_t then considers only the (N)x(n1) matrix with
#### the values of the height field. Y_obs_t are the (n2)x1 data at time t.
#### THE RESULT IS A MATRIX (N)x(n1)
function Calculate_Mean(FH_t, th, st, Yobs_t, Sobs, gr)

  ######### Offline Calculations #########

  p = gr.p; q = gr.q; n1 = (p+1)*(q+1)

  st_ij = st.st_ij; st_p = st.st_p

  n2 = size(st_ij)[1]; N = size(FH_t)[1]
  ## No of Stations + No of Particles

  ##########################

  R22 = f_R22(th, st, Sobs, gr)

  invR22 = inv(R22)

  ##########################

  mu20 = Sobs^(-2)*(R22 - Diagonal(fill(Sobs^2,n2)))
  #### (n2)x(n2)

  mu21 = f_Sigma(f_A1_T(invR22*(mu20*Yobs_t), st_p, n1), th, gr)
  #### (n1)x1

  mu22 = f_Sigma((Sobs^(-2))*f_A1_T(Yobs_t, st_p, n1), th, gr)
  #### (n1)x1

  mu2 = mu22 - mu21
  #### (n1)x1
  ##########################

  mu1 = zeros(N, n1)

  mu  = zeros(N, n1)

  ##### the required mean is different for each of the N particles
  for i in 1:N

    mu10 = invR22*f_A1(FH_t[i,:],st_p)
    #### (n2)x1

    mu11 = f_Sigma(f_A1_T(mu10, st_p, n1), th, gr)
    #### (n1)x1

    mu1[i,:] = FH_t[i,:] - mu11
    #### (n1)x1

    mu[i,:] = mu1[i,:] + mu2
    #### (n1)x1

  end

  return(mu)

end

###################################################
### Example #######################################

####N = 10; FH_t = rand(N, n1)

####Calculate_Mean(FH_t, th, st, Yobs_t, Sobs, gr)

###################################################
###################################################
#
# INPUT:
#### FH_t: an (N)x(n1) matrix where N is the number of particles
#### and n1 the size of the grid on the 2D-domain.
#### The i-th row of FH_t (referring to the i-th particle)
#### is the n1-dimensional Height field calculated after solving the PDE,
#### from time instance (t-1) to time t, using the i-th particle at time t-1 as initial condition.
#### th: Parameters in Covariance Function
#### st: Information about Stations
#### Yobs_t: n2-dimensional vector corresponding to the observations from the n2 stations
#### Sobs: Standard Deviation of Observations at Stations
#### gr: Parameters Specifying Regular Grid
# OUTPUT:
#### The numbers of particles N must be EVEN!!
function Sample_Height_Proposal(FH_t, th, st, Yobs_t, Sobs, gr)

  ######### Offline Calculations: #########
  ######### i.e. will be the same for 0->1, 1->2, t-1->t, etc. #########

  st_ij = st.st_ij; st_p = st.st_p

  n2 = size(st_ij)[1]; N = size(FH_t)[1]
  ## No of Stations + No of Particles

  p = gr.p; q = gr.q; P = gr.P; Q = gr.Q

  bar_n1 = 4*p*q; n1 = (p+1)*(q+1)

  Dx = P/p; Dy = Q/q

  bar_rho = f_bar_rho(th, gr)

  Lambda = sqrt(bar_n1)*Diagonal(real.(f_W(bar_rho, gr)))
  #### THEORY CHECKPOINT:
  #### All values in Lambda must be POSITIVE;
  #### OTHERWISE, the METHOD FAILS!!!

  bar_R21 = f_bar_R21(th, st, gr); bar_R12 = bar_R21'

  WHbar_R12 = zeros(Complex, bar_n1, n2)

  for i in 1:n2

    WHbar_R12[:,i] = f_WH(bar_R12[:,i], gr)

  end

  KH = Lambda^(-1/2)*WHbar_R12; K = KH' ## this gives Conjugate Transpose

  R22 = f_R22(th, st, Sobs, gr)

  L = cholesky(real.(R22-K*KH)).L
  #### THEORY CHECKPOINT:
  #### R22-K*KH must have POSITIVE EIGENVALUES,
  #### i.e. all(eigen(R22-K*KH).values)>0,
  #### OTHERWISE, the METHOD FAILS!!!

  R12 = f_R12(th, st, gr);

  R12_invR22 = R12*inv(R22)

  ##muZ1_C_Yobs = R12_invR22*Yobs_t

  Proj_Index = fill(0, n1)

  for j in 0:q

    Proj_Index[(j*(p+1)+1):((j+1)*(p+1))] = collect((j*2*p+1):(j*2*p+p+1))

  end

  ####

  Samples = zeros(N,n1)
  N0 = Int(N/2) ########## N must be even!

  ########################################
  ######### Online Calculations: ##########
  ######### here we need info from the data at time t, or/and a separate ##########
  ######### calculation for each particle. ##########
  means = Calculate_Mean(FH_t, th, st, Yobs_t, Sobs, gr)
  rng = Random.MersenneTwister(123)

  e1 = Vector{ComplexF64}(undef, bar_n1)
  e2 = Vector{ComplexF64}(undef, n2)

  for i in 0:(N0-1)

    @. e1 = complex(randn(rng), randn(rng))
    @. e2 = complex(randn(rng), randn(rng))

    bar_z1 = f_W(Lambda^(1/2)*e1, gr)
    z2 = (K*e1) + (L*e2)

    Rbar_z = real.([bar_z1; z2])
    Ibar_z = imag.([bar_z1; z2])

    Proj1 = [ Rbar_z[Proj_Index]; Rbar_z[(bar_n1+1):end] ]

    Proj2 = [ Ibar_z[Proj_Index]; Ibar_z[(bar_n1+1):end] ]

    ##Z1_C_Yobs_1 = muZ1_C_Yobs + Proj1[1:n1] - R12_invR22*Proj1[(n1+1):end]
    Z11 = Proj1[1:n1] - R12_invR22*Proj1[(n1+1):end]

    ##Z1_C_Yobs_2 = muZ1_C_Yobs + Proj2[1:n1] - R12_invR22*Proj2[(n1+1):end]
    Z12 = Proj2[1:n1] - R12_invR22*Proj2[(n1+1):end]

    ##Samples[2*i+1,:] = means[2*i+1,:] + Z1_C_Yobs_1
    Samples[2*i+1,:] = means[2*i+1,:] + Z11

    ##Samples[2*i+2,:] = means[2*i+2,:] + Z1_C_Yobs_2
    Samples[2*i+2,:] = means[2*i+2,:] + Z12

  end

  return(Samples)

end

####################################################################
####################################################################

###################################################
#### Example ######################################

N = 10; FH_t = rand(N, n1)

Ss = Sample_Height_Proposal(FH_t, th, st, Yobs_t, Sobs, gr)

####

Ss[1,st_p] - Yobs_t #### all the differences should be close to 0 if Sobs<<1

Ss[2,st_p] - Yobs_t #### all the differences should be close to 0 if Sobs<<1

Ss[3,st_p] - Yobs_t #### all the differences should be close to 0 if Sobs<<1

Ss[4,st_p] - Yobs_t #### all the differences should be close to 0 if Sobs<<1

####

Fs = zeros(N, q+1, p+1)

for i in 1:N

  Fs[i,:,:] = reshape(Ss[i,:], (p+1,q+1))'

end

i = 2

#heatmap(1:size(Fs[i,:,:],2),
#           1:size(Fs[i,:,:],1), Fs[i,:,:],
#           c=cgrad([:blue, :white,:red, :yellow]),
#           xlabel="x-axis", ylabel="y-axis",
#           title="")

#####!!!!!! Need to Check if I Combine in the Correct Way
#####!!!!!! the algorithm for the Simulation of the Conditional Field
#####!!!!!! and the Use of the Mean.

###################################################
###################################################


text:
@testset "@add_reactions! helper" begin
    mod = StandardModel()
    A = Metabolite("A")
    B = Metabolite("B")
    C = Metabolite("C")
    add_metabolites!(mod, [A, B, C])

    @add_reactions! mod begin
        "v1", nothing ↔ A
        "v2", nothing ↔ B, -500
        "v3", nothing ↔ C, -500, 500
    end

    rxn = mod.reactions["v1"]
    @test rxn.lb == -1000.0
    @test rxn.ub == 1000.0

    rxn = mod.reactions["v2"]
    @test rxn.lb == -500
    @test rxn.ub == 1000.0

    rxn = mod.reactions["v3"]
    @test rxn.lb == -500
    @test rxn.ub == 500
end


text:
export convolution, convolution!, ∇convolution!
export CUDNN_CONVOLUTION, CUDNN_CROSS_CORRELATION

#=
type ConvolutionDesc{N}
  ptr::Ptr{Void}
  padding::NTuple{N,Int}
  stride::NTuple{N,Int}
end

function ConvolutionDesc{T,N}(::Type{T}, padding::NTuple{N,Int}, stride::NTuple{N,Int};
  mode=CUDNN_CROSS_CORRELATION)
  p = Ptr{Void}[0]
  cudnnCreateConvolutionDescriptor(p)
  cpadding = Cint[padding[i] for i=N:-1:1]
  cstride = Cint[stride[i] for i=N:-1:1]
  cupscale = fill(Cint(1), N)
  cudnnSetConvolutionNdDescriptor(p[1], N, cpadding, cstride, cupscale, mode, datatype(T))
  desc = ConvolutionDesc(p[1], padding, stride)
  finalizer(desc, cudnnDestroyConvolutionDescriptor)
  desc
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::ConvolutionDesc) = desc.ptr
=#

function convolution_desc(padding, stride)
  N = length(padding)
  p = Ptr{Void}[0]
  cudnnCreateConvolutionDescriptor(p)
  cpadding = Cint[padding[i] for i=N:-1:1]
  cstride = Cint[stride[i] for i=N:-1:1]
  cupscale = fill(Cint(1), N)
  cudnnSetConvolutionNdDescriptor(p[1], N, cpadding, cstride, cupscale, mode, datatype(T))
  p[1]
end

function convolution!{T}(x::CuArray{T}, w::CuArray{T}, padding, stride, y::CuArray{T};
  mode=CUDNN_CROSS_CORRELATION, alpha=1.0, beta=0.0)

  h = gethandle(x.dev)
  xdesc = tensor_desc(x)
  wdesc = filter_desc(w)
  convdesc = convolution_desc(padding, stride)
  ydesc = tensor_desc(y)

  algo_p = cudnnConvolutionFwdAlgo_t[0]
  cudnnGetConvolutionForwardAlgorithm(h, xdesc, wdesc, convdesc, ydesc,
    CUDNN_CONVOLUTION_FWD_PREFER_FASTEST, 0, algo_p)
  algo = algo_p[1]

  worksize_p = Cint[0]
  cudnnGetConvolutionForwardWorkspaceSize(h, xdesc, wdesc, desc, ydesc, algo, worksize_p)
  worksize = worksize_p[1]
  workspace = CuArray(Int8, Int(worksize))

  cudnnConvolutionForward(h, T[alpha], xdesc, x, wdesc, w, convdesc,
    algo, workspace, worksize, T[beta], ydesc, y)

  cudnnDestroyTensorDescriptor(xdesc)
  cudnnDestroyFilterDescriptor(wdesc)
  cudnnDestroyConvolutionDescriptor(convdesc)
  cudnnDestroyTensorDescriptor(ydesc)
  y
end

function convolution(x::CuArray, w::CuArray, padding, stride;
  mode=CUDNN_CROSS_CORRELATION, alpha=1.0)

  N = length(padding)
  outdims = Array(Int, N)
  for i = 1:N
    outdims[i] = (size(x,i) + 2*padding[i] - size(w,i)) ÷ stride[i] + 1
  end
  y = similar(x, outdims..., size(w,N+2), size(x,N+2))
  convolution!(x, w, padding, stride, y, mode=mode, alpha=alpha)
end

function ∇convolution_bias!{T}(dy::CuArray{T}, db::CuArray{T}; alpha=1.0, beta=0.0)
  h = gethandle(x.dev)
  dydesc = tensor_desc(dy)
  dbdesc = tensor_desc(db)
  cudnnConvolutionBackwardBias(h, T[alpha], dydesc, dy, T[beta], dbdesc, db)

  cudnnDestroyTensorDescriptor(dydesc)
  cudnnDestroyTensorDescriptor(dbdesc)
end

function ∇convolution_filter!{T}(x::CuArray{T}, dy::CuArray{T}, convdesc, dw::CuArray{T};
  alpha=1.0, beta=0.0)
  h = gethandle(x.dev)
  xdesc = tensor_desc(x)
  dydesc = tensor_desc(dy)
  dwdesc = filter_desc(dw)

  algo_p = cudnnConvolutionBwdAlgo_t[0]
  cudnnFindConvolutionBackwardFilterAlgorithm(h, xdesc, dydesc, convdesc, dwdesc,
    CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST, 0, algo_p)
  algo = algo_p[1]

  worksize_p = Cint[0]
  cudnnGetConvolutionBackwardFilterWorkspaceSize(h, xdesc, dydesc, convdesc, dwdesc, algo, worksize_p)
  worksize = worksize_p[1]
  workspace = CuArray(Int8, Int(worksize))

  cudnnConvolutionBackwardFilter(h, T[alpha], xdesc, x, dydesc, dy, convdesc,
    algo, workspace, worksize, T[beta], dwdesc, dw)

end

function ∇convolution_data!{T}(w::CuArray{T}, dy::CuArray{T}, convdesc, dx::CuArray{T};
  alpha=1.0, beta=0.0)
  wdesc = FilterDesc(w)
  dydesc = TensorDesc(dy)
  dxdesc = TensorDesc(dx)

  algo_p = cudnnConvolutionBwdDataAlgo_t[0]
  cudnnGetConvolutionBackwardDataAlgorithm(h, wdesc, dydesc, convdesc, dxdesc,
    CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST, 0, algo_p)
  algo = algo_p[1]

  worksize_p = Cint[0]
  cudnnGetConvolutionBackwardDataWorkspaceSize(h, wdesc, dydesc, convdesc,
    dxdesc, algo, worksize_p)
  worksize = worksize_p[1]
  workspace = CuArray(Int8, Int(worksize))

  cudnnConvolutionBackwardData(h, T[alpha], wdesc, w, dydesc, dy, convdesc,
    algo, workspace, T[beta], dxdesc, dx)
end


text:
module ReinforcementLearningEnvironmentViZDoom
using Reexport
@reexport using ReinforcementLearning
import ReinforcementLearning:interact!, reset!, getstate
using ViZDoom
const vz = ViZDoom

function listconsts(typ)
    allconsts = names(vz, all = true)
    idx = findall(x -> typeof(getfield(vz, x)) == getfield(vz, typ), allconsts)
    allconsts[idx]
end
list_available_buttons() = listconsts(:Button)
list_screen_resolution() = listconsts(:ScreenResolution)
list_screen_format() = listconsts(:ScreenFormat)
list_mode() = listconsts(:Mode)
function list_options()
    allconsts = names(vz, all = true)
    idx = findall(x -> length(string(x)) > 2 && string(x)[1:3] == "set", allconsts)
    map(x -> Symbol(string(x)[5:end]), allconsts[idx])
end

"""
    struct ViZDoomEnvironment
        game::vz.DoomGameAllocated
        actions::Array{Array{Float64, 1}, 1}
        sleeptime::Float64
"""
struct ViZDoomEnvironment
    game::vz.DoomGameAllocated
    actions::Array{Array{Float64, 1}, 1}
    sleeptime::Float64
end
export ViZDoomEnvironment
"""
    ViZDoomEnvironment(; add_game_args = "", kw...)

Creates a new ViZDoomEnvironment. Use `list_options()` for possible `kw` and
`list_available_buttons()`, `list_screen_resolution()`, `list_screen_format()`,
`list_mode()` for possible values.
"""
function ViZDoomEnvironment(; add_game_args = "", kw...)
    defaults = (screen_format = :GRAY8, screen_resolution = :RES_160X120,
                window_visible = false, living_reward = 0,
                episode_timeout = 500)
    config = Dict(pairs(merge(defaults, kw)))
    for (k, v) in config
        if typeof(v) == Symbol
            config[k] = getfield(vz, v)
        elseif typeof(v) <: AbstractArray && typeof(v[1]) == Symbol
            config[k] = map(x -> getfield(vz, x), v)
        end
    end
    game = vz.basic_game(; config...)
    vz.add_game_args(game, add_game_args)
    if config[:window_visible]
        sleeptime = 1.0 / vz.DEFAULT_TICRATE
    else
        sleeptime = 0.
    end
    na = haskey(config, :available_buttons) ? length(config[:available_buttons]) : 3
    env = ViZDoomEnvironment(game,
                             [Float64[i == j for i in 1:na] for j in 1:na],
                             sleeptime)
    init!(env)
    env
end

function interact!(a, env::ViZDoomEnvironment)
    r = vz.make_action(env.game, env.actions[a])
    done = vz.is_episode_finished(env.game)
    if done
        state = zeros(UInt8, vz.get_screen_size(env.game))
    else
        state = vz.get_screen_buffer(env.game)
    end
    if env.sleeptime > 0 sleep(env.sleeptime) end
    return state, r, done
end
function reset!(env::ViZDoomEnvironment)
    vz.new_episode(env.game)
    vz.get_screen_buffer(env.game)
end
function getstate(env::ViZDoomEnvironment)
    vz.get_screen_buffer(env.game), vz.is_episode_finished(env.game)
end

close!(env::ViZDoomEnvironment) = vz.close(env.game)
init!(env::ViZDoomEnvironment) = vz.init(env.game)

export list_available_buttons, list_screen_resolution,
list_screen_format, list_mode, list_options, close!, init!
end # module

max_stars_repo_name: jbrea/RLEnvVizDoom

max_issues_repo_name: jbrea/RLEnvVizDoom

max_forks_repo_name: jbrea/RLEnvVizDoom


text:
using FactCheck

function pushcleanedargs!(arglist,inputstring)
    for arg in split(strip(inputstring),",")
        cleaned_arg = strip(arg)
        if cleaned_arg != ""
            cleaned_arg = split(cleaned_arg,"(")[1]
            if ismatch(r"\)",cleaned_arg)
                continue
            end
            push!(arglist,cleaned_arg);
        end
    end
end

test_clean1 = "a,b,c"
test_clean2 ="typeboun(*),inpc(*)"
test_output = []
test_output2 = []
facts("Testing pushcleanedargs!(arglist,inputstring)") do
    pushcleanedargs!(test_output,test_clean1)
    @fact test_output --> ["a","b","c"]
    pushcleanedargs!(test_output2,test_clean2)
    @fact test_output2 --> ["typeboun","inpc"]
end

function find_subroutine_args(file_lines)
    activeline = false
    subroutineargs = []
    #fil = open(file_name)

    for line in file_lines
        if activeline & ismatch(r"\)",line)
            activeline = false
            line_middle_part = split(split(line,"&")[2],")")[1]
            pushcleanedargs!(subroutineargs,line_middle_part)
        end
        if activeline
            line_end_part = split(line,"&")[2]
            pushcleanedargs!(subroutineargs,line_end_part)
        end
        if ismatch(r"subroutine",line)
            activeline = true
            line_end_part = split(line,"(")[2]
            pushcleanedargs!(subroutineargs,line_end_part)
        end
    end
    #close(fil)
    return subroutineargs
end

test_data1 = ["      subroutine calinput(co,nk,"
              "     &  ne,nodeboun)"
              "! comment line"]

facts("Testing find_subroutine_args(file_lines)") do
    @fact find_subroutine_args(test_data1) --> ["co","nk","ne","nodeboun"]
end

function findalldefinitions(file_lines)
    activeline = false
    continuecomp = "empty"
    #fil = open(file_name)
    output = Dict()
    keyw = ["logical" "character" "integer" "real"]
    #for key in keyw
    #    output[key] = []
    #end
    for line in file_lines
        if activeline
            if ismatch(r"\&",line)
                #println(line)
                m = match(r"(\s*&\s+)(.*)",line)
                pushcleanedargs!(output[continuecomp],m.captures[2])
            else
                activeline = false
            end
        end
        for comp in keyw
            if ismatch(Regex(comp),line)
                activeline = true
                m = match(r"(\s*\w+\s+)(.*)",line)
                if ~(comp in keys(output))
                    output[comp] = []
                end
                pushcleanedargs!(output[comp],m.captures[2])
                continuecomp = comp
            end
        end
    end
    #close(fil)
    return output
end

test_def1 = ["      logical boun_flag,cload_flag,"
             "     &  nodeprint_flag"]

test_def2 = ["      character*1 typeboun(*),inpc(*)"
             "      character*3 output"]

test_def3 = ["      integer kon(*),nodempc(3,*),"
             "     &  nodeforc(2,*),iaxial,"
             "     &  istartset(*)"]

facts("Testing findalldefinitions(file_lines)") do
    @fact findalldefinitions(test_def1) -->
        Dict{Any,Any}("logical" => Any["boun_flag","cload_flag","nodeprint_flag"])
    @fact findalldefinitions(test_def2) -->
        Dict{Any,Any}("character" => Any["typeboun","inpc","output"])
    @fact findalldefinitions(test_def3) -->
        Dict{Any,Any}("integer" => Any["kon","nodempc","nodeforc","iaxial","istartset"])
end

function ptrtypes(arg,dictoftypes)
    for key in keys(dictoftypes)
        if arg in dictoftypes[key]
            if key == "real"
                return "Ptr{Float64}"
            elseif key == "integer"
                return "Ptr{Int64}"
            elseif key == "logical"
                return "Ptr{Bool}"
            elseif key == "character"
                return "Ptr{ASCIIString}"
            end
        end
    end
end

function writeccall(args,dictoftypes)
    out = @sprintf "ccall((:%s_, ../deps/usr/lib/libcalculix.so), Int64,\n" "calinput"
    final = length(args)
    for (num,arg) in enumerate(args)
        if num == final
            out = out * @sprintf "%s),\n       " ptrtypes(arg,dictoftypes)
            break
        end
        if num == 1
            out = out * "       ("
        elseif (num % 5) == 0
            out = out * "\n        "
        else
            out = out * @sprintf "%s," ptrtypes(arg,dictoftypes)
        end
    end
    for (num,arg) in enumerate(args)
        if num == final
            out = out * @sprintf "%s)\n" arg
            break
        end
        if (num % 4) == 0
            out = out * "\n       "
        end
        out = out * @sprintf "%s," arg
    end
    println(out)
end
fil = open("../deps/src/CalculiX-cmake/src/fortran/calinput.f")
lines = readlines(fil)
close(fil)
arglist = find_subroutine_args(lines)
argdefs = findalldefinitions(lines)
writeccall(arglist,argdefs)


text:
@testset "Deque" begin

    # empty dequeue

    q = Deque{Int}()
    @test length(q) == 0
    @test isempty(q)
    @test q.blksize == DataStructures.DEFAULT_DEQUEUE_BLOCKSIZE
    @test_throws ArgumentError front(q)
    @test_throws ArgumentError back(q)
    @test length(sprint(dump,q)) >= 0

    @test typeof(deque(Int)) === typeof(Deque{Int}())

    q = DataStructures.DequeBlock{Int}(0,0)
    @test length(q) == 0
    @test capacity(q) == 0
    @test isempty(q)
    @test length(sprint(show,q)) >= 0

    q = Deque{Int}(3)
    @test length(q) == 0
    @test isempty(q)
    @test q.blksize == 3
    @test num_blocks(q) == 1
    @test_throws ArgumentError front(q)
    @test_throws ArgumentError back(q)
    @test isa(collect(q), Vector{Int})
    @test collect(q) == Int[]

    # push back

    n = 10

    for i = 1 : n
        push!(q, i)
        @test length(q) == i
        @test isempty(q) == false
        @test num_blocks(q) == div(i-1, 3) + 1

        @test front(q) == 1
        @test back(q) == i

        k = 1
        for j in q
            @test j == k
            k += 1
        end
        cq = collect(q)
        @test isa(cq, Vector{Int})
        @test cq == collect(1:i)
        @test length(sprint(show,q)) >= 0
    end

    # pop back

    for i = 1 : n
        x = pop!(q)
        @test length(q) == n - i
        @test isempty(q) == (i == n)
        @test num_blocks(q) == div(n-i-1, 3) + 1
        @test x == n - i + 1

        if !isempty(q)
            @test front(q) == 1
            @test back(q) == n - i
        else
            @test_throws ArgumentError front(q)
            @test_throws ArgumentError back(q)
        end

        cq = collect(q)
        @test cq == collect(1:n-i)
    end

    # push front

    q = Deque{Int}(3)

    for i = 1 : n
        unshift!(q, i)
        @test length(q) == i
        @test isempty(q) == false
        @test num_blocks(q) == div(i-1, 3) + 1

        @test front(q) == i
        @test back(q) == 1

        cq = collect(q)
        @test isa(cq, Vector{Int})
        @test cq == collect(i:-1:1)
    end

    # pop front

    for i = 1 : n
        x = shift!(q)
        @test length(q) == n - i
        @test isempty(q) == (i == n)
        @test num_blocks(q) == div(n-i-1, 3) + 1
        @test x == n - i + 1

        if !isempty(q)
            @test front(q) == n - i
            @test back(q) == 1
        else
            @test_throws ArgumentError front(q)
            @test_throws ArgumentError back(q)
        end

        cq = collect(q)
        @test cq == collect(n-i:-1:1)
    end

    # random operations

    q = Deque{Int}(5)
    r = Int[]
    m = 100

    for k = 1 : m
        la = rand(1:20)
        x = rand(1:1000, la)

        for i = 1 : la
            if rand(Bool)
                push!(r, x[i])
                push!(q, x[i])
            else
                unshift!(r, x[i])
                unshift!(q, x[i])
            end
        end

        @test length(q) == length(r)
        @test collect(q) == r

        lr = rand(1:length(r))
        for i = 1 : lr
            if rand(Bool)
                pop!(r)
                pop!(q)
            else
                shift!(r)
                shift!(q)
            end
        end

        @test length(q) == length(r)
        @test collect(q) == r
    end

    # hash and ==
    a = Deque{Int64}(2)
    b = Deque{Int64}(3)
    # Note: blksize does not distinguish == or hash
    @test a == b
    @test hash(a) === hash(b)
    push!(a, 1)
    @test a != b
    push!(a, 2)
    push!(b, 2, 1)
    @test a != b
    @test hash(a) !== hash(b)
    shift!(b)
    push!(b, 2)
    @test a == b
    @test hash(a) == hash(b)

    # issue #38

    q = Deque{Int}(1)
    push!(q,1)
    @test !isempty(q)
    empty!(q)
    @test isempty(q)

    #empty!
    q = Deque{Int}(1)
    push!(q,1)
    push!(q,2)
    @test length(sprint(dump,q)) >= 0
    @test typeof(empty!(q)) === typeof(Deque{Int}())
    @test isempty(q)

end # @testset Deque


text:
export lse, lomax, randomindex

# Log-sum exp.
function lse(x::AbstractMatrix{<:Real}; dims = 2)
    m = maximum(x, dims = dims)
    v = exp.(x .- m)
    return log.(sum(v, dims = dims)) + m
end

# see https://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf
function tpdf(μ::T, σ::T, ν::T, x::V) where {T<:Real,V<:Real}
    c = gamma(ν/2 + 1/2) / gamma(ν/2)
    c *= 1/(sqrt(ν*π)*σ)
    return c * (1 + 1/ν * ((x-μ)/σ)^2)^((ν+1)/2)
end

function tlogpdf(μ::T, σ::T, ν::T, x::V) where {T<:Real,V<:Real}
    c = lgamma(ν/2 + 1/2) - lgamma(ν/2)
    c -= log(sqrt(ν*π)*σ)
    return c + ((ν+1)/2)*log(1 + 1/ν * ((x-μ)/σ)^2)
end

@inline lomax(x, α, λ) = log(α) - log(λ) + -(α + 1) * log(1 + (x / λ))

function randomindex(p::AbstractVector{T}) where T<: Real
    return randomindex(p, rand(T))
end

function randomindex(p::AbstractVector{T}, rnd::T) where T<: Real

    csum = T(0.0)
    thres = rnd * sum(p)

    for i in 1:length(p)
        @inbounds begin
            csum += p[i]
            if csum >= thres
                return Int(i)
            end
        end
    end

    return Int(length(p))
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

p = [1.0 2.0 3.0]

bp = CutCellDG.extend_to_face(p, 1)
testbp = vcat(p, -ones(3)')
@test allapprox(bp, testbp)

rp = CutCellDG.extend_to_face(p, 2)
testrp = vcat(ones(3)', p)
@test allapprox(rp, testrp)

tp = CutCellDG.extend_to_face(p, 3)
testtp = vcat(p, ones(3)')
@test allapprox(tp, testtp)

lp = CutCellDG.extend_to_face(p, 4)
testlp = vcat(-ones(3)', p)
@test allapprox(testlp, lp)

quad1d = tensor_product_quadrature(1, 4)
points = quad1d.points
bp = vcat(points, -ones(4)')
rp = vcat(ones(4)', points)
tp = vcat(points, ones(4)')
lp = vcat(-ones(4)', points)
bq = CutCellDG.extend_to_face(quad1d, 1)
rq = CutCellDG.extend_to_face(quad1d, 2)
tq = CutCellDG.extend_to_face(quad1d, 3)
lq = CutCellDG.extend_to_face(quad1d, 4)
@test allapprox(bq.points, bp)
@test allapprox(rq.points, rp)
@test allapprox(tq.points, tp)
@test allapprox(lq.points, lp)
@test allapprox(bq.weights, quad1d.weights)
@test allapprox(rq.weights, quad1d.weights)
@test allapprox(tq.weights, quad1d.weights)
@test allapprox(lq.weights, quad1d.weights)

facequads = CutCellDG.face_quadratures(4)
@test allapprox(facequads[1].points, bp)
@test allapprox(facequads[2].points, rp)
@test allapprox(facequads[3].points, tp)
@test allapprox(facequads[4].points, lp)
@test allapprox(facequads[1].weights, quad1d.weights)
@test allapprox(facequads[2].weights, quad1d.weights)
@test allapprox(facequads[3].weights, quad1d.weights)
@test allapprox(facequads[4].weights, quad1d.weights)

polyorder = 1
numqp = 2
levelsetbasis = HermiteTensorProductBasis(2)

elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
dim, levelsetnf = size(interpolation_points(levelsetbasis))
points = interpolation_points(elasticitybasis)
levelset = InterpolatingPolynomial(1, levelsetbasis)
x0 = [0.0, 0.0]
meshwidths = [2.0, 1.0]
nelements = [2, 1]
dgmesh = CutCellDG.DGMesh(x0, meshwidths, nelements, points)
cgmesh = CutCellDG.CGMesh(x0, meshwidths, nelements, levelsetnf)

normal = [1.0, 0.0]
xI = [0.5, 0.0]

levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    cgmesh,
    levelsetbasis,
)

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)
idx = [5, 9, 6, 10, 7, 11, 8, 12, 1, 0, 2, 0, 3, 0, 4, 0]
testfacetoquad = reshape(idx, 2, 4, 2)
@test allequal(testfacetoquad, facequads.facetoquad)


text:
function testNLEQ1()
    refSol = Dict{Int32,Any}();

    refSol[2] = [0.21132486540517281, 0.78867513459482719];

    refSol[3] = [0.14644660863174364, 0.5, 0.85355339136825636];

    refSol[4] = [0.10267276384825177, 0.40620376278203324, 0.59379623721796682,
                    0.89732723615174825];

    refSol[5] = [0.083751256499240576, 0.31272929522379189, 0.5, 0.68727070477620811,
                    0.91624874350075947];

    refSol[6] = [0.066876587754941208, 0.28874077956671557, 0.36668218955839393,
                    0.63331781044160607, 0.71125922043328438, 0.93312341224505879];

    refSol[7] = [0.058069149620871681, 0.23517161236154796, 0.33804409473596703, 0.5,
                    0.66195590526403303, 0.76482838763845207, 0.94193085037912838];

    refSol[8] = [0.0475500084536796, 0.238025867229826, 0.243145948398381,
                    0.502027800864952, 0.497972199135048, 0.756854051601463,
                    0.76197413277033,0.95244999154632];

    refSol[9] = [0.044205346135779103, 0.19949067230774983, 0.23561910847320636,
                    0.41604690789257154, 0.50000000000000011, 0.58395309210742818,
                    0.76438089152679367, 0.80050932769225003, 0.95579465386422091];

    testResult = false;
    for dim = 2:9
        n1 = dim + 1

        # Initialize the options
        opt = OptionsNLEQ(OPT_JACGEN            => 1,
                          OPT_JACFCN            => chebyQuadJac,
                          OPT_MSTOR             => 0,
                          OPT_NOROWSCAL         => 0,
                          OPT_NITMAX            => 10,
                          OPT_RTOL              => 1e-5)

        x0    = collect(1:dim)./n1
        xScal = zeros(x0)

        (x0, _, retCode) = nleq1(chebyQuad,x0,xScal,opt)

        relNormDiff = norm(x0-refSol[dim],Inf)/norm(refSol[dim],Inf)

        # Here we check < 1e-12 due to the fact that we use the Julia LU Decomposition
        # rather than the same code for factorization and solution as with NLEQ1
        if dim == 2
            testResult = relNormDiff < 1e-12
        else
            testResult &= relNormDiff < 1e-12
        end
    end
    return testResult
end

testNLEQ1()


text:
# The [Neherlab COVID-19](https://neherlab.org/covid19/) forecast model
include("COVID-19-utils.jl")
include("COVID-19-model.jl")
include("COVID-19-run-model.jl")
include("COVID-19-data.jl")

using BlackBoxOptim

# If running for the first time, or no updates for a long time
# updateData()

COUNTRY_LIST = [("France",                      :north),
                ("Germany",                     :north),
                ("Italy",                       :north),
                ("Spain",                       :north),
                ("Switzerland",                 :north),
                ("United Kingdom",              :north),
                ("United States of America",    :north)]

countryData = Dict( c => populateCountryDate(c, h) for (c, h) in COUNTRY_LIST)

#--------------------------------------------------------------------------------------------------
#--
#-- Find optimal parameters for every country
#--
for (c, _) in COUNTRY_LIST
    global country = c

    # Determine optimal parameters - 15 seconds per country
    result = bboptimize(deathsLoss, SearchRange = countryData[country][:range], MaxTime = 15)
    global countryData[country][:params] = best_candidate(result)
end

# Select a country
country = "Spain"

#--------------------------------------------------------------------------------------------------
#--
#-- Default parameters
#--
sol = calculateSolution(country, countryData[country][:params])
calculateTotalDeaths(sol)

#--------------------------------------------------------------------------------------------------
#-- Plot efault parameters
#--
pyplot();
clf();
ioff();

fig, ax = PyPlot.subplots();

ax.plot(countryData[country][:cases].time,
        calculateTotalDeaths(calculateSolution(country, countryData[country][:params])),
        label = "Forecast");

ax.plot(countryData[country][:cases].time, countryData[country][:cases].deaths, "ro", label = "Actual", alpha = 0.3);

ax.legend(loc="lower right");
ax.set_xlabel("time");
ax.set_ylabel("Individuals");
ax.set_yscale("log");

gcf()

#--------------------------------------------------------------------------------------------------
#-- Optimal parameters
#--
result = bboptimize(deathsLoss, SearchRange = countryData[country][:range], MaxTime = 15)
countryData[country][:params] = best_candidate(result)

#-------------------------------------------------------------------------------------------------
#--
#-- Global optimisation across all countries
#--

# One run is optimising the disease, then optimising the countries.
for runs in 1:5

    # First freeze the countries and optimise the epidemiology
    freezeCountryRange()

    # Optimise across all countries
    result = bboptimize(allCountriesLoss, SearchRange = countryData["France"][:range], MaxTime = 30)

    #Update the parameters of each country
    for (c, _) in COUNTRY_LIST
        global countryData[c][:params] = (diseaseMask .* best_candidate(result)) .+ (countryMask .* countryData[c][:params])
    end

    #-- Freeze the parameters related to the disease
    freezeDiseaseRange()

    # Then optimise each country
    for (c, _) in COUNTRY_LIST
        # Make a note of the disease parameters
        p = countryData[c][:params]

        # Determine optimal parameters - 15 seconds per country
        result = bboptimize(deathsLoss, SearchRange = countryData[country][:range], MaxTime = 15)
        global countryData[c][:params] = (diseaseMask .* p) .+ (countryMask .* best_candidate(result))
    end
end

allParams = [countryData[c][:params] for (c, _) in COUNTRY_LIST]

using Printf

# Do what I say, not what I do
Base.show(io::IO, f::Float64) = @printf io "%1.3f" f

@show allParams


text:
#using MatrixNetworks

function load_netalign_problem(probname)

if probname == "lcsh2wiki-small"
    mainloc = joinpath(dirname(pathof(NetworkAlign)),"../data/lcsh2wiki-small")
    # mainloc = "../data/lcsh2wiki-small"
    location = join([mainloc,"_A.smat"])
    A = MatrixNetworks.readSMAT(location)

    location = join([mainloc,"_B.smat"])
    B = MatrixNetworks.readSMAT(location)

    location = join([mainloc,"_L.smat"])
    (rows,header) = readdlm(location;header=true)
    L = sparse(
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1,
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),2]).+1,
               rows[1:parse(Int,header[3]),3],
               parse(Int,header[1]),
               parse(Int,header[2])
               )

    location = join([mainloc,"_S.smat"])
    S = MatrixNetworks.readSMAT(location)

    location = join([mainloc,"_li.smat"])
    (rows,header) = readdlm(location;header=true)
    li = zeros(Int64,parse(Int,header[1]))
    ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1
    li[ri] = convert(Array{Int64,1},rows[1:parse(Int,header[3]),3])

    location = join([mainloc,"_lj.smat"])
    (rows,header) = readdlm(location;header=true)
    lj = zeros(Int64,parse(Int,header[1]))
    ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1
    lj[ri] = convert(Array{Int64,1},rows[1:parse(Int,header[3]),3])

    location = join([mainloc,"_lw.smat"])
    (rows,header) = readdlm(location;header=true)
    w = zeros(Float64,parse(Int,header[1]))
    ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1
    w[ri] = rows[1:parse(Int,header[3]),3]

    return (S,w,li,lj,A,B,L)

elseif probname == "example-overlap"

    mainloc = joinpath(dirname(pathof(NetworkAlign)),"../data/example-overlap")
    # mainloc = "../data/example-overlap"
    location = join([mainloc,"_A.smat"])
    A = MatrixNetworks.readSMAT(location)

    location = join([mainloc,"_B.smat"])
    B = MatrixNetworks.readSMAT(location)

    location = join([mainloc,"_L.smat"])
    (rows,header) = readdlm(location;header=true)
    #println(rows[1:parse(Int,header[3]),1])
    L = sparse(
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1,
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),2]).+1,
               rows[1:parse(Int,header[3]),3],
               parse(Int,header[1]),
               parse(Int,header[2])
               )

    location = join([mainloc,"_S.smat"])
    S = MatrixNetworks.readSMAT(location)

    location = join([mainloc,"_li.smat"])
    (rows,header) = readdlm(location;header=true)
    li = zeros(Int64,parse(Int,header[1]))
    ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1
    li[ri] = convert(Array{Int64,1},rows[1:parse(Int,header[3]),3])

    location = join([mainloc,"_lj.smat"])
    (rows,header) = readdlm(location;header=true)
    lj = zeros(Int64,parse(Int,header[1]))
    ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1
    lj[ri] = convert(Array{Int64,1},rows[1:parse(Int,header[3]),3])

    location = join([mainloc,"_lw.smat"])
    (rows,header) = readdlm(location;header=true)
    w = zeros(Float64,parse(Int,header[1]))
    ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1]).+1
    w[ri] = rows[1:parse(Int,header[3]),3]

    return (S,w,li,lj,A,B,L)
elseif probname == "lcsh2wiki-full"
  mainloc = joinpath(dirname(pathof(NetworkAlign)),"data",probname)

    location = join([mainloc,"_S.smat"])
    S = MatrixNetworks.readSMAT(location)

    # used vmat storage because the smat files were too big for github
    location = join([mainloc,"_li.vmat"])
    rows = readdlm(location)
    li = convert(Array{Int64,1},rows[:])

    location = join([mainloc,"_lj.vmat"])
    rows = readdlm(location)
    lj = convert(Array{Int64,1},rows[:])

    location = join([mainloc,"_lw.vmat"])
    rows = readdlm(location)
    w = rows[:]

    # location = join([mainloc,"_li.smat"])
    # (rows,header) = readdlm(location;header=true)
    # li = zeros(Int64,parse(Int,header[1]))
    # ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1])+1
    # li[ri] = convert(Array{Int64,1},rows[1:parse(Int,header[3]),3])
    #
    # location = join([mainloc,"_lj.smat"])
    # (rows,header) = readdlm(location;header=true)
    # lj = zeros(Int64,parse(Int,header[1]))
    # ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1])+1
    # lj[ri] = convert(Array{Int64,1},rows[1:parse(Int,header[3]),3])
    #
    # location = join([mainloc,"_lw.smat"])
    # (rows,header) = readdlm(location;header=true)
    # w = zeros(Float64,parse(Int,header[1]))
    # ri = convert(Array{Int64,1},rows[1:parse(Int,header[3]),1])+1
    # w[ri] = rows[1:parse(Int,header[3]),3]

    return (S,w,li,lj,[],[],[])

else
    error("NetAlign: problem name not recognized")
end
end


text:
@testset "serialize" begin

    @testset "MVHistory" begin
        expected = MVHistory()
        push!(expected, :a, 1, 1.2)
        push!(expected, :a, 2, 2.2)
        push!(expected, :b, 1, 3)
        push!(expected, :b, 2, 4)
        push!(expected, :c, 1, OneClassActiveLearning.ConfusionMatrix(1,1,1,1))
        push!(expected, :c, 2, OneClassActiveLearning.ConfusionMatrix(2,2,2,2))

        actual = Unmarshal.unmarshal(MVHistory, JSON.parse(JSON.json(expected)))

        for k in keys(expected)
            @test expected[k] == actual[k]
        end
    end

    @testset "al summary" begin
        expected = Dict(:m1 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]),
                           :m2 => Dict(:a => 1.0, :b => 1, :c => [1.0], :d => [1, 2]))

        actual = OneClassActiveLearning.unmarshal_al_summary(JSON.parse(JSON.json(expected)))

        for (k_metrics, v_metrics) in expected
            for (k_summary, v_summary) in v_metrics
                @test expected[k_metrics][k_summary] == actual[k_metrics][k_summary]
            end
        end
    end

    @testset "Result" begin
        id = 42
        experiment = Dict{Symbol, Any}(
                :hash => id,
                :data_file => "$(@__DIR__)/../data/input/dummy.csv",
                :output_file => "OneClassActiveLearning.jl/data/output/scenarioA/data_qs_model_id.tmp",
                :model => Dict( :type => :(SVDD.RandomOCClassifier),
                                :init_strategy => :(SVDD.FixedParameterInitialization(GaussianKernel(2), 0.5))),
                :query_strategy => Dict(:type => :(OneClassActiveLearning.QueryStrategies.RandomQs), :param => Dict()),
                :split_strategy => OneClassActiveLearning.DataSplits(trues(123), OneClassActiveLearning.FullSplitStrat()),
                :param => Dict(:num_al_iterations => 5,
                               :solver => :(IpoptSolver),
                               :initial_pools => fill(:U, 123)))

       al_history = MVHistory()
       push!(al_history, :a, 1, 1.2)
       push!(al_history, :a, 2, 2.2)
       push!(al_history, :b, 1, 3)
       push!(al_history, :b, 2, 4)
       push!(al_history, :query_label, 1, true)
       push!(al_history, :query_label, 2, false)
       push!(al_history, :c, 1, OneClassActiveLearning.ConfusionMatrix(1,1,1,1))
       push!(al_history, :c, 2, OneClassActiveLearning.ConfusionMatrix(2,2,2,2))

        al_summary = Dict(:m1 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]),
                   :m2 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]))

        worker_info = OneClassActiveLearning.get_worker_info()
        expected = OneClassActiveLearning.Result(id, experiment, worker_info, DataStats(1, 1, 0.1, 0.1, [1,2,3,4], [5,6,7]), al_history, al_summary)

        actual = Unmarshal.unmarshal(OneClassActiveLearning.Result, JSON.parse(JSON.json(expected)))
        @test actual.id == 42
        @test length(actual.experiment[:param][:initial_pools]) == 123
        @test haskey(actual.worker_info, :hostname)
        @test haskey(actual.al_history, :c)
        @test haskey(actual.al_summary, :m1)
        @test haskey(actual.status, :exit_code)
    end
end


text:
include("abstract_explorer.jl")
include("batch_explorer.jl")
include("epsilon_greedy_explorer.jl")
include("UCB_explorer.jl")
include("weighted_explorer.jl")
include("gumbel_softmax_explorer.jl")
# TODO
# https://arxiv.org/abs/1509.02971
# Boltzmann exploration
# additive Gaussian noise
# additive Ornstein-Uhlenbeck noise


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""

use Scenario struct to create reopt.jl model inputs
"""

struct REoptInputs <: AbstractInputs
    techs::Array{String, 1}
    pvtechs::Array{String, 1}
    gentechs::Array{String,1}
    elec_techs::Array{String, 1}
    segmented_techs::Array{String, 1}
    pbi_techs::Array{String, 1}
    techs_no_turndown::Array{String, 1}
    min_sizes::DenseAxisArray{Float64, 1}  # (techs)
    max_sizes::DenseAxisArray{Float64, 1}  # (techs)
    existing_sizes::DenseAxisArray{Float64, 1}  # (techs)
    cap_cost_slope::Dict{String, Any}  # (techs)
    om_cost_per_kw::DenseAxisArray{Float64, 1}  # (techs)
    elec_load::ElectricLoad
    time_steps::UnitRange
    time_steps_with_grid::Array{Int, 1}
    time_steps_without_grid::Array{Int, 1}
    hours_per_timestep::Float64
    months::UnitRange
    production_factor::DenseAxisArray{Float64, 2}  # (techs, time_steps)
    levelization_factor::Dict{String, Float64}  # (techs)
    VoLL::Array{R, 1} where R<:Real #default set to 1 US dollar per kwh
    pwf_e::Float64
    pwf_om::Float64
    two_party_factor::Float64
    owner_tax_pct::Float64
    offtaker_tax_pct::Float64
    microgrid_premium_pct::Float64
    pvlocations::Array{Symbol, 1}
    maxsize_pv_locations::DenseAxisArray{Float64, 1}  # indexed on pvlocations
    pv_to_location::DenseAxisArray{Int, 2}  # (pvtechs, pvlocations)
    etariff::ElectricTariff
    ratchets::UnitRange
    techs_by_exportbin::Dict{Symbol, AbstractArray}  # keys can include [:NEM, :WHL, :CUR]
    storage::Storage
    generator::Generator
    elecutil::ElectricUtility
    min_resil_timesteps::Int
    mg_tech_sizes_equal_grid_sizes::Bool
    node::Int
    export_bins_by_tech::Dict
    n_segs_by_tech::Dict{String, Int}
    seg_min_size::Dict{String, Dict{Int, Float64}}
    seg_max_size::Dict{String, Dict{Int, Float64}}
    seg_yint::Dict{String, Dict{Int, Float64}}
    pbi_pwf::Dict{String, Any}  # (pbi_techs)
    pbi_max_benefit::Dict{String, Any}  # (pbi_techs)
    pbi_max_kw::Dict{String, Any}  # (pbi_techs)
    pbi_benefit_per_kwh::Dict{String, Any}  # (pbi_techs)
end

"""
    REoptInputs(fp::String)

Use `fp` to load in JSON scenario:
```
function REoptInputs(fp::String)
    s = Scenario(JSON.parsefile(fp))
    REoptInputs(s)
end
```
Useful if you want to manually modify REoptInputs before solving the model.
"""
function REoptInputs(fp::String)
    s = Scenario(JSON.parsefile(fp))
    REoptInputs(s)
end

"""
    REoptInputs(s::Scenario)

Constructor for REoptInputs
"""
function REoptInputs(s::Scenario)

    time_steps = 1:length(s.electric_load.loads_kw)
    hours_per_timestep = 1 / s.settings.time_steps_per_hour
    techs, pvtechs, gentechs, elec_techs, segmented_techs, pv_to_location, maxsize_pv_locations, pvlocations,
        production_factor, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, n_segs_by_tech,
        seg_min_size, seg_max_size, seg_yint, techs_by_exportbin, export_bins_by_tech = setup_tech_inputs(s)
    techs_no_turndown = copy(pvtechs)
    if "Wind" in techs
        append!(techs_no_turndown, ["Wind"])
    end

    pbi_techs, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh = setup_pbi_inputs(s, techs)

    months = 1:length(s.electric_tariff.monthly_demand_rates)

    levelization_factor, pwf_e, pwf_om, two_party_factor = setup_present_worth_factors(s, techs, pvtechs)
    # the following hardcoded value for levelization_factor matches the public REopt API value
    # for test_with_cplex (test_time_of_export_rate) and makes the test values match.
    # the REopt code herein uses the Desktop method for levelization_factor, which is more accurate
    # (Desktop has non-linear degradation vs. linear degradation in API)
    # levelization_factor = Dict("PV" => 0.9539)
    # levelization_factor = Dict("PV" => 0.9539, "Generator" => 1.0)  # w/generator
    time_steps_with_grid, time_steps_without_grid, = setup_electric_utility_inputs(s)

    if any(pv.existing_kw > 0 for pv in s.pvs)
        adjust_load_profile(s, production_factor)
    end

    REoptInputs(
        techs,
        pvtechs,
        gentechs,
        elec_techs,
        segmented_techs,
        pbi_techs,
        techs_no_turndown,
        min_sizes,
        max_sizes,
        existing_sizes,
        cap_cost_slope,
        om_cost_per_kw,
        s.electric_load,
        time_steps,
        time_steps_with_grid,
        time_steps_without_grid,
        hours_per_timestep,
        months,
        production_factor,
        levelization_factor,
        typeof(s.financial.VoLL) <: Array{<:Real, 1} ? s.financial.VoLL : fill(s.financial.VoLL, length(time_steps)),
        pwf_e,
        pwf_om,
        two_party_factor,
        s.financial.owner_tax_pct,
        s.financial.offtaker_tax_pct,
        s.financial.microgrid_premium_pct,
        pvlocations,
        maxsize_pv_locations,
        pv_to_location,
        s.electric_tariff,
        1:length(s.electric_tariff.tou_demand_ratchet_timesteps),  # ratchets
        techs_by_exportbin,
        s.storage,
        s.generator,
        s.electric_utility,
        s.site.min_resil_timesteps,
        s.site.mg_tech_sizes_equal_grid_sizes,
        s.site.node,
        export_bins_by_tech,
        n_segs_by_tech,
        seg_min_size,
        seg_max_size,
        seg_yint,
        pbi_pwf,
        pbi_max_benefit,
        pbi_max_kw,
        pbi_benefit_per_kwh
    )
end

"""
    function setup_tech_inputs(s::Scenario)

Create data arrays associated with techs necessary to build the JuMP model.
"""
function setup_tech_inputs(s::Scenario)

    pvtechs = String[pv.name for pv in s.pvs]
    if length(Base.Set(pvtechs)) != length(pvtechs)
        error("PV names must be unique, got $(pvtechs)")
    end

    techs = copy(pvtechs)
    gentechs = String[]
    segmented_techs = String[]
    if s.wind.max_kw > 0
        push!(techs, "Wind")
    end
    if s.generator.max_kw > 0
        push!(techs, "Generator")
        push!(gentechs, "Generator")
    end

    elec_techs = copy(techs)  # only modeling electric loads/techs so far

    time_steps = 1:length(s.electric_load.loads_kw)

    # REoptInputs indexed on techs:
    max_sizes = DenseAxisArray{Float64}(undef, techs)
    min_sizes = DenseAxisArray{Float64}(undef, techs)
    existing_sizes = DenseAxisArray{Float64}(undef, techs)
    cap_cost_slope = Dict{String, Any}()
    om_cost_per_kw = DenseAxisArray{Float64}(undef, techs)
    production_factor = DenseAxisArray{Float64}(undef, techs, time_steps)

    # export related inputs
    techs_by_exportbin = Dict(k => [] for k in s.electric_tariff.export_bins)
    export_bins_by_tech = Dict{String, Array{Symbol, 1}}()

    #REoptInputs indexed on segmented_techs
    n_segs_by_tech = Dict{String, Int}()
    seg_min_size = Dict{String, Any}()
    seg_max_size = Dict{String, Any}()
    seg_yint = Dict{String, Any}()

    # PV specific arrays
    pvlocations = [:roof, :ground, :both]
    pv_to_location = DenseAxisArray{Int}(undef, pvtechs, pvlocations)
    maxsize_pv_locations = DenseAxisArray([1.0e5, 1.0e5, 1.0e5], pvlocations)
    # default to large max size per location. Max size by roof, ground, both

    if !isempty(pvtechs)
        setup_pv_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor,
                        pvlocations, pv_to_location, maxsize_pv_locations, segmented_techs, n_segs_by_tech,
                        seg_min_size, seg_max_size, seg_yint, techs_by_exportbin)
    end

    if "Wind" in techs
        setup_wind_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor,
            techs_by_exportbin)
    end

    if "Generator" in techs
        setup_gen_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor,
            techs_by_exportbin)
    end

    # filling export_bins_by_tech MUST be done after techs_by_exportbin has been filled in
    for t in elec_techs
        export_bins_by_tech[t] = [bin for (bin, ts) in techs_by_exportbin if t in ts]
    end

    return techs, pvtechs, gentechs, elec_techs, segmented_techs, pv_to_location, maxsize_pv_locations, pvlocations,
    production_factor, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, n_segs_by_tech,
    seg_min_size, seg_max_size, seg_yint, techs_by_exportbin, export_bins_by_tech
end

"""
    setup_pbi_inputs(s::Scenario, techs::Array{String, 1})

Create data arrays for production based incentives.
All arrays can be empty if no techs have production_incentive_per_kwh > 0.
"""
function setup_pbi_inputs(s::Scenario, techs::Array{String, 1})

    pbi_techs = String[]
    pbi_pwf = Dict{String, Any}()
    pbi_max_benefit = Dict{String, Any}()
    pbi_max_kw = Dict{String, Any}()
    pbi_benefit_per_kwh = Dict{String, Any}()

    for tech in techs
        T = typeof(eval(Meta.parse(tech)))
        if :production_incentive_per_kwh in fieldnames(T)
            if eval(Meta.parse("s."*tech*".production_incentive_per_kwh")) > 0
                push!(pbi_techs, tech)
                pbi_pwf[tech], pbi_max_benefit[tech], pbi_max_kw[tech], pbi_benefit_per_kwh[tech] =
                    production_incentives(eval(Meta.parse("s."*tech)), s.financial)
            end
        end
    end
    return pbi_techs, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh
end

function setup_pv_inputs(s::Scenario, max_sizes, min_sizes,
    existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor,
    pvlocations, pv_to_location, maxsize_pv_locations,
    segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs_by_exportbin)

    pv_roof_limited, pv_ground_limited, pv_space_limited = false, false, false
    roof_existing_pv_kw, ground_existing_pv_kw, both_existing_pv_kw = 0.0, 0.0, 0.0
    roof_max_kw, land_max_kw = 1.0e5, 1.0e5

    for pv in s.pvs
        production_factor[pv.name, :] = prodfactor(pv, s.site.latitude, s.site.longitude)
        for location in pvlocations
            if pv.location == location
                pv_to_location[pv.name, location] = 1
            else
                pv_to_location[pv.name, location] = 0
            end
        end

        beyond_existing_kw = pv.max_kw
        if pv.location == "both"
            both_existing_pv_kw += pv.existing_kw
            if !(s.site.roof_squarefeet === nothing) && !(s.site.land_acres === nothing)
                # don"t restrict unless both land_area and roof_area specified,
                # otherwise one of them is "unlimited"
                roof_max_kw = s.site.roof_squarefeet * pv.kw_per_square_foot
                land_max_kw = s.site.land_acres / pv.acres_per_kw
                beyond_existing_kw = min(roof_max_kw + land_max_kw, beyond_existing_kw)
                pv_space_limited = true
            end
        elseif pv.location == "roof"
            roof_existing_pv_kw += pv.existing_kw
            if !(s.site.roof_squarefeet === nothing)
                roof_max_kw = s.site.roof_squarefeet * pv.kw_per_square_foot
                beyond_existing_kw = min(roof_max_kw, beyond_existing_kw)
                pv_roof_limited = true
            end

        elseif pv.location == "ground"
            ground_existing_pv_kw += pv.existing_kw
            if !(s.site.land_acres === nothing)
                land_max_kw = s.site.land_acres / pv.acres_per_kw
                beyond_existing_kw = min(land_max_kw, beyond_existing_kw)
                pv_ground_limited = true
            end
        end

        existing_sizes[pv.name] = pv.existing_kw
        min_sizes[pv.name] = pv.existing_kw + pv.min_kw
        max_sizes[pv.name] = pv.existing_kw + beyond_existing_kw

        cost_slope, cost_curve_bp_x, cost_yint, n_segments = cost_curve(pv, s.financial)
        cap_cost_slope[pv.name] = cost_slope[1]
        if n_segments > 1
            cap_cost_slope[pv.name] = cost_slope
            push!(segmented_techs, pv.name)
            seg_max_size[pv.name] = Dict{Int,Float64}()
            seg_min_size[pv.name] = Dict{Int,Float64}()
            n_segs_by_tech[pv.name] = n_segments
            seg_yint[pv.name] = Dict{Int,Float64}()
            for s in 1:n_segments
                seg_min_size[pv.name][s] = cost_curve_bp_x[s]
                seg_max_size[pv.name][s] = cost_curve_bp_x[s+1]
                seg_yint[pv.name][s] = cost_yint[s]
            end
        end

        om_cost_per_kw[pv.name] = pv.om_cost_per_kw
        fillin_techs_by_exportbin(techs_by_exportbin, pv, pv.name)
    end

    if pv_roof_limited
        maxsize_pv_locations[:roof] = float(roof_existing_pv_kw + roof_max_kw)
    end
    if pv_ground_limited
        maxsize_pv_locations[:ground] = float(ground_existing_pv_kw + land_max_kw)
    end
    if pv_space_limited
        maxsize_pv_locations[:both] = float(both_existing_pv_kw + roof_max_kw + land_max_kw)
    end

    return nothing
end

function setup_wind_inputs(s::Scenario, max_sizes, min_sizes, existing_sizes,
    cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin)
    # TODO add incentives to Wind and use cost_curve function
    max_sizes["Wind"] = s.wind.max_kw
    min_sizes["Wind"] = s.wind.min_kw
    existing_sizes["Wind"] = 0.0
    cap_cost_slope["Wind"] = s.wind.installed_cost_per_kw
    om_cost_per_kw["Wind"] = s.wind.om_cost_per_kw
    production_factor["Wind", :] = prodfactor(s.wind, s.site.latitude, s.site.longitude, s.settings.time_steps_per_hour)
    fillin_techs_by_exportbin(techs_by_exportbin, s.wind, "Wind")
    return nothing
end

function setup_gen_inputs(s::Scenario, max_sizes, min_sizes, existing_sizes,
    cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin)
    # TODO add incentives to Generator and use cost_curve function
    max_sizes["Generator"] = s.generator.max_kw
    min_sizes["Generator"] = s.generator.existing_kw + s.generator.min_kw
    existing_sizes["Generator"] = s.generator.existing_kw
    cap_cost_slope["Generator"] = s.generator.installed_cost_per_kw
    om_cost_per_kw["Generator"] = s.generator.om_cost_per_kw
    production_factor["Generator", :] = prodfactor(s.generator)
    fillin_techs_by_exportbin(techs_by_exportbin, s.generator, "Generator")
    return nothing
end

function setup_present_worth_factors(s::Scenario, techs::Array{String, 1}, pvtechs::Array{String, 1})

    lvl_factor = Dict(t => 1.0 for t in techs)  # default levelization_factor of 1.0
    for (i, tech) in enumerate(pvtechs)  # replace 1.0 with actual PV levelization_factor (only tech with degradation)
        lvl_factor[tech] = levelization_factor(
            s.financial.analysis_years,
            s.financial.elec_cost_escalation_pct,
            s.financial.offtaker_discount_pct,
            s.pvs[i].degradation_pct  # TODO generalize for any tech (not just pvs)
        )
    end

    pwf_e = annuity(
        s.financial.analysis_years,
        s.financial.elec_cost_escalation_pct,
        s.financial.offtaker_discount_pct
    )

    pwf_om = annuity(
        s.financial.analysis_years,
        s.financial.om_cost_escalation_pct,
        s.financial.owner_discount_pct
    )

    if s.financial.third_party_ownership
        pwf_offtaker = annuity(s.financial.analysis_years, 0.0, s.financial.offtaker_discount_pct)
        pwf_owner = annuity(s.financial.analysis_years, 0.0, s.financial.owner_discount_pct)
        two_party_factor = (pwf_offtaker * (1 - s.financial.offtaker_tax_pct)) /
                           (pwf_owner * (1 - s.financial.owner_tax_pct))
    else
        two_party_factor = 1.0
    end

    return lvl_factor, pwf_e, pwf_om, two_party_factor
end

function setup_electric_utility_inputs(s::AbstractScenario)
    if s.electric_utility.outage_end_timestep > 0 &&
            s.electric_utility.outage_end_timestep > s.electric_utility.outage_start_timestep
        time_steps_without_grid = Int[i for i in range(s.electric_utility.outage_start_timestep,
                                                    stop=s.electric_utility.outage_end_timestep)]
        if s.electric_utility.outage_start_timestep > 1
            time_steps_with_grid = append!(
                Int[i for i in range(1, stop=s.electric_utility.outage_start_timestep - 1)],
                Int[i for i in range(s.electric_utility.outage_end_timestep + 1,
                                     stop=length(s.electric_load.loads_kw))]
            )
        else
            time_steps_with_grid = Int[i for i in range(s.electric_utility.outage_end_timestep + 1,
                                       stop=length(s.electric_load.loads_kw))]
        end
    else
        time_steps_without_grid = Int[]
        time_steps_with_grid = Int[i for i in range(1, stop=length(s.electric_load.loads_kw))]
    end
    return time_steps_with_grid, time_steps_without_grid
end

function adjust_load_profile(s::Scenario, production_factor::DenseAxisArray)
    if s.electric_load.loads_kw_is_net
        for pv in s.pvs if pv.existing_kw > 0
            s.electric_load.loads_kw .+= pv.existing_kw * production_factor[pv.name, :].data
        end end
    end

    if s.electric_load.critical_loads_kw_is_net
        for pv in s.pvs if pv.existing_kw > 0
            s.electric_load.critical_loads_kw .+= pv.existing_kw * production_factor[pv.name, :].data
        end end
    end
end

"""
    production_incentives(tech::AbstractTech, financial::Financial)

Intermediate function for building the PBI arrays in REoptInputs
"""
function production_incentives(tech::AbstractTech, financial::Financial)
    pwf_prod_incent = 0.0
    max_prod_incent = 0.0
    max_size_for_prod_incent = 0.0
    production_incentive_rate = 0.0
    T = typeof(tech)
    # TODO should Generator be excluded? (v1 has the PBI inputs for Generator)
    if !(nameof(T) in [:Generator, :Boiler, :Elecchl, :Absorpchl])
        if :degradation_pct in fieldnames(T)  # PV has degradation
            pwf_prod_incent = annuity_escalation(tech.production_incentive_years, -1*tech.degradation_pct,
                                                 financial.owner_discount_pct)
        else
            # prod incentives have zero escalation rate
            pwf_prod_incent = annuity(tech.production_incentive_years, 0, financial.owner_discount_pct)
        end
        max_prod_incent = tech.production_incentive_max_benefit
        max_size_for_prod_incent = tech.production_incentive_max_kw
        production_incentive_rate = tech.production_incentive_per_kwh
    end

    return pwf_prod_incent, max_prod_incent, max_size_for_prod_incent, production_incentive_rate
end

function fillin_techs_by_exportbin(techs_by_exportbin::Dict, tech::AbstractTech, tech_name::String)
    if tech.can_net_meter && :NEM in keys(techs_by_exportbin)
        push!(techs_by_exportbin[:NEM], tech_name)
        if tech.can_export_beyond_nem_limit && :EXC in keys(techs_by_exportbin)
            push!(techs_by_exportbin[:EXC], tech_name)
        end
    end

    if tech.can_wholesale && :WHL in keys(techs_by_exportbin)
        push!(techs_by_exportbin[:WHL], tech_name)
    end
    return nothing
end

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite


text:
function from_file(filename, sep = '\n')
    s = open(filename) do file
        read(file, String)
    end
    str_sudokus = split(strip(s), sep)
    grids = AbstractArray[]
    for str_sudoku in str_sudokus
        str_sudoku = replace(str_sudoku, "." => "0")
        one_line_grid = parse.(Int, split(str_sudoku, ""))
        nvals = length(one_line_grid)
        side_len = isqrt(nvals)
        grid = reshape(one_line_grid, side_len, side_len)
        push!(grids, grid)
    end
    return grids
end

function solve_sudoku(grid)
    m = CS.Optimizer(logging = [])
    side_len = size(grid, 1)
    block_size = isqrt(side_len)

    x = [[MOI.add_constrained_variable(m, MOI.Integer()) for i in 1:side_len] for j in 1:side_len]
    for r in 1:side_len, c in 1:side_len
        MOI.add_constraint(m, x[r][c][1], MOI.GreaterThan(1.0))
        MOI.add_constraint(m, x[r][c][1], MOI.LessThan(convert(Float64, side_len)))
    end

    # set variables
    for r in 1:side_len, c in 1:side_len
        if grid[r, c] != 0
            sat = [MOI.ScalarAffineTerm(1.0, x[r][c][1])]
            MOI.add_constraint(
                m,
                MOI.ScalarAffineFunction{Float64}(sat, 0.0),
                MOI.EqualTo(convert(Float64, grid[r, c])),
            )
        end
    end
    # sudoku constraints
    for r in 1:side_len
        MOI.add_constraint(
            m,
            MOI.VectorOfVariables([x[r][c][1] for c in 1:side_len]),
            CS.CPE.AllDifferent(side_len),
        )
    end
    for c in 1:side_len
        MOI.add_constraint(
            m,
            MOI.VectorOfVariables([x[r][c][1] for r in 1:side_len]),
            CS.CPE.AllDifferent(side_len),
        )
    end
    variables = [MOI.VariableIndex(0) for _ in 1:side_len]
    for br in 0:block_size-1
        for bc in 0:block_size-1
            variables_i = 1
            for i in (br * block_size + 1):((br + 1) * block_size), j in (bc * block_size + 1):((bc + 1) * block_size)
                variables[variables_i] = x[i][j][1]
                variables_i += 1
            end
            MOI.add_constraint(m, variables, CS.CPE.AllDifferent(side_len))
        end
    end

    MOI.optimize!(m)
    @assert MOI.get(m, MOI.TerminationStatus()) == MOI.OPTIMAL
end


text:
export FastNAU

struct FastNAU{F} <: DiffEqFlux.FastLayer
    in::Int
    out::Int
    initial_params::F
    function FastNAU(in::Int, out::Int; init=Flux.glorot_uniform)
        initial_params() = vec(init(out,in))
        new{typeof(initial_params)}(in, out, initial_params)
    end
end

DiffEqFlux.paramlength(f::FastNAU) = f.out * f.in
DiffEqFlux.initial_params(f::FastNAU) = f.initial_params()

(f::FastNAU)(x,p) = reshape(p, f.out, f.in) * x

Zygote.@adjoint function (f::FastNAU)(x,p)
    W = reshape(p, f.out, f.in)
    f(x,p), ȳ -> (nothing, W'ȳ, vec(ȳ * x'))
end

Base.show(io::IO, l::FastNAU) = print(io,"FastNAU(in=$(l.in), out=$(l.out))")


text:
DeveloperBuilder(name::String) = (name, 700)
SpecializedDeveloperBuilder(name::String, language::String) =
        tuplejoin(DeveloperBuilder(name), language)

@testset "Constructors" begin
    junior = Developer("Junior")

    @test junior.name == "Junior"
    @test junior.salary == 700

    artur = SpecializedDeveloper("Artur", "Julia")

    @test artur.name == "Artur"
    @test artur.salary == 700
    @test artur.language == "Julia"
end


text:
function SDPPrimal_cplx2real(sdp::SDPPrimal{T}) where T<:Complex

    block_to_vartype = Dict{String, Symbol}()
    sdpblocks = Dict{Tuple{Moment, String, Exponent, Exponent}, Float64}()
    sdplinsym = Dict{Tuple{Moment, String, Exponent}, Float64}()
    sdplin = Dict{Tuple{Moment, Exponent}, Float64}()
    sdpcst = Dict{Moment, Float64}()

    matrix_terms = Dict{String, Set{Tuple{Exponent, Exponent}}}()

    ## Complex blocks to real
    for ((moment, block_name, γ, δ), coeff) in sdp.blocks
        ctr_re, ctr_im = cplx2real_sdpctr(moment)
        γ_re, γ_im = cplx2real_sdpctr(γ)
        δ_re, δ_im = cplx2real_sdpctr(δ)

        @assert ctr_re != ctr_im

        println("$moment  is $ctr_re, $ctr_im")

        !haskey(matrix_terms, block_name) && (matrix_terms[block_name] = Set{Tuple{Exponent, Exponent}}())
        push!(matrix_terms[block_name], (min(γ, δ), max(γ, δ)))

        # Convert complex linear term to real ones
        ## Real part
        sdpblocks[(ctr_re, block_name, γ_re, δ_re)] = real(coeff)
        sdpblocks[(ctr_re, block_name, γ_im, δ_re)] = -imag(coeff)

        ## Imag part
        if product(moment.conj_part, moment.expl_part) != Exponent()

            sdpblocks[(ctr_im, block_name, γ_re, δ_re)] = imag(coeff)
            sdpblocks[(ctr_im, block_name, γ_im, δ_re)] = real(coeff)
        end
    end

    for (block_name, coords) in matrix_terms
        info(LOGGER, block_name)
        for (γ, δ) in coords
            @show γ, δ
            γ_re, γ_im = cplx2real_sdpctr(γ)
            δ_re, δ_im = cplx2real_sdpctr(δ)

            Xi_ctrname_re = get_Xictrname_re(block_name, γ, δ)

            sdpblocks[(Xi_ctrname_re, block_name, γ_re, δ_re)] =  1
            sdpblocks[(Xi_ctrname_re, block_name, γ_im, δ_im)] = -1
            info(LOGGER, "sdpblocks[($Xi_ctrname_re, $block_name, $γ_re, $δ_re)] =  1")
            info(LOGGER, "sdpblocks[($Xi_ctrname_re, $block_name, $γ_im, $δ_im)] = -1")

            Xi_ctrname_im = get_Xictrname_im(block_name, γ, δ)

            sdpblocks[(Xi_ctrname_im, block_name, γ_im, δ_re)] = 1
            sdpblocks[(Xi_ctrname_im, block_name, δ_im, γ_re)] = 1
            info(LOGGER, "sdpblocks[($Xi_ctrname_im, $block_name, $γ_im, $δ_re)] = 1")
            info(LOGGER, "sdpblocks[($Xi_ctrname_im, $block_name, $δ_im, $γ_re)] = 1")

        end
    end

    ## Complex symetric blocks to real
    for ((moment, block_name, var), coeff) in sdp.linsym
        ctr_re, ctr_im = cplx2real_sdpctr(moment)
        var_re, var_im = cplx2real_sdpctr(var)

        sdplinsym[(ctr_re, block_name, var_re)] = real(coeff)
        sdplinsym[(ctr_re, block_name, var_im)] = -imag(coeff)

        sdplinsym[(ctr_im, block_name, var_re)] = imag(coeff)
        sdplinsym[(ctr_im, block_name, var_im)] = real(coeff)
    end

    ## Complex vectors to real
    for ((moment, var), coeff) in sdp.lin
        ctr_re, ctr_im = cplx2real_sdpctr(moment)
        var_re, var_im = cplx2real_sdpctr(var)

        sdplin[(ctr_re, var_re)] = real(coeff)
        sdplin[(ctr_re, var_im)] = -imag(coeff)

        sdplin[(ctr_im, var_re)] = imag(coeff)
        sdplin[(ctr_im, var_im)] = real(coeff)
    end

    ## Complex coeff to real
    for (moment, coeff) in sdp.cst
        ctr_re, ctr_im = cplx2real_sdpctr(moment)

        sdpcst[ctr_re] = real(coeff)
        sdpcst[ctr_im] = imag(coeff)
    end

    ## Convert variable types
    for (block, vartype) in sdp.block_to_vartype
        if vartype == :SDPC
            block_to_vartype[block] = :SDP
        elseif vartype == :SymC
            block_to_vartype[block] = :Sym
        else
            error(LOGGER, "SDPPrimal_cplx2real(): Unhandled matrix type $vartype for $block")
        end
    end

    return SDPPrimal(block_to_vartype, sdpblocks, sdplinsym, sdplin, sdpcst)
end

function cplx2real_sdpctr(moment::Moment)
    ctr_re = Moment(moment.conj_part, moment.expl_part, moment.clique*"_Re")
    ctr_im = Moment(moment.conj_part, moment.expl_part, moment.clique*"_Im")

    return ctr_re, ctr_im
end

function cplx2real_sdpctr(expo::Exponent)
    expo_expl = expo

    (expo.degree.conjvar > 0) && (expo_expl = conj(expo))
    @assert expo_expl.degree.conjvar == 0

    expo_re, expo_im = Exponent(), Exponent()
    expo_re = Exponent(Variable(string(expo_expl, "_Re"), Real))
    expo_im = Exponent(Variable(string(expo_expl, "_Im"), Real))
    return expo_re, expo_im
end

function get_Xictrname_re(block_name::String, γ::Exponent, δ::Exponent)
    δ_ctr = δ
    if product(γ, δ) == Exponent()
        δ_ctr = Exponent(Variable("1_ctr", Real))
    end
    @show Moment(γ, δ_ctr, block_name*"_ReCtr")
    return Moment(γ, δ_ctr, block_name*"_ReCtr")
    # return Exponent(Variable(block_name*"_Re", Real)), Exponent(Variable(string(γ, "_", δ), Real))
end

function get_Xictrname_im(block_name::String, γ::Exponent, δ::Exponent)
    δ_ctr = δ
    if product(γ, δ) == Exponent()
        δ_ctr = Exponent(Variable("1_ctr", Real))
    end
    @show Moment(γ, δ_ctr, block_name*"_ImCtr")
    return Moment(γ, δ_ctr, block_name*"_ImCtr")
    # return Exponent(Variable(block_name*"_Im", Real)), Exponent(Variable(string(γ, "_", δ), Real))
end


text:
using GLExperiments
using Documenter

makedocs(;
    modules=[GLExperiments],
    authors="Bradley Lyman <lyman.brad3211@gmail.com> and contributors",
    repo="https://github.com/BradLyman/GLExperiments.jl/blob/{commit}{path}#L{line}",
    sitename="GLExperiments.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://BradLyman.github.io/GLExperiments.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/BradLyman/GLExperiments.jl",
    devbranch="main"
)


text:
module FromPython
using PyCall

include("ast.jl")

export transpile

const METHODS = Dict()

function register!(k, v)
    METHODS[k] = v
end

function transpile(ctx, code)
    ast = pyimport("ast")
    py_ast = ast.parse(code)
    AST.dispatch(ctx, py_ast)
end

function trace()
end

end