text:
#
# Warn about the order of the ITensor after contractions
#

const default_warn_order = 14

const warn_order =
  Ref{Union{Int, Nothing}}(default_warn_order)

"""
    get_warn_order()

Return the threshold for the order of an ITensor above which
ITensors will emit a warning.

You can set the threshold with the function `set_warn_order!(N::Int)`.
"""
get_warn_order() = warn_order[]

"""
    set_warn_order!(N::Int)

After this is called, ITensor will warn about ITensor contractions
that result in ITensors above the order `N`.

This function returns the initial warning threshold (what it was
set to before this function was called).

You can get the current threshold with the function `get_warn_order(N::Int)`. You can reset to the default value with
`reset_warn_order!()`.
"""
function set_warn_order!(N::Union{Int, Nothing})
  N_init = get_warn_order()
  warn_order[] = N
  return N_init
end

"""
    reset_warn_order!()

After this is called, ITensor will warn about ITensor contractions
that result in ITensors above the default order
$default_warn_order.

This function returns the initial warning threshold (what it was
set to before this function was called).
"""
reset_warn_order!() =
  set_warn_order!(default_warn_order)

"""
    disable_warn_order!()

After this is called, ITensor will not warn about ITensor
contractions that result in large ITensor orders.

This function returns the initial warning threshold (what it was
set to before this function was called).
"""
disable_warn_order!() =
  set_warn_order!(nothing)

"""
    @disable_warn_order

Disable warning about the ITensor order in a block of code.

# Examples
```julia
@disable_warn_order A * B
```
"""
macro disable_warn_order(block)
  quote
    local old_order = disable_warn_order!()
    r = $(esc(block))
    set_warn_order!(old_order)
    r
  end
end

"""
    @set_warn_order

Temporarily set the order threshold for warning about the ITensor
order in a block of code.

# Examples
```julia
@set_warn_order 12 A * B

@set_warn_order 15 begin
  C = A * B
  E = C * D
end
```
"""
macro set_warn_order(new_order, block)
  quote
    local old_order = set_warn_order!($(esc(new_order)))
    r = $(esc(block))
    set_warn_order!(old_order)
    r
  end
end

"""
    @reset_warn_order

Temporarily reset the order threshold for warning about the ITensor
order in a block of code to the default value $default_warn_order.

# Examples
```julia
@reset_warn_order A * B
```
"""
macro reset_warn_order(block)
  quote
    local old_order = reset_warn_order!()
    r = $(esc(block))
    set_warn_order!(old_order)
    r
  end
end

#
# A global timer used with TimerOutputs.jl
#

using NDTensors: timer

#
# Turn enable or disable combining QN ITensors before contracting
#

const _use_combine_contract = Ref(false)

use_combine_contract() = _use_combine_contract[]

function enable_combine_contract!()
  _use_combine_contract[] = true
  return nothing
end

function disable_combine_contract!()
  _use_combine_contract[] = false
  return nothing
end


text:
"""
For an initial guess from the index of a Hopf bifurcation point located in ContResult.bifpoint, returns a point which will be refined using `newtonHopf`.
"""
function HopfPoint(br::ContResult, index::Int64)
	@assert br.bifpoint[index].type == :hopf "The provided index does not refer to a Hopf point"
	bifpoint = br.bifpoint[index]							# Hopf point
	eigRes   = br.eig										# eigenvector at the Hopf point
	p = bifpoint.param										# parameter value at the Hopf point
	ω = abs(imag(eigRes[bifpoint.idx].eigenvals[bifpoint.ind_bif]))	# frequency at the Hopf point
	return BorderedArray(bifpoint.u, [p, ω] )
end

@with_kw struct HopfProblemMinimallyAugmented{TF, TJ, TJa, Td2f, vectype, S <: AbstractLinearSolver, Sbd <: AbstractBorderedLinearSolver, Sbda <: AbstractBorderedLinearSolver}
	F::TF 								# Function F(x, p) = 0
	J::TJ 								# Jacobian of F wrt x
	Jadjoint::TJa						# Adjoint of the Jacobian of F
	d2F::Td2f = nothing					# Hessian of F
	a::vectype							# close to null vector of (J - iω I)^*
	b::vectype							# close to null vector of J - iω I
	linsolver::S						# linear solver
	linbdsolver::Sbd					# linear bordered solver
	linbdsolverAdjoint::Sbda			# linear bordered solver for the jacobian adjoint
end

hasHessian(pb::HopfProblemMinimallyAugmented{TF, TJ, TJa, Td2f, vectype, S, Sbd, Sbda}) where {TF, TJ, TJa, Td2f, vectype, S, Sbd, Sbda} = Td2f !== Nothing

HopfProblemMinimallyAugmented(F, J, Ja, a, b, linsolve) = HopfProblemMinimallyAugmented(F, J, Ja, nothing, a, b, linsolve, BorderingBLS(linsolve), BorderingBLS(linsolve))

HopfProblemMinimallyAugmented(F, J, Ja, d2F, a, b, linsolve) = HopfProblemMinimallyAugmented(F, J, Ja, d2F, a, b, linsolve, BorderingBLS(linsolve), BorderingBLS(linsolve))

function (fp::HopfProblemMinimallyAugmented)(x, p::T, ω::T) where {T}
	# These are the equations of the minimally augmented (MA) formulation of the Hopf bifurcation point
	# input:
	# - x guess for the point at which the jacobian has a purely imaginary eigenvalue
	# - p guess for the parameter for which the jacobian has a purely imaginary eigenvalue
	# The jacobian of the MA problem is solved with a bordering method
	a = fp.a
	b = fp.b

	# we solve (J+iω)v + a σ1 = 0 with <b, v> = n
	n = T(1)
	v, σ1, flag, it = fp.linbdsolver(fp.J(x, p),
							a, b,
							T(0), zero(x), n; shift = Complex{T}(0, ω))

	# we solve (J+iω)'w + b σ2 = 0 with <a, w> = n
	# we find sigma2 = conj(sigma1)
	# w, σ2, _ = fp.linbdsolver(fp.Jadjoint(x, p) - Complex(0, ω) * I, b, a, 0., zeros(N), n)

	# the constraint is σ = <w, Jv> / n
	# σ = -dot(w, apply(fp.J(x, p) + Complex(0, ω) * I, v)) / n
	# we should have σ = σ1

	return fp.F(x, p), real(σ1), imag(σ1)
end

function (hopfpb::HopfProblemMinimallyAugmented{TF, TJ, TJa, Td2f, vectypeC, S, Sbd, Sbda})(x::BorderedArray{vectypeR, T}) where {TF, TJ, TJa, Td2f, vectypeC, vectypeR, S, Sbd, Sbda, T}
	res = hopfpb(x.u, x.p[1], x.p[2])
	return BorderedArray(res[1], [res[2], res[3]])
end

# Struct to invert the jacobian of the Hopf MA problem.
struct HopfLinearSolveMinAug <: AbstractLinearSolver; end

"""
The function solve the linear problem associated with a linearization of the minimally augmented formulation of the Hopf bifurcation point. The keyword `debug_` is used to debug the routine by returning several key quantities.
"""
function hopfMALinearSolver(x, p::T, ω::T, pbMA::HopfProblemMinimallyAugmented,
	 						duu, dup, duω;
							debug_ = false) where T
	# N = length(du) - 2
	# The jacobian should be passed as a tuple as Jac_hopf_MA(u0, pb::HopfProblemMinimallyAugmented) = (return (u0, pb, d2F::Bool))
	# The Jacobian J of the vector field is expressed at (x, p)
	# the jacobian expression of the hopf problem Jhopf is
	#					[ J dpF   0
	#					 σx  σp  σω]
	########## Resolution of the bordered linear system ########
	# J * dX	  + dpF * dp		   = du => dX = x1 - dp * x2
	# The second equation
	#	<σx, dX> +  σp * dp + σω * dω = du[end-1:end]
	# thus becomes
	#   (σp - <σx, x2>) * dp + σω * dω = du[end-1:end] - <σx, x1>
	# This 2x2 system is then solved to get (dp, dω)
	############### Extraction of function names #################
	Fhandle = pbMA.F
	J = pbMA.J
	Jadjoint = pbMA.Jadjoint
	d2F = pbMA.d2F
	a = pbMA.a
	b = pbMA.b

	# we define the following jacobian. It is used at least 3 times below. This avoid doing 3 times the possibly costly building of J(x, p)
	J_at_xp = J(x, p)

	δ = T(1e-9)
	ϵ1, ϵ2, ϵ3 = T(δ), T(δ), T(δ)

	# we solve Jv + a σ1 = 0 with <b, v> = n
	n = T(1)
	v, σ1, _, _ = pbMA.linbdsolver(J_at_xp, a, b, T(0), zero(x), n; shift = Complex{T}(0, ω))
	w, σ2, _, _ = pbMA.linbdsolverAdjoint(Jadjoint(x, p), b, a, T(0), zero(x), n; shift = -Complex{T}(0, ω))

	################### computation of σx σp ####################
	dpF   = (Fhandle(x, p + ϵ1)	 - Fhandle(x, p - ϵ1)) / T(2ϵ1)
	dJvdp = (apply(J(x, p + ϵ3), v) - apply(J(x, p - ϵ3), v)) / T(2ϵ3)
	σp = -dot(w, dJvdp) / n

	# case of sigma_omega
	σω = -dot(w, Complex{T}(0, 1) * v) / n

	x1, x2, _, (it1, it2) = pbMA.linsolver(J_at_xp, duu, dpF)

	# the case of ∂_xσ is a bit more involved
	# we first need to compute the value of ∂_xσ written σx
	# σx = zeros(Complex{T}, length(x))
	σx = similar(x, Complex{T})

	if hasHessian(pbMA) == false
		# We invert the jacobian of the Hopf problem when the Hessian of x -> F(x, p) is not known analytically. We thus rely on finite differences which can be slow for large dimensions
		prod(size(x)) > 1e4 && @warn "You might want to pass the Hessian, finite differences with $(prod(size(x))) unknowns"
		e = zero(x)
		for ii in CartesianIndices(x)
			e[ii] = T(1)
			d2Fve = (apply(J(x + ϵ2 * e, p), v) - apply(J(x - ϵ2 * e, p), v)) / T(2ϵ2)
			σx[ii] = -dot(w, d2Fve) / n
			e[ii] = T(0)
		end
		σxx1 = dot(σx, x1)
		σxx2 = dot(σx, x2)
	else
		d2Fv = d2F(x, p, v, x1)
		σxx1 = -dot(w, d2Fv) / n
		d2Fv = d2F(x, p, v, x2)
		σxx2 = -dot(w, d2Fv) / n
	end
	# we need to be carefull here because the dot produce conjugates. Hence the + dot(σx, x2) and + imag(dot(σx, x1) and not the opposite
	dp, dω = [real(σp - σxx2) real(σω);
			  imag(σp + σxx2) imag(σω) ] \
			  [dup - real(σxx1), duω + imag(σxx1)]

	if debug_
		return x1 - dp * x2, dp, dω, true, it1 + it2, (σx, σp, σω, dpF)
	else
		return x1 - dp * x2, dp, dω, true, it1 + it2
	end
end

function (hopfl::HopfLinearSolveMinAug)(Jhopf, du::BorderedArray{vectype, T}; debug_ = false)  where {vectype, T}
	out = hopfMALinearSolver((Jhopf[1]).u,
				(Jhopf[1]).p[1],
				(Jhopf[1]).p[2],
				Jhopf[2],
				du.u, du.p[1], du.p[2];
				debug_ = debug_)
	if debug_ == false
		return BorderedArray(out[1], [out[2], out[3]]), out[4], out[5]
	else
		return BorderedArray(out[1], [out[2], out[3]]), out[4], out[5], out[6]
	end
end

################################################################################################### Newton / Continuation functions
"""
This function turns an initial guess for a Hopf point into a solution to the Hopf problem based on a Minimally Augmented formulation. The arguments are as follows
- `F  = (x, p) -> F(x, p)` where `p` is the parameter associated to the Hopf point
- `J  = (x, p) -> d_xF(x, p)` associated jacobian
- `Jt = (x, p) -> transpose(d_xF(x, p))` associated jacobian transpose
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x,p)[v1, v2]`.
- `hopfpointguess` initial guess (x_0, p_0) for the Hopf point. It should a `BorderedArray` as given by the function HopfPoint.
- `eigenvec` guess for the  iω eigenvector
- `eigenvec_ad` guess for the -iω eigenvector
- `options::NewtonPar`
"""
function newtonHopf(F, J, Jt, d2F, hopfpointguess::BorderedArray{vectypeR, T}, eigenvec, eigenvec_ad, options::NewtonPar; normN = norm) where {vectypeR, T}
	hopfvariable = HopfProblemMinimallyAugmented(
		(x, p) -> F(x, p),
		(x, p) -> J(x, p),
		(x, p) -> Jt(x, p),
		d2F,
		copy(eigenvec),
		copy(eigenvec_ad),
		options.linsolver)
	hopfPb = u -> hopfvariable(u)

	# Jacobian for the Hopf problem
	Jac_hopf_MA(u0, pb::HopfProblemMinimallyAugmented) = (return (u0, pb))

	# options for the Newton Solver
	opt_hopf = @set options.linsolver = HopfLinearSolveMinAug()

	# solve the hopf equations
	return newton(x ->  hopfPb(x),
				x -> Jac_hopf_MA(x, hopfvariable),
				hopfpointguess,
				opt_hopf, normN = normN)
end

"""
call when hessian is unknown, finite differences are then used
"""
newtonHopf(F, J, Jt, hopfpointguess::BorderedArray{vectype, T}, eigenvec, eigenvec_ad, options::NewtonPar; kwargs...) where {T,vectype} = newtonHopf(F, J, Jt, nothing, hopfpointguess, eigenvec, eigenvec_ad, options; kwargs...)

newtonHopf(F, J, hopfpointguess::BorderedArray{vectype, T}, eigenvec, eigenvec_ad, options::NewtonPar; kwargs...) where {T,vectype} = newtonHopf(F, J, (x, p) -> transpose(J(x, p)), hopfpointguess, eigenvec, eigenvec_ad, options; kwargs...)

"""
Simplified call to refine an initial guess for a Hopf point. More precisely, the call is as follows

	`newtonHopf(F, J, Jt, br::ContResult, index::Int64, options)`

or

	`newtonHopf(F, J, Jt, d2F, br::ContResult, index::Int64, options)`

when the Hessian d2F is known. The parameters are as usual except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine.

!!! warning "Eigenvectors"
    This simplified call has been written when the eigenvectors are organised in a 2d Array `evec` where `evec[:,2]` is the second eigenvector in the list.
"""
function newtonHopf(F, J, Jt, d2F, br::ContResult, ind_hopf::Int64, options::NewtonPar ; normN = norm)
	hopfpointguess = HopfPoint(br, ind_hopf)
	bifpt = br.bifpoint[ind_hopf]
	options.verbose && println("--> Newton Hopf, the eigenvalue considered here is ", br.eig[bifpt.idx].eigenvals[bifpt.ind_bif])
	@assert bifpt.idx == bifpt.step + 1 "Error, the bifurcation index does not refer to the correct step"
	eigenvec = geteigenvector(options.eigsolver ,br.eig[bifpt.idx].eigenvec, bifpt.ind_bif)
	eigenvec_ad = conj.(eigenvec)

	# solve the hopf equations
	outhopf, hist, flag =  newtonHopf(F, J, Jt, d2F, hopfpointguess, eigenvec_ad, eigenvec, options; normN = normN)
	return outhopf, hist, flag
end

newtonHopf(F, J, Jt, br::ContResult, ind_hopf::Int64, options::NewtonPar;kwargs...) =  newtonHopf(F, J, Jt, nothing, br, ind_hopf, options::NewtonPar ;kwargs...)

newtonHopf(F, J, br::ContResult, ind_hopf::Int64, options::NewtonPar; kwargs...) = newtonHopf(F, J, (x, p) -> transpose(J(x, p)), br, ind_hopf, options; kwargs...)

"""
codim 2 continuation of Hopf points. This function turns an initial guess for a Hopf point into a curve of Hopf points based on a Minimally Augmented formulation. The arguments are as follows
- `(x, p1, p2)-> F(x, p1, p2)` where `p` is the parameter associated to the hopf point
- `J = (x, p1, p2)-> d_xF(x, p1, p2)` associated jacobian
- `Jt = (x, p1, p2) -> transpose(d_xF(x, p1, p2))` associated jacobian transpose
- `d2F = p2 -> ((x, p1, v1, v2) -> d2F(x, p1, p2, v1, v2))` this is the hessian of `F` computed at `(x, p1, p2)` and evaluated at `(v1, v2)`.
- `hopfpointguess` initial guess (x_0, p1_0) for the Hopf point. It should be a `Vector` or a `BorderedArray`
- `p2` parameter p2 for which hopfpointguess is a good guess
- `eigenvec` guess for the iω eigenvector at p1_0
- `eigenvec_ad` guess for the -iω eigenvector at p1_0
- `options::NewtonPar`
"""
function continuationHopf(F, J, Jt, d2F, hopfpointguess::BorderedArray{vectype, Tb}, p2_0::T, eigenvec, eigenvec_ad, options_cont::ContinuationPar ; kwargs...) where {T,Tb,vectype}
	# Bad way it creates a struct for every p2?
	# Jacobian for the hopf problem
	Jac_hopf_MA(u0, pb) = (return (u0, pb))

	# options for the Newton Solver inheritated from the ones the user provided
	options_newton = options_cont.newtonOptions

	hopfvariable = p2 -> HopfProblemMinimallyAugmented(
		(x, p1) -> F(x, p1, p2),
		(x, p1) -> J(x, p1, p2),
		(x, p1) -> Jt(x, p1, p2),
		d2F(p2),
		copy(eigenvec),
		copy(eigenvec_ad),
		options_newton.linsolver)

	hopfPb = (u, p2) -> hopfvariable(p2)(u)

	opt_hopf_cont = @set options_cont.newtonOptions.linsolver = HopfLinearSolveMinAug()

	# solve the hopf equations
	return continuation((x, p2) -> hopfPb(x, p2),
		(x, p2) -> Jac_hopf_MA(x, hopfvariable(p2)),
		hopfpointguess, p2_0,
		opt_hopf_cont,
		plot = true,
		printSolution = (u, p) -> u.p[1],
		plotSolution = (x;kwargs...) -> (xlabel!("p2", subplot=1); ylabel!("p1", subplot=1)  ) ; kwargs...)
end

continuationHopf(F, J, Jt, hopfpointguess::BorderedArray{vectype, Tb}, p2_0::T, eigenvec, eigenvec_ad, options_cont::ContinuationPar ; kwargs...) where {T, Tb, vectype} = continuationHopf(F, J, Jt, p2 -> nothing, hopfpointguess, p2_0, eigenvec, eigenvec_ad, options_cont ; kwargs...)

continuationHopf(F, J, Jt, hopfpointguess::BorderedArray{vectype, Tb}, p2_0::T, eigenvec, eigenvec_ad, options_cont::ContinuationPar ; kwargs...) where {T, Tb, vectype} = continuationHopf(F, J, (x, p1, p2) -> transpose(J(x, p1, p2)), hopfpointguess, p2_0, eigenvec, eigenvec_ad, options_cont ; kwargs...)

"""
Simplified call for continuation of Hopf point. More precisely, the call is as follows `continuationHopf(F, J, Jt, d2F, br::ContResult, index::Int64, options)` where the parameters are as for `continuationHopf` except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine.

Simplified calls are also provided but at the cost of using finite differences.

!!! warning "Eigenvectors"
    This simplified call has been written when the eigenvectors are organised in a 2d Array `evec` where `evec[:,2]` is the second eigenvector in the list.
"""
function continuationHopf(F, J, Jt, d2F, br::ContResult, ind_hopf::Int64, p2_0::Real, options_cont::ContinuationPar ; kwargs...)
	hopfpointguess = HopfPoint(br, ind_hopf)
	bifpt = br.bifpoint[ind_hopf]
	eigenvec = geteigenvector(options_cont.newtonOptions.eigsolver ,br.eig[bifpt.idx].eigenvec, bifpt.ind_bif)
	eigenvec_ad = conj.(eigenvec)
	return continuationHopf(F, J, Jt, d2F, hopfpointguess, p2_0, eigenvec_ad, eigenvec, options_cont ; kwargs...)
end

continuationHopf(F, J, Jt, br::ContResult, ind_hopf::Int64, p2_0::Real, options_cont::ContinuationPar ; kwargs...) = continuationHopf(F, J, Jt, p2 -> nothing, br, ind_hopf, p2_0, options_cont ; kwargs...)

continuationHopf(F, J, br::ContResult, ind_hopf::Int64, p2_0::Real, options_cont::ContinuationPar ; kwargs...) = continuationHopf(F, J, (x, p1, p2) -> transpose(J(x, p1, p2)), br, ind_hopf, p2_0, options_cont ; kwargs...)

ext: jl


text:
version = readchomp(joinpath(@__DIR__, "..", "ag-grid.version"))

isdir(joinpath(@__DIR__, "ag-grid-$(version)")) || mkdir(joinpath(@__DIR__, "ag-grid-$(version)"))

ag_grid_base = joinpath(@__DIR__, "ag-grid-$(version)", "ag-grid.js")
isfile(ag_grid_base) || download("https://unpkg.com/ag-grid-community@$(version)/dist/ag-grid-community.min.noStyle.js", ag_grid_base)

ag_grid_base_style = joinpath(@__DIR__, "ag-grid-$(version)", "ag-grid.css")
isfile(ag_grid_base_style) || download("https://unpkg.com/ag-grid-community@$(version)/dist/styles/ag-grid.css", ag_grid_base_style)

ag_grid_light = joinpath(@__DIR__, "ag-grid-$(version)", "ag-grid-light.css")
isfile(ag_grid_light) || download("https://unpkg.com/ag-grid-community@$(version)/dist/styles/ag-theme-balham.css", ag_grid_light)

ag_grid_dark = joinpath(@__DIR__, "ag-grid-$(version)", "ag-grid-dark.css")
isfile(ag_grid_dark) || download("https://unpkg.com/ag-grid-community@$(version)/dist/styles/ag-theme-balham-dark.css", ag_grid_dark)


text:
using IndexTransforms
using Test

@testset "IndexTransforms.jl" begin
    # Write your tests here.
end


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.
# *********************************************************************************
function prodfactor(pv::PV, latitude::Real, longitude::Real; timeframe="hourly")

    if !(ismissing(pv.prod_factor_series_kw))
        return pv.prod_factor_series_kw
    end

    url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
        "&lat=", latitude , "&lon=", longitude, "&tilt=", pv.tilt,
        "&system_capacity=1", "&azimuth=", pv.azimuth, "&module_type=", pv.module_type,
        "&array_type=", pv.array_type, "&losses=", round(pv.losses*100, digits=3), "&dc_ac_ratio=", pv.dc_ac_ratio,
        "&gcr=", 0.4, "&inv_eff=", pv.inv_eff*100, "&timeframe=", timeframe, "&dataset=nsrdb",
        "&radius=", 100
    )

    try
        @info "Querying PVWatts for prodfactor with " pv.name
        r = HTTP.get(url)
        response = JSON.parse(String(r.body))
        if r.status != 200
            error("Bad response from PVWatts: $(response["errors"])")
            # julia does not get here even with status != 200 b/c it jumps ahead to CIDER/reopt/src/core/reopt_inputs.jl:114
            # and raises ArgumentError: indexed assignment with a single value to many locations is not supported; perhaps use broadcasting `.=` instead?
        end
        @info "PVWatts success."
        watts = get(response["outputs"], "ac", []) / 1000  # scale to 1 kW system (* 1 kW / 1000 W)

        return collect(watts)
    catch e
        return "Error occurred : $e"
    end
end

function prodfactor(g::AbstractGenerator; ts_per_hour::Int=1)
    return ones(8760 * ts_per_hour)
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
# This implements an efficient Gibbs sampler for the multivariate truncated
# normal distribution subject to linear constraints from:

# Y. Li and S. Ghosh (2015)
# "Efficient Sampling Methods for Truncated Multivariate Normal and Student-t
#  Distributions Subject to Linear Inequality Constraints"
# Journal of Statistical Theory & Practice
# doi: 10.1080/15598608.2014.996690
# http://dx.doi.org/10.1080/15598608.2014.996690

function GibbsTMVN(μ, Σ, c, d, R̃, NumSamples = 10000,
                   x_init = zeros(length(μ)))
  p = length(μ)
  m = size(R̃, 1)
  samples = zeros(p, NumSamples)
  Σhalf = chol(Σ)'
  Σhalfinv = Σhalf \ eye(m)
  R = R̃ * Σhalf
  a = c - R̃ * μ
  b = d - R̃ * μ
  x = Σhalfinv * (x_init - μ)
  Rx = R * x

  # Make sure starting value is feasible:
  if !(all(Rx .>= a) && all(Rx .<= b))
    x = getFeasiblePoint(a, b, R)
  end

  bounds = [-Inf, Inf]
  tmp_l = zeros(m)
  tmp_u = zeros(m)

  for iter = 1:NumSamples
    for i = 1:p
      bounds[1] = -Inf
      bounds[2] = Inf
      GetBounds!(bounds, tmp_l, tmp_u, a, b, R, x, i)
      Spl = TruncatedNormalSampler(bounds[1], bounds[2])
      x[i] = rand(Spl)
    end
    samples[:, iter] = x
  end

  for iter = 1:NumSamples
     samples[:, iter] = Σhalf * samples[:, iter] + μ
  end
  return samples
end

function getFeasiblePoint(a::Vector{Float64}, b::Vector{Float64}, R::Matrix{Float64})
  set_default_solver(SCSSolver(verbose=false))
  m = size(R, 2)
  x = Variable(m)
  # Convex solver doesn't work with Infinity constraints
  finiteidx_a = (1:m)[a .> -Inf]
  finiteidx_b = (1:m)[b .< Inf]
  problem = minimize(0)
  for i in finiteidx_a
    problem.constraints += R[i, :] * x >= a
  end
  for i in finiteidx_b
    problem.constraints += R[i, :] * x <= b
  end
  solve!(problem)
  return vec(x.value)
end

function GetBounds!(bounds::Vector{Float64}, tmp_l::Vector{Float64},
                    tmp_u::Vector{Float64}, a::Array{Float64},
                    b::Vector{Float64}, R::Matrix{Float64},
                    x::Vector{Float64}, i::Int)
  m = size(R, 1)
  p = size(R, 2)
  tmp_l[:] = 0.0
  tmp_u[:] = 0.0

  for j = 1:m
    for k = 1:m
      if k != i
        tmp_l[j] -= R[j, k] * x[k]
      end
    end
    tmp_u[j] = copy(tmp_l[j])
    tmp_l[j] += a[j]
    tmp_l[j] /= R[j, i]
    tmp_u[j] += b[j]
    tmp_u[j] /= R[j, i]
  end

  for j = 1:m
    if R[j, i] > 0
      bounds[1] = max(tmp_l[j], bounds[1])
      bounds[2] = min(tmp_u[j], bounds[2])
    elseif R[j, i] < 0
      bounds[1] = max(tmp_u[j], bounds[1])
      bounds[2] = min(tmp_l[j], bounds[2])
    end
  end
  nothing
end


text:
"""Insert the projectors computed from the boundary MPS"""

default_projector_center(tn::Matrix) = (:, (size(tn, 2) + 1) ÷ 2)

# Split the links of the 2D tensor network in preperation for
# inserting MPS projectors.
function split_network(
  tn::Matrix{ITensor}; projector_center=default_projector_center(tn), rotation=false
)
  tn = rotation ? rotr90(tn) : tn
  @assert (projector_center[1] == :)
  tn_split = copy(tn)
  nrows, ncols = size(tn)
  @assert length(projector_center) == 2
  @assert projector_center[1] == Colon()
  projector_center_cols = projector_center[2]
  for ncol in 1:ncols
    if ncol ∉ projector_center_cols
      tn_split[:, ncol] = split_links(tn_split[:, ncol])
    end
  end
  return rotation ? rotr90(tn_split, -1) : tn_split
end

# From an MPS, create a 1-site projector onto the MPS basis
function projector(x::MPS, projector_center)
  # Gauge the boundary MPS towards the projector_center column
  x = orthogonalize(x, projector_center)

  l = commoninds(get_itensor(x, projector_center - 1), x[projector_center])
  r = commoninds(get_itensor(x, projector_center + 1), x[projector_center])

  uₗ = x[1:(projector_center - 1)]
  uᵣ = reverse(x[(projector_center + 1):end])
  nₗ = length(uₗ)
  nᵣ = length(uᵣ)

  uₗᴴ = dag.(uₗ)
  uᵣᴴ = dag.(uᵣ)

  uₗ′ = reverse(prime.(uₗ))
  uᵣ′ = reverse(prime.(uᵣ))

  if !isempty(uₗ′)
    uₗ′[1] = replaceinds(uₗ′[1], l' => l)
  end
  if !isempty(uᵣ′)
    uᵣ′[1] = replaceinds(uᵣ′[1], r' => r)
  end

  Pₗ = vcat(uₗᴴ, uₗ′)
  Pᵣ = vcat(uᵣᴴ, uᵣ′)
  return Pₗ, Pᵣ
end

function insert_projectors(
  tn::Matrix{ITensor},
  boundary_mps::Vector{MPS},
  dirs::BoundaryMPSDir"top_to_bottom";
  projector_center,
)
  nrows, ncols = size(tn)
  @assert length(boundary_mps) == nrows - 1
  @assert all(x -> length(x) == ncols, boundary_mps)
  @assert length(projector_center) == 2
  @assert (projector_center[1] == :)
  projector_center_cols = projector_center[2]
  projectors = [projector(x, projector_center_cols) for x in boundary_mps]
  projectors_left = first.(projectors)
  projectors_right = last.(projectors)
  tn_split = split_network(tn; projector_center=projector_center)
  return tn_split, projectors_left, projectors_right
end

function insert_projectors(
  tn::Matrix{ITensor}, boundary_mps::Vector{MPS}; dirs, projector_center
)
  return insert_projectors(
    tn, boundary_mps, BoundaryMPSDir(dirs); projector_center=projector_center
  )
end

function insert_projectors(tn, boundary_mps::BoundaryMPS; center, projector_center=nothing)
  top_bottom = true
  if (center[1] == :)
    center = reverse(center)
    tn = rotr90(tn)
    top_bottom = false
  end

  @assert (center[2] == :)
  center_row = center[1]

  if isnothing(projector_center)
    projector_center = default_projector_center(tn)
  end
  @assert (projector_center[1] == :)

  # Approximately contract the tensor network.
  # Outputs a Vector of boundary MPS.
  #boundary_mps_top = contract_approx(tn; alg="boundary_mps", dirs="top_to_bottom", cutoff=cutoff, maxdim=maxdim)
  #boundary_mps_bottom = contract_approx(tn; alg="boundary_mps", dirs="bottom_to_top", cutoff=cutoff, maxdim=maxdim)

  boundary_mps_top = top_bottom ? boundary_mps.top : boundary_mps.left
  boundary_mps_bottom = top_bottom ? boundary_mps.bottom : boundary_mps.right

  # Insert approximate projectors into rows of the network
  boundary_mps_tot = vcat(
    boundary_mps_top[1:(center_row - 1)], dag.(boundary_mps_bottom[center_row:end])
  )
  tn_split, projectors_left, projectors_right = insert_projectors(
    tn, boundary_mps_tot; dirs="top_to_bottom", projector_center=projector_center
  )

  if !top_bottom
    tn_split = rotl90(tn_split)
  end
  return tn_split, projectors_left, projectors_right
end

function contraction_cache_top(tn, boundary_mps::Vector{MPS}, n)
  tn_cache = fill(ITensor(1.0), size(tn))
  for nrow in 1:size(tn, 1)
    if nrow == n
      tn_cache[nrow, :] .= boundary_mps[nrow]
    elseif nrow > n
      tn_cache[nrow, :] = tn[nrow, :]
    end
  end
  return tn_cache
end

# Return a network that when contracted equals the
# squared norm of the input network.
# TODO: rename norm2_network
function sqnorm(ψ::Matrix{ITensor})
  # Square the tensor network
  ψᴴ = addtags(linkinds, dag.(ψ), "bra")
  ψ′ = addtags(linkinds, ψ, "ket")
  return mapreduce(vec, vcat, (ψᴴ, ψ′))
end

# Return a network that approximates the squared norm of the
# input network by inserting projectors determined from approximately
# contracting the network with boundary MPS.
# TODO: rename norm2_network_approx
function sqnorm_approx(ψ::Matrix{ITensor}; center, cutoff, maxdim)
  # Square the tensor network
  ψᴴ = addtags(linkinds, dag.(ψ), "bra")
  ψ′ = addtags(linkinds, ψ, "ket")
  # TODO: implement contract(commoninds, ψ′, ψᴴ)
  tn = ψ′ .* ψᴴ
  # Insert projectors horizontally (to measure e.g. properties
  # in a row of the network)
  bmps = boundary_mps(tn; cutoff=cutoff, maxdim=maxdim)
  tn_projected = insert_projectors(tn, bmps; center=center)
  tn_split, Pl, Pr = tn_projected
  ψᴴ_split = split_network(ψᴴ)
  ψ′_split = split_network(ψ′)
  Pl_flat = reduce(vcat, Pl)
  Pr_flat = reduce(vcat, Pr)
  return mapreduce(vec, vcat, (ψᴴ_split, ψ′_split, Pl_flat, Pr_flat))
end


text:
rng = StableRNGs.StableRNG(123)

table = load_iris()
y, X = unpack(table, ==(:target), _->true, rng=rng)

@testset_accelerated "regularization has an effect" accel begin

    model = MLJFlux.NeuralNetworkClassifier(acceleration=accel,
                                    builder=MLJFlux.Linear(),
                                    rng=rng)
    model2 = deepcopy(model)
    model3 = deepcopy(model)
    model3.lambda = 0.1

    e = evaluate(model, X, y, resampling=Holdout(), measure=LogLoss())
    loss1 = e.measurement[1]

    e = evaluate(model2, X, y, resampling=Holdout(), measure=LogLoss())
    loss2 = e.measurement[1]

    e = evaluate(model3, X, y, resampling=Holdout(), measure=LogLoss())
    loss3 = e.measurement[1]

    @test loss1 ≈ loss2
    @test !(loss2 ≈ loss3)
end


text:
function aero_wigner(a, b,   wigner_A, wigner_B)

    # Constants
    # μ  = 0.3
    # σ  = 6.82
    wl = 0.55
    FT = Float64
    # DiscreteUniform(a, b)
    size_distribution = Uniform(a, b)
    maxi = b

    # Generate aerosol:
    aero = Scattering.UnivariateAerosol(size_distribution, maxi, 5, 1.3, 0.0)
    r, wᵣ = Scattering.gauleg(aero.nquad_radius, a, b ; norm=true)
    wₓ = pdf.(aero.size_distribution, r)

    # pre multiply with wᵣ to get proper means eventually:
    wₓ .*= wᵣ

    # normalize (could apply a check whether cdf.(aero.size_distribution,r_max) is larger than 0.99:
    wₓ /= sum(wₓ)
    # @show wₓ, r
    greek_coefs = compute_B(aero, wigner_A, wigner_B, wl, r, wₓ)
    N_max = Scattering.get_n_max(2 * π * maxi / 0.55)
    n_mu = 5 * N_max - 1;

    μ, w_μ = Scattering.gausslegendre(n_mu)

    f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄ = Scattering.reconstruct_phase(greek_coefs, μ)
    # @show greek_coefs.β
    return f₁₁, f₁₂, f₂₂, f₃₃, μ
end

# plot(μ,f₁₁,yscale=:log10 )

p1 = plot(μ, f₁₁, yscale=:log10, title="f₁₁")

p2 = plot(μ, f₁₂, yscale=:log10, title="f₁₂")

p3 = plot(μ, f₂₂, yscale=:log10, title="f₂₂")

p4 = plot(μ, f₃₃, yscale=:log10, title="f₃₃")

plot(p1, p2, p3, p4,layout=(2, 2), legend=false)

@gif for a = 0.001:0.035:5
    # @show a
    f₁₁, f₁₂, f₂₂, f₃₃, μ =  aero_wigner(a, a + 0.05, wigner_A, wigner_B)
    # @show f₁₁[1]
    p1 = plot(μ, f₁₁, yscale=:log10, title="f₁₁", label="r(μm)=$a")
    ylims!(1e-3, 1e3)
    p2 = plot(μ, f₁₂ ./ f₁₁,  title="f₁₂/f₁₁", label="Q/I")
    ylims!(-1.1, 1.1)
    plot(p1, p2, layout=(2, 1))
end

anim = @animate for a = 0.001:0.01:1
    # @show a
    f₁₁, f₁₂, f₂₂, f₃₃, μ =  aero_wigner(a, a + 0.01, wigner_A, wigner_B)
    # @show f₁₁[1]
    p1 = plot(μ, f₁₁, yscale=:log10, title="f₁₁", label="r(μm)=$a")
    ylims!(1e-3, 1e3)
    p2 = plot(μ, f₁₂ ./ f₁₁,  title="f₁₂/f₁₁", label="Q/I")
    ylims!(-1.1, 1.1)
    plot(p1, p2, layout=(2, 1))
end


text:
"""
    rm2lspm(NUM, DEN; contr = false, obs = false, atol = 0, rtol) -> (A, B, C, D, blkdims)

Construct a representation of the rational matrix `R(λ) := NUM(λ)./DEN(λ)` as the sum of its strictly proper part,
realized as a structured linearization

              | A-λI | B |
              |------|---|
              |  C   | 0 |

and its polynomial part `D(λ)`. The resulting representation satisfies

     R(λ) = C*inv(λI-A)*B + D(λ).

The linearization `(A-λI,B,C,0)` of the strictly proper part is in general non-minimal.
A controllable realization is constructed if `contr = true`, while an observable realization results if `obs = true`.
The matrix `A` results block diagonal and the vector `blkdims` contains the sizes of the diagonal blocks of `A`,
such that the size of the `i`-th diagonal block is provided in `blkdims[i]`. See [1] for details.

`NUM(λ)` is a polynomial matrix of the form `NUM(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`, for which
the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `NUM`,
where `NUM[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

`DEN(λ)` is a polynomial matrix of the form `DEN(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`, for which
the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `DEN`,
where `DEN[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `NUM(λ)` and `DEN(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The polynomial part `D(λ)` is contained in `D`, which is a 2-dimensional array if `D(λ)` has degree 0
or a 3-dimensional array if `D(λ)` has degree  greater than 0. In the latter case `D[:,:,i]` contains the `i`-th
coefficient matrix multiplying `λ**(i-1)`.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances for the
nonzero coefficients of `DEN(λ)`, respectively.

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function rm2lspm(NUM::AbstractArray{T1,3},DEN::AbstractArray{T2,3}; contr::Bool = false, obs::Bool = false,
                 atol::Real = 0, rtol::Real = Base.rtoldefault(float(real(T1)))) where {T1,T2}
    p, m, knum = size(NUM)
    p1, m1, kden = size(DEN)
    (p,m) == (p1, m1) || error("Numerator and denominator polynomial matrices must have the same size")

    degQ1 = 0
    degD1 = 0
    n = 0
    for j = 1:m
        for i = 1:p
            n1 = poldeg1(NUM[i,j,1:knum])
            d1 = poldeg1(DEN[i,j,1:kden])
            d1 == 0 && error("DivideError: zero denominator polynomial")
            n1 < d1 || (degQ1 = max(degQ1,n1-d1+1))
            degD1 = max(degD1,d1)
            n += d1-1
        end
    end
    T = eltype(one(T1)/one(T2))
    degQ1 ≤ 1 ? D = zeros(T,p,m) : D = zeros(T,p,m,max(degQ1,1))
    A = zeros(T,n,n)
    B = zeros(T,n,m)
    C = zeros(T,p,n)
    ONE = one(T)
    contr && obs && (m >= p ? obs = false : contr = false)
    if contr
        k = 1
        nb = 0
        blkdims = zeros(Int,m)
        for j = 1:m
            d1 = poldeg1(DEN[1,j,1:kden])
            d1 == 1 ? mc = [one(T)] : (mc = DEN[1,j,1:d1])
            for i = 2:p
                d1 = poldeg1(DEN[i,j,1:kden])
                d1 > 1 && ( mc = pollcm(mc, DEN[i,j,1:d1], atol = atol, rtol = rtol) )
            end
            d1 = poldeg1(mc)
            n1 = d1-1
            i2 = k+n1-1
            if n1 > 0
                nb += 1
                blkdims[nb] = n1
                mc = mc/mc[d1]
                A[k, k:i2] = -mc[n1:-1:1]
                for ii = k:i2-1
                    A[ii+1, ii] = ONE
                end
                B[k,j] = ONE
                for i = 1:p
                    q, r = poldivrem(NUM[i,j,1:knum],DEN[i,j,1:kden])
                    nq1 = poldeg1(q)
                    nr1 = poldeg1(r)
                    degQ1 ≤ 1 ? D[i,j] = q[1] : D[i,j,1:nq1] = copy_oftype(q[1:nq1],T)
                    if nr1 > 0
                        pc = poldiv(conv(mc,r),DEN[i,j,1:kden])
                        npc = length(pc)
                        C[i,i2-npc+1:i2] = pc[npc:-1:1]
                    end
                end
                k += n1
            else
                for i = 1:p
                    nq1 = poldeg1(NUM[i,j,1:knum])
                    if nq1 > 0
                       q = NUM[i,j,1:nq1]/DEN[i,j,1]
                       degQ1 ≤ 1 ? D[i,j] = q[1] : D[i,j,1:nq1] = copy_oftype(q[1:nq1],T)
                    end
                end
            end
        end
        n = k-1
        return A[1:n,1:n], B[1:n,:], C[:,1:n], D, blkdims[1:nb]
    elseif obs
        k = 1
        nb = 0
        blkdims = zeros(Int,p)
        for i = 1:p
            d1 = poldeg1(DEN[i,1,1:kden])
            d1 == 1 ? mc = [one(T)] : (mc = DEN[i,1,1:d1])
            for j = 2:m
                d1 = poldeg1(DEN[i,j,1:kden])
                d1 > 1 && ( mc = pollcm(mc, DEN[i,j,1:d1], atol = atol, rtol = rtol) )
            end
            d1 = poldeg1(mc)
            n1 = d1-1
            i2 = k+n1-1
            if n1 > 0
                nb += 1
                blkdims[nb] = n1
                mc = mc/mc[d1]
                A[k:i2,i2] = -mc[1:n1]
                for ii = k:i2-1
                    A[ii+1,ii] = ONE
                end
                C[i,i2] = ONE
                for j = 1:m
                    q, r = poldivrem(NUM[i,j,1:knum],DEN[i,j,1:kden])
                    nq1 = poldeg1(q)
                    nr1 = poldeg1(r)
                    degQ1 ≤ 1 ? D[i,j] = q[1] : D[i,j,1:nq1] = copy_oftype(q[1:nq1],T)
                    if nr1 > 0
                        pc = poldiv(conv(mc,r),DEN[i,j,1:kden])
                        npc = length(pc)
                        B[k:k+npc-1,j] = pc[1:npc]
                    end
                end
                k += n1
            else
                for j = 1:m
                    nq1 = poldeg1(NUM[i,j,1:knum])
                    if nq1 > 0
                       q = NUM[i,j,1:nq1]/DEN[i,j,1]
                       degQ1 ≤ 1 ? D[i,j] = q[1] : D[i,j,1:nq1] = copy_oftype(q[1:nq1],T)
                    end
                end
            end
        end
        n = k-1
        return A[1:n,1:n], B[1:n,:], C[:,1:n], D, blkdims[1:nb]
    else
        k = 1
        blkdims = zeros(Int,m*p)
        nb = 0
        for j = 1:m
            for i = 1:p
                q, r = poldivrem(NUM[i,j,1:knum],DEN[i,j,1:kden])
                nq1 = poldeg1(q)
                nr1 = poldeg1(r)
                degQ1 ≤ 1 ? D[i,j] = q[1] : D[i,j,1:nq1] = copy_oftype(q[1:nq1],T)
                d1 = poldeg1(DEN[i,j,1:kden])
                n1 = d1-1
                i2 = k+n1-1
                if n1 > 0
                    nb += 1
                    blkdims[nb] = n1
                    s = -DEN[i,j,d1]
                    A[k, k:i2] = DEN[i,j,n1:-1:1]/s
                    for ii = k:i2-1
                        A[ii+1, ii] = ONE
                    end
                    B[k,j] = ONE
                    C[i,i2-nr1+1:i2] = r[nr1:-1:1]/(-s)
                    k += n1
                end
            end
        end
        return A, B, C, D, blkdims[1:nb]
    end
end
function rm2lspm(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}};kwargs...)
    return rm2lspm(poly2pm(NUM),poly2pm(DEN);kwargs...)
end
function rm2lspm(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}};kwargs...)
   return rm2lspm(poly2pm(NUM);kwargs...)
end
function rm2lspm(NUM::AbstractArray{T,3}; kwargs...) where T
    sys = rm2lspm(NUM,ones(T,size(NUM)[1:2]...,1); kwargs...)
    return sys[1:4]..., [0]
end

"""
    rmeval(NUM,DEN,val) -> Rval

Evaluate the rational matrix `R(λ) := NUM(λ)./DEN(λ)` for `λ = val`.

`NUM(λ)` is a polynomial matrix of the form `NUM(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`, for which
the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

`DEN(λ)` is a polynomial matrix of the form `DEN(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`, for which
the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `NUM(λ)` and `DEN(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.
"""
function rmeval(NUM::AbstractArray{T1,3},DEN::AbstractArray{T2,3},val::Number) where {T1,T2}
    return pmeval(NUM,val) ./ pmeval(DEN,val)
end
function rmeval(NUM::AbstractArray{T,3},val::Number) where {T}
    return pmeval(NUM,val)
end
rmeval(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
       DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},val::Number) =
       rmeval(poly2pm(NUM),poly2pm(DEN),val::Number)
rmeval(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},val::Number) =
       pmeval(poly2pm(NUM),val::Number)
"""
      rm2ls(NUM, DEN; contr = false, obs = false, noseig = false, minimal = false,
            fast = true, atol = 0, rtol) -> (A, E, B, C, D, blkdims)

Build a structured linearization as a system matrix `S(λ)` of the form

               | A-λE | B |
        S(λ) = |------|---|
               |  C   | D |

for the rational matrix `R(λ) = NUM(λ) ./ DEN(λ)`, such that `S(λ)` preserves a part of the Kronecker structure of `R(λ)`.
The regular pencil `A-λE` has the block diagonal form

             | Af-λI |   0   |
      A-λE = |-------|-------|
             |  0    | I-λEi |

and the dimensions of the diagonal blocks `Af` and `Ei` are provided in `blkdims[1]` and `blkdims[2]`, respectively.

`NUM(λ)` is a polynomial matrix of the form `NUM(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`, for which
the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `NUM`,
where `NUM[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

`DEN(λ)` is a polynomial matrix of the form `DEN(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`, for which
the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `DEN`,
where `DEN[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `NUM(λ)` and `DEN(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

If `n` is the order of `A-λE`, then the computed linearization satisfies:

(1) `A-λE` is regular and `R(λ) = C*inv(λE-A)*B+D`;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `obs = true`, in which case
the left Kronecker structure is preserved;

(4) `A-λE` has no non-dynamic modes if `minimal = true` or `noseig = true`.

If conditions (1)-(4) are satisfied, the linearization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the linearization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
For an irreducible linearization `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The descriptor system based linearization is built using the methods described in [1] in conjunction with
pencil manipulation algorithms [2] and [3] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances, respectively, for the
nonzero coefficients of `NUM(λ)` and `DEN(λ)`.

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[2] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[3] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function rm2ls(NUM::AbstractArray{T1,3},DEN::AbstractArray{T2,3}; minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
                fast::Bool = true, atol::Real = zero(real(T1)),
                rtol::Real = (min(size(NUM)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
    Ap, Bp, Cp, D = rm2lspm(NUM, DEN, contr = contr, obs = obs, atol = atol, rtol = rtol)
    minimal && (contr = true; obs = true; noseig = true;)
    if ndims(D) == 2
        A, B, C  = lsminreal(Ap, Bp, Cp, contr = contr, obs = obs, fast = fast, atol = atol, rtol = rtol)
        E = I
        nf = size(A,1)
        ni = 0
    else
        T = typeof(one(T1)/one(T2))
        Ai, Ei, Bi, Ci, Di  = pm2ls(D, minimal = minimal, contr = contr, obs = obs, fast = fast, atol = atol, rtol = rtol)
        Af, Ef, Bf, Cf, D  = lsminreal2(Ap,I,Bp,Cp,Di,contr = contr, obs = obs, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
        nf = size(Af,1)
        ni = size(Ai,1)
        A = [Af zeros(T,nf,ni); zeros(T,ni,nf) Ai]
        E = [Ef zeros(T,nf,ni); zeros(T,ni,nf) Ei]
        B = [Bf; Bi]
        C = [Cf Ci]
    end
    return A, E, B, C, D, [nf,ni]
end
rm2ls(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
      DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       rm2ls(poly2pm(NUM),poly2pm(DEN); kwargs...)
function rm2ls(NUM::AbstractArray{T,3}; kwargs...) where T
    sys = pm2ls(NUM; kwargs...)
    return sys..., [0;size(sys[1],1)]
end
function rm2ls(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    sys = pm2ls(poly2pm(NUM); kwargs...)
    return sys..., [0;size(sys[1],1)]
end
"""
    ls2rm(A, E, B, C, D; fast = true, atol1 = 0, atol2 = 0, gaintol = 0, rtol = min(atol1,atol2) > 0 ? 0 : n*ϵ, val) -> (NUM, DEN)

Build the rational matrix `R(λ) = C*inv(λE-A)*B+D := NUM(λ) ./ DEN(λ)` corresponding to its structured linearization

            | A-λE | B |
     S(λ) = |------|---|
            |  C   | D |

by explicitly determining for each rational entry of `R(λ)`, its numerator coefficients from its finite zeros and gain, and
its denominator coefficients from its finite poles (see [1]).

The keyword arguments `atol1` and `atol2` specify the absolute tolerances for the elements of `A`, `B`, `C`, `D`, and,
respectively, of `E`, and `rtol` specifies the relative tolerances for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `(n+1)*ϵ`, where `n` is the dimension of `A`, and `ϵ` is the
machine epsilon of the element type of coefficients of `A`.

The keyword argument `gaintol` specifies the threshold for the magnitude of the nonzero elements of the gain matrix
`C*inv(γE-A)*B+D`, where `γ = val` if `val` is a number or `γ` is a randomly chosen complex value of unit magnitude,
if `val = missing`. Generally, `val` should not be a zero of any of entries of `R`.

`NUM(λ)` is a polynomial matrix of the form `NUM(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`, for which
the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `NUM`,
where `NUM[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

`DEN(λ)` is a polynomial matrix of the form `DEN(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`, for which
the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `DEN`,
where `DEN[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

[1] A. Varga Computation of transfer function matrices of generalized state-space models.
    Int. J. Control, 50:2543–2561, 1989.
"""
function ls2rm(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}}, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix;
                fast::Bool = true, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                gaintol::Real = zero(real(eltype(A))), val::Union{Number,Missing} = missing,
                rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)))
   emat = (typeof(E) <: AbstractMatrix)
   eident = !emat || isequal(E,I)
   n = LinearAlgebra.checksquare(A)
   emat && (n,n) != size(E) && throw(DimensionMismatch("A and E must have the same dimensions"))
   p, m = size(D)
   (n,m) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions"))
   (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions"))

   T = promote_type(eltype(A),eltype(E), eltype(B),eltype(C),eltype(D))
   T <: BlasFloat || (T = promote_type(Float64,T))
   compl = T <: Complex

   ismissing(val) && (val = exp(rand()*im))
   Rval = lseval(A, E, B, C, D, val)
   unimodular = emat && isunimodular(A, E, atol1 = atol1, atol2 = atol2, rtol = rtol)
   NUM = zeros(T,p,m,n+1)
   if unimodular
      DEN = ones(T,p,m,1)
   else
      DEN = zeros(T,p,m,n+1)
      DEN[:,:,1] = ones(T,p,m)
   end
   for i = 1:p
       Ao, Eo, Bo, Co, Do = lsminreal2(A, E, B, C[i:i,:], D[i:i,:], contr = false, infinite = false, noseig = false,
                                       atol1 = atol1, atol2 = atol2, rtol = rtol)
       for j = 1:m
           if abs(Rval[i,j]) > gaintol
              A1, E1, B1, C1, D1 = lsminreal2(Ao, Eo, Bo[:,j:j], Co, Do[:,j:j], obs = false, infinite = false, noseig = false,
                                              atol1 = atol1, atol2 = atol2, rtol = rtol)
              zer, iz, _ = spzeros(A1, E1, B1, C1, D1; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
              cz, zval = polcoeffval(zer[1:(length(zer)-sum(iz))],val)
              if unimodular
                 NUM[i,j,1:length(cz)] = compl ? cz*(Rval[i,j]/zval) : real(cz*(Rval[i,j]/zval))
              else
                 pol, ip, _ = pzeros(A1, E1; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
                 cp, pval = polcoeffval(pol[1:(length(pol)-sum(ip))],val)
                 NUM[i,j,1:length(cz)] = compl ? cz*(Rval[i,j]*pval/zval) : real(cz*(Rval[i,j]*pval/zval))
                 DEN[i,j,1:length(cp)] = compl ? cp : real(cp)
              end
            end
       end
   end
   return NUM[:,:,1:pmdeg(NUM)+1], DEN[:,:,1:pmdeg(DEN)+1]
end
ls2rm(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix; kwargs...) =
      ls2rm(A::AbstractMatrix, I, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix; kwargs...)
"""
     rm2lps(NUM, DEN; contr = false, obs = false) -> (A, E, B, F, C, G, D, H, blkdims)

Build a structured pencil based linearization as a system matrix `S(λ)` of the form

            | A-λE | B-λF |
     S(λ) = |------|------|
            | C-λG | D-λH |

for the rational matrix `R(λ) = NUM(λ) ./ DEN(λ)` such that `S(λ)` preserves a part of the Kronecker structure of `R(λ)`.
The regular pencil `A-λE` has the block diagonal form

             | Af-λI |   0    |
      A-λE = |-------|--------|
             |  0    | Ai-λEi |

and the dimensions of the diagonal blocks `Af` and `Ai` are provided in `blkdims[1]` and `blkdims[2]`, respectively.

`NUM(λ)` is a polynomial matrix of the form `NUM(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`, for which
the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `NUM`,
where `NUM[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

`DEN(λ)` is a polynomial matrix of the form `DEN(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`, for which
the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `DEN`,
where `DEN[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `NUM(λ)` and `DEN(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

If `n` is the order of the matrix `A`, then the computed linearization satisfies:

(1) `A-λE` is a `n x n` regular pencil;

(2) `R(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH`;

(3) `rank[B-λF A-λE] = n` for any finite and infinite `λ` (strong controllability) if `contr = true`, in which case
the right Kronecker structure is preserved;

(4) `rank[A-λE; C-λG] = n` for any finite and infinite `λ` (strong observability)  if `obs = true`, in which case
the left Kronecker structure is preserved.

If conditions (1)-(4) are satisfied, the linearization is called `strongly minimal`, the resulting order `n`
is the least achievable order and `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The pencil based linearization is built using the methods described in [1] in conjunction with
pencil manipulation algorithms [2] and [3] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances for the
nonzero coefficients of `NUM(λ)` and `DEN(λ)`, respectively.

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[2] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[3] F.M. Dopico, M.C. Quintana and P. Van Dooren, Linear system matrices of rational transfer functions,
in "Realization and Model Reduction of Dynamical Systems, A Festschrift to honor the 70th birthday of Thanos Antoulas",
Eds. C. Beattie, P. Benner, M. Embree, S. Gugercin and S. Lefteriu, Springer-Verlag, 2020.
[arXiv:1903.05016](https://arxiv.org/pdf/1903.05016.pdf)
"""
function rm2lps(NUM::AbstractArray{T1,3},DEN::AbstractArray{T2,3}; minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
                fast::Bool = true, atol::Real = zero(real(T1)),
                rtol::Real = (min(size(NUM)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
    T = typeof(one(T1)/one(T2))
    Ap, Bp, Cp, D = rm2lspm(NUM,DEN)
    minimal && (contr = true; obs = true; noseig = true;)
    if ndims(D) == 2
        A, B, C, _  = lsminreal(Ap, Bp, Cp, contr = contr, obs = obs, fast = fast, atol = atol, rtol = rtol)
        F = zeros(T,size(B)...)
        G = zeros(T,size(C)...)
        H = zeros(T,size(D)...)
        E = Matrix{T}(I,size(A)...)
        nf = size(A,1)
        ni = 0
    else
        Ai, Ei, Bi, Fi, Ci, Gi, D, H  = lpsminreal(pm2lps(D)...; contr = contr, obs = obs,
                                                    fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
        Ap, Bp, Cp, _  = lsminreal(Ap, Bp, Cp, contr = contr, obs = obs, fast = fast, atol = atol,  rtol = rtol)
        nf = size(Ap,1)
        ni = size(Ai,1)
        A = [Ap zeros(T,nf,ni); zeros(T,ni,nf) Ai]
        E = [I zeros(T,nf,ni); zeros(T,ni,nf) Ei]
        B = [Bp; Bi]
        F = [zeros(T,size(Bp)...); Fi]
        C = [Cp Ci]
        G = [zeros(T,size(Cp)...) Gi]
    end
    return A, E, B, F, C, G, D, H, [nf, ni]
end
rm2lps(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
       DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       rm2lps(poly2pm(NUM),poly2pm(DEN); kwargs...)
function rm2lps(NUM::AbstractArray{T,3}; kwargs...) where T
    sys = pm2lps(NUM; kwargs...)
    return sys..., [0;size(sys[1],1)]
end
function rm2lps(NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    sys = pm2lps(poly2pm(NUM); kwargs...)
    return sys..., [0;size(sys[1],1)]
end
"""
    lps2rm(A, E, B, F, C, G, D, H; fast = true, atol1 = 0, atol2 = 0, gaintol = 0, rtol = min(atol1,atol2) > 0 ? 0 : n*ϵ, val) -> (NUM, DEN)

Build the rational matrix `R(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH := NUM(λ) ./ DEN(λ)` corresponding to its structured linearization

            | A-λE | B-λF |
     S(λ) = |------|------|
            | C-λG | D-λH |

by explicitly determining for each rational entry of `R(λ)`, its numerator coefficients from its finite zeros and gain, and
its denominator coefficients from its finite poles. An extension of the approach of [1] is employed, relying on the procedures
of [2] to compute strongly irreducible linearizations.

The keyword arguments `atol1` and `atol2` specify the absolute tolerances for the elements of `A`, `B`, `C`, `D`, and of
`E`, `F`, `G`, `H`, respectively,  and `rtol` specifies the relative tolerances for the nonzero elements of
`A`, `B`, `C`, `D`, `E`, F`, `G`, `H`.
The default relative tolerance is `(n+2)*ϵ`, where `n` is the size of the size dimension of `A`, and `ϵ` is the
machine epsilon of the element type of coefficients of `A`.

The keyword argument `gaintol` specifies the threshold for the magnitude of the nonzero elements of the gain matrix
`C*inv(γE-A)*B+D`, where `γ = val` if `val` is a number or `γ` is a randomly chosen complex value of unit magnitude,
if `val = missing`. Generally, `val` should not be a root of any of entries of `P`.

[1] A. Varga Computation of transfer function matrices of generalized state-space models.
    Int. J. Control, 50:2543–2561, 1989.

[2] F.M. Dopico, M.C. Quintana and P. Van Dooren, Linear system matrices of rational transfer functions,
in "Realization and Model Reduction of Dynamical Systems, A Festschrift to honor the 70th birthday of Thanos Antoulas",
Eds. C. Beattie, P. Benner, M. Embree, S. Gugercin and S. Lefteriu, Springer-Verlag, 2020.
[arXiv:1903.05016](https://arxiv.org/pdf/1903.05016.pdf)
"""
function lps2rm(A::AbstractMatrix, E::AbstractMatrix, B::AbstractMatrix, F::AbstractMatrix,
                C::AbstractMatrix, G::AbstractMatrix, D::AbstractMatrix, H::AbstractMatrix;
                fast::Bool = true, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                gaintol::Real = zero(real(eltype(A))), val::Union{Number,Missing} = missing,
                rtol::Real = ((size(A,1)+2)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)))
   n = LinearAlgebra.checksquare(A)
   (n,n) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   p, m = size(D)
   (n,m) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions"))
   (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions"))
   (n,m) == size(F) || throw(DimensionMismatch("B and F must have the same dimensions"))
   (p,n) == size(G) || throw(DimensionMismatch("C and G must have the same dimensions"))
   (p,m) == size(H) || throw(DimensionMismatch("D and H must have the same dimensions"))
   T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D), eltype(E), eltype(F), eltype(G), eltype(H))
   T <: BlasFloat || (T = promote_type(Float64,T))

   compl = T <: Complex
   ismissing(val) && (val = exp(rand()*im))
   Rval = lpseval(A, E, B, F, C, G, D, H, val)
   unimodular = isunimodular(A, E, atol1 = atol1, atol2 = atol2, rtol = rtol)
   NUM = zeros(T,p,m,n+2)
   if unimodular
      DEN = ones(T,p,m,1)
   else
      DEN = zeros(T,p,m,n+2)
      DEN[:,:,1] = ones(T,p,m)
   end
   for i = 1:p
       Ao, Eo, Bo, Fo, Co, Go, Do, Ho, Vo, Wo = lpsminreal(A, E, B, F, C[i:i,:], G[i:i,:], D[i:i,:], H[i:i,:],
                                                           contr = false, atol1 = atol1, atol2 = atol2, rtol = rtol)
       for j = 1:m
           A1, E1, B1, F1, C1, G1, D1, H1, _ = lpsminreal(Ao, Eo, (Bo/Wo)[:,j:j], (Fo/Wo)[:,j:j], Co, Go, (Do/Wo)[:,j:j], (Ho/Wo)[:,j:j],
                                           obs = false, atol1 = atol1, atol2 = atol2, rtol = rtol)
           if abs(Rval[i,j]) > gaintol
              zer, iz, _ = pzeros([A1 B1; C1 D1],[E1 F1; G1 H1]; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
              cz, zval = polcoeffval(zer[1:(length(zer)-sum(iz))],val)
              if unimodular
                 NUM[i,j,1:length(cz)] = compl ? cz*(Rval[i,j]/zval) : real(cz*(Rval[i,j]/zval))
              else
                 n1 = size(A1,1)
                 M1 = [A1 zeros(T,n1,2); zeros(T,1,n1+1) one(T); zeros(T,1,n1) one(T)  zero(T)]
                 N1 = [E1 F1 zeros(T,n1,1); G1 H1 zero(T); zeros(T,1,n1+2)]
                 pol, ip, _ = pzeros(M1,N1; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
                 cp, pval = polcoeffval(pol[1:(length(pol)-sum(ip))],val)
                 NUM[i,j,1:length(cz)] = compl ? cz*(Rval[i,j]*pval/zval) : real(cz*(Rval[i,j]*pval/zval))
                 DEN[i,j,1:length(cp)] = compl ? cp : real(cp)
              end
           end
       end
   end
   return NUM[:,:,1:pmdeg(NUM)+1], DEN[:,:,1:pmdeg(DEN)+1]
end
"""
     lpmfd2ls(DEN, NUM; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, C, D)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B |
     S(λ) = |------|---|
            |  C   | D |

of the left polynomial matrix fractional description `R(λ) = inv(DEN(λ))*NUM(λ)`,
such that `R(λ) = C*inv(λE-A)*B+D`. The resulting linearization `S(λ)` preserves a part,
if `minimal = false`, or the complete Kronecker structure, if `minimal = true`, of `R(λ)`.
In the latter case, the order `n` of `A-λE` is the least possible one.

`DEN(λ)` and `NUM(λ)` can be specified as polynomial matrices of the form `X(λ) = X_1 + λ X_2 + ... + λ**k X_(k+1)`,
for `X = DEN` or `X = NUM`, for which the coefficient matrices `X_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrices `X`, where `X[:,:,i]` contains the `i`-th coefficient matrix `X_i` (multiplying `λ**(i-1)`).

`DEN(λ)` and `NUM(λ)` can also be specified as matrices, vectors or scalars of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `obs = true`, in which case
the left Kronecker structure is preserved;

(4) `S(λ)` has no simple infinite eigenvalues if `minimal = true`.

If conditions (1)-(4) are satisfied, the linearization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the linearization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
For an irreducible linearization `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrices `DEN(λ)` and `NUM(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `DEN(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The structured linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function lpmfd2ls(DEN::Union{AbstractArray{T1,3},AbstractArray{T1,2}},NUM::Union{AbstractArray{T2,3},AbstractArray{T2,2}};
                  contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                  fast::Bool = true, atol::Real = zero(real(T1)),
                  rtol::Real = size(DEN,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1, T2}
    p = size(NUM,1)
    m = size(NUM,2)
    p1 =  size(DEN,1)
    p1 == size(DEN,2) || throw(DimensionMismatch("DEN(λ) must be a square polynomial matrix"))
    p == p1 || throw(DimensionMismatch("DEN(λ) and NUM(λ) must have the same number of rows"))

    return spm2ls(DEN, NUM, Matrix{T1}(I,p,p), zeros(T1,p,m), contr = contr, obs = obs, minimal = minimal,
                  fast = fast, atol = atol, rtol = rtol)
end
function lpmfd2ls(DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                  NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    return lpmfd2ls(poly2pm(DEN),poly2pm(NUM); kwargs...)
end
"""
     rpmfd2ls(DEN, NUM; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, C, D)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B |
     S(λ) = |------|---|
            |  C   | D |

of the right polynomial matrix fractional description `R(λ) = NUM(λ)*inv(DEN(λ))`,
such that `R(λ) = C*inv(λE-A)*B+D`. The resulting linearization `S(λ)` preserves a part,
if `minimal = false`, or the complete Kronecker structure, if `minimal = true`, of `R(λ)`.
In the latter case, the order `n` of `A-λE` is the least possible one.

`DEN(λ)` and `NUM(λ)` can be specified as polynomial matrices of the form `X(λ) = X_1 + λ X_2 + ... + λ**k X_(k+1)`,
for `X = DEN` or `X = NUM`, for which the coefficient matrices `X_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrices `X`, where `X[:,:,i]` contains the `i`-th coefficient matrix `X_i` (multiplying `λ**(i-1)`).

`DEN(λ)` and `NUM(λ)` can also be specified as matrices, vectors or scalars of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `obs = true`, in which case
the left Kronecker structure is preserved;

(4) `S(λ)` has no simple infinite eigenvalues if `minimal = true`.

If conditions (1)-(4) are satisfied, the linearization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the linearization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
For an irreducible linearization `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrices `DEN(λ)` and `NUM(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `DEN(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The structured linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function rpmfd2ls(DEN::Union{AbstractArray{T1,3},AbstractArray{T1,2}},NUM::Union{AbstractArray{T2,3},AbstractArray{T2,2}};
                  contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                  fast::Bool = true, atol::Real = zero(real(T1)),
                  rtol::Real = size(DEN,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1, T2}
    p = size(NUM,1)
    m = size(NUM,2)
    m1 =  size(DEN,1)
    m1 == size(DEN,2) || throw(DimensionMismatch("DEN(λ) must be a square polynomial matrix"))
    m == m1 || throw(DimensionMismatch("DEN(λ) and NUM(λ) must have the same number of rows"))

    return spm2ls(DEN, Matrix{T1}(I,m,m), NUM, zeros(T1,p,m), contr = contr, obs = obs, minimal = minimal,
                  fast = fast, atol = atol, rtol = rtol)
end
function rpmfd2ls(DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                  NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    return rpmfd2ls(poly2pm(DEN),poly2pm(NUM); kwargs...)
end
"""
     lpmfd2lps(DEN, NUM; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, F, C, G, D, H)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B-λF |
     S(λ) = |------|------|
            | C-λG | D-λH |

of the left polynomial matrix fractional description `R(λ) = inv(DEN(λ))*NUM(λ)`,
such that `R(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH`. The resulting linearization `S(λ)` preserves a part,
if `minimal = false`, or the complete Kronecker structure, if `minimal = true`, of `R(λ)`.
In the latter case, the order `n` of `A-λE` is the least possible one.

`DEN(λ)` and `NUM(λ)` can be specified as polynomial matrices of the form `X(λ) = X_1 + λ X_2 + ... + λ**k X_(k+1)`,
for `X = DEN` or `X = NUM`, for which the coefficient matrices `X_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrices `X`, where `X[:,:,i]` contains the `i`-th coefficient matrix `X_i` (multiplying `λ**(i-1)`).

`DEN(λ)` and `NUM(λ)` can also be specified as matrices, vectors or scalars of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B-λF A-λE] = n` (strong controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C-λG] = n` (strong observability)  if `minimal = true` or `obs = true`, in which case
the left Kronecker structure is preserved.

If conditions (1)-(3) are satisfied, the linearization is called `strongly minimal`, the resulting order `n`
is the least achievable order and `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrices `DEN(λ)` and `NUM(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `DEN(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The structured linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function lpmfd2lps(DEN::Union{AbstractArray{T1,3},AbstractArray{T1,2}},NUM::Union{AbstractArray{T2,3},AbstractArray{T2,2}};
                  contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                  fast::Bool = true, atol::Real = zero(real(T1)),
                  rtol::Real = size(DEN,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1, T2}
    p = size(NUM,1)
    m = size(NUM,2)
    p1 =  size(DEN,1)
    p1 == size(DEN,2) || throw(DimensionMismatch("DEN(λ) must be a square polynomial matrix"))
    p == p1 || throw(DimensionMismatch("DEN(λ) and NUM(λ) must have the same number of rows"))

    return spm2lps(DEN, NUM, Matrix{T1}(I,p,p), zeros(T1,p,m), contr = contr, obs = obs, minimal = minimal,
                  fast = fast, atol = atol, rtol = rtol)
end
function lpmfd2lps(DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                   NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    return lpmfd2lps(poly2pm(DEN),poly2pm(NUM); kwargs...)
end
"""
     rpmfd2lps(DEN, NUM; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, F, C, G, D, H)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B-λF |
     S(λ) = |------|------|
            | C-λG | D-λH |

of the right polynomial matrix fractional description `R(λ) = NUM(λ)*inv(DEN(λ))`,
such that `R(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH`. The resulting linearization `S(λ)` preserves a part,
if `minimal = false`, or the complete Kronecker structure, if `minimal = true`, of `R(λ)`.
In the latter case, the order `n` of `A-λE` is the least possible one.

`DEN(λ)` and `NUM(λ)` can be specified as polynomial matrices of the form `X(λ) = X_1 + λ X_2 + ... + λ**k X_(k+1)`,
for `X = DEN` or `X = NUM`, for which the coefficient matrices `X_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrices `X`, where `X[:,:,i]` contains the `i`-th coefficient matrix `X_i` (multiplying `λ**(i-1)`).

`DEN(λ)` and `NUM(λ)` can also be specified as matrices, vectors or scalars of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B-λF A-λE] = n` (strong controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C-λG] = n` (strong observability)  if `minimal = true` or `obs = true`, in which case
the left Kronecker structure is preserved.

If conditions (1)-(3) are satisfied, the linearization is called `strongly minimal`, the resulting order `n`
is the least achievable order and `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrices `DEN(λ)` and `NUM(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `DEN(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The structured linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function rpmfd2lps(DEN::Union{AbstractArray{T1,3},AbstractArray{T1,2}},NUM::Union{AbstractArray{T2,3},AbstractArray{T2,2}};
                  contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                  fast::Bool = true, atol::Real = zero(real(T1)),
                  rtol::Real = size(DEN,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1, T2}
    p = size(NUM,1)
    m = size(NUM,2)
    m1 =  size(DEN,1)
    m1 == size(DEN,2) || throw(DimensionMismatch("DEN(λ) must be a square polynomial matrix"))
    m == m1 || throw(DimensionMismatch("DEN(λ) and NUM(λ) must have the same number of rows"))

    return spm2lps(DEN, Matrix{T1}(I,m,m), NUM, zeros(T1,p,m), contr = contr, obs = obs, minimal = minimal,
                  fast = fast, atol = atol, rtol = rtol)
end
function rpmfd2lps(DEN::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                   NUM::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    return rpmfd2lps(poly2pm(DEN),poly2pm(NUM); kwargs...)
end
"""
     pminv2ls(P; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, C, D)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B |
     S(λ) = |------|---|
            |  C   | D |

of the inverse `R(λ) = inv(P(λ))`, such that `R(λ) = C*inv(λE-A)*B+D`.
The resulting linearization `S(λ)` preserves a part, if `minimal = false`, or
the complete Kronecker structure, if `minimal = true`, of `R(λ)`.
In the latter case, the order `n` of `A-λE` is the least possible one.

`P(λ)` can be specified as a polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrix `P`, where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `obs = true`, in which case
the left Kronecker structure is preserved;

(4) `S(λ)` has no simple infinite eigenvalues if `minimal = true`.

If conditions (1)-(4) are satisfied, the linearization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the linearization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
For an irreducible linearization `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrix `P(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `P(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The structured linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function pminv2ls(P::Union{AbstractArray{T1,3},AbstractArray{T1,2}};
                  contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                  fast::Bool = true, atol::Real = zero(real(T1)),
                  rtol::Real = size(P,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1}
    n = size(P,1)
    n == size(P,2) || throw(DimensionMismatch("P(λ) must be a square polynomial matrix"))

    return spm2ls(P, Matrix{T1}(I,n,n), Matrix{T1}(I,n,n), zeros(T1,n,n), contr = contr, obs = obs, minimal = minimal,
                  fast = fast, atol = atol, rtol = rtol)
end
function pminv2ls(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    return pminv2ls(poly2pm(P); kwargs...)
end
"""
     pminv2lps(P; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, F, C, G, D, H)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B-λF |
     S(λ) = |------|------|
            | C-λG | D-λH |

of the inverse `R(λ) = inv(P(λ))`, such that `R(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH`.
The resulting linearization `S(λ)` preserves a part, if `minimal = false`, or
the complete Kronecker structure, if `minimal = true`, of `R(λ)`.
In the latter case, the order `n` of `A-λE` is the least possible one.

`P(λ)` can be specified as a polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrix `P`, where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B-λF A-λE] = n` (strong controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C-λG] = n` (strong observability)  if `minimal = true` or `obs = true`, in which case
the left Kronecker structure is preserved.

If conditions (1)-(3) are satisfied, the linearization is called `strongly minimal`, the resulting order `n`
is the least achievable order and `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `R(λ)`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrix `P(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `P(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The structured linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function pminv2lps(P::Union{AbstractArray{T1,3},AbstractArray{T1,2}};
                  contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                  fast::Bool = true, atol::Real = zero(real(T1)),
                  rtol::Real = size(P,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1}
    n = size(P,1)
    n == size(P,2) || throw(DimensionMismatch("P(λ) must be a square polynomial matrix"))

    return spm2lps(P, Matrix{T1}(I,n,n), Matrix{T1}(I,n,n), zeros(T1,n,n), contr = contr, obs = obs, minimal = minimal,
                  fast = fast, atol = atol, rtol = rtol)
end
function pminv2lps(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    return pminv2lps(poly2pm(P); kwargs...)
end


text:
# -------------------------------
#  Cyclic data exchange example - Master
# -------------------------------

# Every worker will have a global variable "mynumber" which holds some data (a random number in this example).
# The goal is to make the workers exchange this data cyclically (forward it to the next worker).
# We will realize this exchange by multiple channels, one on each worker.
# Every worker will push his number to the next workers channel and check his own channel to receive a new number.

using Distributed

N = 4

# Create workers
addprocs(N)

# Create a channel on every worker for receiving data.
# Store RemoteChannel handles to those channels in a dictoniary pid => RemoteChannel
const channels = Dict(p=>RemoteChannel(() -> Channel{Int}(1), p) for p in workers())

# Make the dictionary a global constant on every worker process
@everywhere workers() channels=$channels # will also be a constant global on the worker processes

# Function for displaying the current data ("mynumber") distribution.
function print_numbers()
    for p in workers()
        println("$p => ", fetch(@spawnat p mynumber))
    end
end

# Will perform one cyclic exchange between the workers.
function makestep()
    # run worker script on every worker
    @everywhere workers() include("worker.jl")
    print_numbers()
end

# # Tell every worker who his neighbor is (create global constant "npid")
# # We could use ParallelDataTransfer.jl here to make this nicer.
# w = workers()
# nw = nworkers()
# @sync for i in 1:length(w)
#     pid = w[i]
#     npid = w[mod1(i+1, nw)]
#     @spawnat pid eval(:(global const npid = $npid))
# end

max_stars_repo_name: crstnbr/julia-workshop


text:
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

sol = 5.0
#Δt, Nt = 0.08/sol,400
Δt, Nt = 0.120000123/sol,400

#D, Δx = 1.0, 0.25
D, Δx = 1.0, 0.1
#Γ = meshsphere(D, Δx)
Γ = meshrectangle(D,D,Δx,3)
#γ = boundary(Γ)
γ1 = meshsegment(D,D,3)
γ2 = CompScienceMeshes.translate(γ1, point(0,D,0))

#X = raviartthomas(Γ,weld(γ1,γ2))
X = rt_ports(Γ,(γ1,γ2))
X.fns = X.fns[1:281]
numfunctions(X)

T = timebasisshiftedlagrange(Δt, Nt, 3)
U = timebasisdelta(Δt, Nt)

V = X ⊗ T
W = X ⊗ U

duration, delay, amplitude = 8.0/sol, 12.0/sol, 1.0
gaussian = creategaussian(duration, delay, duration)

direction, polarisation = z, x
E = planewave(polarisation, direction, derive(gaussian), sol)
T = MWSingleLayerTDIO(sol,-1/sol,-sol,2,0)

#Z = assemble(T,W,V);

@hilbertspace j
@hilbertspace j′
tdefie = @discretise T[j′,j] == -1E[j′]   j∈V  j′∈W
xefie = solve(tdefie)

using Plotly
# include(Pkg.dir("CompScienceMeshes","examples","matlab_patches.jl"))

Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
ω = collect(ω0 + (0:Nt-1)*Δω)
_, i1 = findmin(abs.(ω.-sol))

ω1 = ω[i1]
ue = Xefie[:,i1] / fouriertransform(gaussian)(ω1)

using LinearAlgebra
fcre, geo = facecurrents(ue, X)
t2 = patch(geo, real.(norm.(fcre)))
Plotly.plot(t2)

fcr_td, geo = facecurrents(xefie[:,end], X)
# fcr_ch, geo = facecurrents(xefie[:,end], divergence(X))

t3 = patch(Γ,real.(norm.(fcr_td)))
Plotly.plot(t3)
# jmatlab_quiver(Γ, fcr_td)


text:
include("Xpress.jl")

using Xpress

println(Xpress.version)

using Base.Test
using MathProgBase
using MathProgBase.SolverInterface

function linprogsolvertest(solver::AbstractMathProgSolver)

    m = LinearQuadraticModel(solver)

    # Min -x
    # s.t. x + y <= 1
    # x, y >= 0
    loadproblem!(m, [1. 1.], [0.,0.], [Inf, Inf], [-1.,0.], [-Inf],[1.], :Min)
    @test numvar(m) == 2
    @test numconstr(m) == 1
    @test getvarLB(m) == [0.,0.]
    @test getvarUB(m)[1] >= 1e20
    @test getvarUB(m)[2] >= 1e20
    @test getconstrLB(m)[1] < -1e30
    @test getconstrUB(m)[1] == 1.0
    @test getobj(m) == [-1.,0.]
    @test getsense(m) == :Min

    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) -1
    @test_approx_eq getsolution(m) [1.0, 0.0]
    @test_approx_eq getconstrsolution(m) [1.0]
    @test_approx_eq getconstrduals(m) [-1.0]
    @test_approx_eq getreducedcosts(m) [0.0, 1.0]

    setsense!(m, :Max)
    # max x
    setobj!(m, [1.0,0.0])
    updatemodel!(m)
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 1
    @test_approx_eq getsolution(m) [1.0, 0.0]
    @test_approx_eq getconstrsolution(m) [1.0]
    @test_approx_eq getconstrduals(m) [1.0]
    @test_approx_eq getreducedcosts(m) [0.0, -1.0]

    # add new variable to get:
    # max x + 2z
    # s.t. x + y + z <= 1
    # x,y,z >= 0
    addvar!(m, [1], [1.0], 0, Inf, 2.0)
    updatemodel!(m)

    @test numvar(m) == 3
    @test numconstr(m) == 1

    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 2
    @test_approx_eq getsolution(m) [0.0, 0.0, 1.0]
    @test_approx_eq getconstrsolution(m) [1.0]
    @test_approx_eq getconstrduals(m) [2.0]
    @test_approx_eq getreducedcosts(m) [-1.0, -2.0, 0.0]

    setvarLB!(m, [-1.0,0.0,0.0])
    updatemodel!(m)
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 3

    # fix z to zero
    setvarLB!(m, [0.0,0.0,0.0])
    setvarUB!(m, [Inf,Inf,0.0])
    updatemodel!(m)
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 1

    setconstrUB!(m, [2.0])
    setconstrLB!(m, [2.0])
    updatemodel!(m)
    optimize!(m)
    @test_approx_eq getobjval(m) 2

    setobj!(m, [1.0, 2.0, 0.0])
    updatemodel!(m)
    optimize!(m)
    @test_approx_eq getobjval(m) 4
    @test_approx_eq getsolution(m) [0.0, 2.0, 0.0]

    # we now have:
    # max x+2y
    # s.t. x + y + z == 2
    # x,y >= 0, z = 0

    # add constraint x - y >= 0
    addconstr!(m, [1,2], [1.0,-1.0], 0.0, Inf)
    updatemodel!(m)
    optimize!(m)

    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) 3
    @test_approx_eq getsolution(m) [1.0, 1.0, 0.0]
    @test_approx_eq getconstrduals(m) [1.5, -0.5]
    @test_approx_eq getreducedcosts(m) [0.0, 0.0, -1.5]

    # test addvar! interface

    m = LinearQuadraticModel(solver)

    # Min -x
    # s.t. x + y <= 1
    # x, y >= 0
    addvar!(m, 0, Inf, -1)
    addvar!(m, 0, Inf, 0)
    updatemodel!(m)
    addconstr!(m, [1, 2], [1.0, 1.0], -Inf, 1.0)
    updatemodel!(m)
    @test numvar(m) == 2
    @test numconstr(m) == 1
    @test getvarLB(m) == [0.,0.]
    @test getvarUB(m)[1] >= 1e20
    @test getvarUB(m)[2] >= 1e20
    @test getconstrLB(m)[1] <= -1e20
    @test getconstrUB(m)[1] == 1.0
    @test getobj(m) == [-1.,0.]
    @test getsense(m) == :Min

    optimize!(m)
    @test status(m) == :Optimal
    @test_approx_eq getobjval(m) -1
    @test_approx_eq getsolution(m) [1.0, 0.0]
    @test_approx_eq getconstrsolution(m) [1.0]
    @test_approx_eq getconstrduals(m) [-1.0]
    @test_approx_eq getreducedcosts(m) [0.0, 1.0]

end

linprogsolvertest( XpressSolver() )


text:
function addPackedFactor(navAbilityClient::NavAbilityClient, client::Client, factor)::String
    response = navAbilityClient.mutate(MutationOptions(
        "addFactor",
        MUTATION_ADDFACTOR,
        Dict(
            "factor" => Dict(
                "client" => client,
                "packedData" => json(factor)
            )
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return "Error" end
    addFactor = get(data,"addFactor","Error")
    return addFactor
end

function addFactor(navAbilityClient::NavAbilityClient, client::Client, factor::Factor)
    return @async addPackedFactor(navAbilityClient, client, factor)
end

function _getFactorEvent(navAbilityClient::NavAbilityClient, client::Client, label::String)::Dict{String,Any}
    response = navAbilityClient.query(QueryOptions(
        "sdk_get_factor",
        """
            $GQL_FRAGMENT_FACTORS
            $GQL_GETFACTOR
        """,
        Dict(
            "label" => label,
            "userId" => client.userId,
            "robotId" => client.robotId,
            "sessionId" => client.sessionId
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return Dict() end
    user = get(data,"users",[])
    if size(user)[1] < 1 return Dict() end
    robots = get(user[1],"robots",[])
    if size(robots)[1] < 1 return Dict() end
    sessions = get(robots[1],"sessions",[])
    if size(sessions)[1] < 1 return Dict() end
    factors = get(sessions[1],"factors",[])
    if size(factors)[1] < 1 return Dict() end
    return factors[1]
end

getFactor(navAbilityClient::NavAbilityClient, client::Client, label::String) = @async _getFactorEvent(navAbilityClient, client, label)

function getFactorsEvent(navAbilityClient::NavAbilityClient, client::Client; detail::QueryDetail = SKELETON)::Vector{Dict{String,Any}}
    response = navAbilityClient.query(QueryOptions(
        "sdk_get_factors",
        """
            $GQL_FRAGMENT_FACTORS
            $GQL_GETFACTORS
        """,
        Dict(
            "userId" => client.userId,
            "robotId" => client.robotId,
            "sessionId" => client.sessionId,
            "fields_summary" => detail === SUMMARY || detail === FULL,
            "fields_full" => detail === FULL,
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return [] end
    user = get(data,"users",[])
    if size(user)[1] < 1 return [] end
    robots = get(user[1],"robots",[])
    if size(robots)[1] < 1 return [] end
    sessions = get(robots[1],"sessions",[])
    if size(sessions)[1] < 1 return [] end
    return get(sessions[1],"factors",[])
end

getFactors( navAbilityClient::NavAbilityClient, client::Client; detail::QueryDetail = SKELETON) = @async getFactorsEvent(navAbilityClient, client; detail )

function listFactors(navAbilityClient::NavAbilityClient, client::Client)
    @async begin
        factors = getFactors(navAbilityClient,client) |> fetch
        map(v -> v["label"], factors)
    end
end

function lsf(navAbilityClient::NavAbilityClient, client::Client)
    return listFactors(navAbilityClient,client)
end


text:
using Distributions, Parameters
import Base.rand
import Distributions: logpdf

struct LNR{T1,T2,T3} <: ContinuousUnivariateDistribution
    μ::T1
    σ::T2
    ϕ::T3
end

Broadcast.broadcastable(x::LNR) = Ref(x)

LNR(;μ, σ, ϕ) = LNR(μ, σ, ϕ)

function rand(dist::LNR)
    @unpack μ,σ,ϕ = dist
    x = @. rand(LogNormal(μ, σ)) + ϕ
    rt,resp = findmin(x)
    return resp,rt
end

rand(dist::LNR, N::Int) = [rand(dist) for i in 1:N]

function logpdf(d::T, r::Int, t::Float64) where {T<:LNR}
    @unpack μ,σ,ϕ = d
    LL = 0.0
    for (i,m) in enumerate(μ)
        if i == r
            LL += logpdf(LogNormal(m, σ), t-ϕ)
        else
            LL += log(1-cdf(LogNormal(m, σ), t-ϕ))
        end
    end
    return LL
end

logpdf(d::LNR, data::Tuple) = logpdf(d, data...)

@model lnrModel(data, idx, N, Nsub, ::Type{T}=Vector{Float64},
    ::Type{T1}=Array{Float64,2}) where {T,T1}  = begin
    Nr = 3
    mu_g = T(undef,Nr)
    sigma_g ~ Truncated(Cauchy(0, 1), 0, Inf)
    mu_g ~ [Normal(0,3)]
    mu = T1(undef,Nsub,Nr)
    for s in 1:Nsub
        for i in 1:Nr
            mu[s,i] ~ Normal(mu_g[i], sigma_g)
        end
    end
    sigma ~ Truncated(Cauchy(0, 1), 0, Inf)
    for i in 1:N
        data[i] ~ LNR(mu[idx[i],:], sigma, 0.0)
    end
end

function simulateLNR(; Nd=1, Nsub, mu_g = [-1.0,-2.0,-.5], sigma_g = 1,
     sigma = 1, kwargs...)
    N = Nd * Nsub
    choice = fill(0, N)
    rt = fill(0.0, N)
    idx = similar(choice)
    i = 0
    for s in 1:Nsub
        mu = @. rand(Normal(mu_g, sigma_g))
        dist = LNR(μ=mu, σ=sigma, ϕ=0.0)
        for _ in 1:Nd
            i += 1
            idx[i] = s
            choice[i],rt[i] = rand(dist)
        end
    end
    return (choice=choice, rt=rt, idx=idx, N=N, Nsub=Nsub)
 end

function AHMClnr(choice, rt, idx, N, Nsub)
    data=[(c,r) for (c,r) in zip(choice,rt)]

    return lnrModel(data, idx, N, Nsub)
end

AHMCconfig = Turing.NUTS(1000,.85)

CmdStanLNR = "
functions{

    real LNR_LL(row_vector mus,real sigma,real ter,real v,int c){
      real LL;
      LL = 0;
      for(i in 1:num_elements(mus)){
        if(i == c){
          LL += lognormal_lpdf(v-ter|mus[i],sigma);
      }else{
        LL += log(1-lognormal_cdf(v-ter,mus[i],sigma));
      }
    }
      return LL;
    }

    }

data{
     int N;
     int Nr;
     int Nsub;
     vector[N] rt;
     vector[N] choice;
}

parameters {
     vector[Nr] mu_g;
     matrix[Nsub,Nr] mu;
     real<lower=0> sigma_g;
     real<lower=0> sigma;
}

model {
     mu_g ~ normal(0, 3);
     sigma_g ~ cauchy(0,1);
     for(s in 1:Nsub){
         for(i in 1:Nr){
            mu[s,i] ~ normal(mu_g[i], sigma_g);
        }
     }
     sigma ~ cauchy(0,1);
     for(i in 1:N){
         target += LNR_LL(mu[idx[i],:], sigma, 0.0, rt[i], choice[i]);
      }
}
"

CmdStanConfig = Stanmodel(name="CmdStanLNR", model=CmdStanLNR, nchains=1,
   Sample(num_samples=1000, num_warmup=1000, adapt=CmdStan.Adapt(delta=0.8),
   save_warmup=true))

   struct LNRProb{T}
      data::T
      N::Int
      Nc::Int
  end

  function (problem::LNRProb)(θ)
      @unpack data=problem
      @unpack v,A,k,tau=θ
      d=LBA(ν=v, A=A, k=k, τ=tau)
      minRT = minimum(x->x[2], data)
      logpdf(d,data)+sum(logpdf.(TruncatedNormal(0, 3, 0, Inf), v)) +
      logpdf(TruncatedNormal(.8, .4, 0, Inf),A)+logpdf(TruncatedNormal(.2, .3 ,0 ,Inf), k)+
      logpdf(TruncatedNormal(.4, .1, 0, minRT), tau)
  end

function sampleDHMC(choice, rt, N, Nc, nsamples)
    data = [(c, r) for (c, r) in zip(choice, rt)]
    return sampleDHMC(data, N, Nc, nsamples)
end

# Define problem with data and inits.
function sampleDHMC(data, N, Nc, nsamples, autodiff)
    p = LBAProb(data, N, Nc)
    p((v=fill(.5, Nc), A=.8, k=.2, tau=.4))
    # Write a function to return properly dimensioned transformation.
    trans = as((v=as(Array, asℝ₊, Nc),A=asℝ₊, k=asℝ₊, tau=asℝ₊))
    # Use Flux for the gradient.
    P = TransformedLogDensity(trans, p)
    ∇P = ADgradient(autodiff, P)
    # FSample from the posterior.
    n = dimension(trans)
    results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, nsamples;
        q = zeros(n), p = ones(n),reporter = NoProgressReport())
    # Undo the transformation to obtain the posterior from the chain.
    posterior = transform.(trans, results.chain)
    chns = nptochain(results, posterior)
    return chns
end

function simulateLBA(;Nd, v=[1.0,1.5,2.0], A=.8, k=.2, tau=.4, kwargs...)
    return (rand(LBA(ν=v, A=A, k=k, τ=tau), Nd)..., N=Nd, Nc=length(v))
end