text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 9f18ccd0-8e25-11eb-20e6-ab941f663550
begin
	import Pkg
    Pkg.activate(".")
    Pkg.add("PlutoUI")
	Pkg.add("Plots")

	using PlutoUI
	using Plots
end

# ╔═╡ 10a9374c-8e1a-11eb-10b4-fb8272dd46d1
md"""
## SIR Model Simulation

#### Intro

The **SIR Model** has been developed to simulate an epidemic over time. The model consists of a system of 4 differential equations that express the rates of change of 4 variables over time. The 4 variables are:

1. **S** - the susceptibles of getting the infection
2. **I** - the infected
3. **R** - the recovered from the infection
4. **D** - the number of dead

#### The Model
Here follow the 4 equations that govern the model dynamics:

${ \begin{align*}
\frac{\mathrm{d}S}{\mathrm{d}t} & = -\beta \cdot \frac{S \cdot I}{N} + \delta \cdot R \\
\frac{\mathrm{d}I}{\mathrm{d}t} & = \beta \cdot \frac{S \cdot I}{N} + \gamma \cdot I \\
\frac{\mathrm{d}R}{\mathrm{d}t} & = \gamma \cdot (1-\rho) \cdot I - \delta \cdot R \\
\frac{\mathrm{d}D}{\mathrm{d}t} & = \gamma \cdot \rho \cdot I
\end{align*} }$

#### Model simulation
The following is a simulation of the model described above. First choose the model parameters
"""

# ╔═╡ 43970fa0-8e1b-11eb-3052-1701966a4478
begin
	T = 38*30 	#period of 300 daysplot
	Δt = 1/400 	#time interval of 6 hours (1/4 of a day)
end;

# ╔═╡ 742d1fd0-907d-11eb-27d1-59ca2af59f97
parameter = Dict(
	"rate of infection" => 0.14,
	"rate of recovery" => 0.07,
	"rate of immuity loss" => 0.00,
	"rate of death of infected" => 0.001,
	"total population" => 6*10^7,
	"initial number of infected" => 100.0
);

# ╔═╡ d59335ac-913b-11eb-228f-273be8d7a86a
md"""
---
**Additional parameter and functions for simulations with two subgroups**
"""

# ╔═╡ 56b9df4c-913b-11eb-180a-eb784c2b5181
begin
	parameter_subpopulation = Dict(
		"rate of infection" => 0.0,
		"rate of recovery" => 0.035,
		"rate of immuity loss" => 0.0,
		"rate of death of infected" => 0.2,
		"total population" => 1*10^6,
		"initial number of infected" => 0.0
		)

	mixing_parameter = Dict(
		"rate of infection 1 -> 2" => 0.0,
		"rate of infection 2 -> 1" => 0.0
	)
end;

# ╔═╡ 6618e836-9080-11eb-39f4-21ac8e8af61a
md"""
---
**Additional Parameter and functions for simulations with an incidence dependent recovery rate with delay.**
"""

# ╔═╡ f8e23e02-9142-11eb-1700-6974da1a09c0
md"""
---
#### Worksheet
Here all the calculations start. Set the initial state, choose the set of parameters and if needed the parameter changing function.
"""

# ╔═╡ da008842-9142-11eb-388c-d7e1f662d9db
md"""
---
#### Plots
"""

# ╔═╡ 93b1fcc4-9145-11eb-3d5e-4b9bd577a545
md"""
---
#### Model Functions

**Euler Algorithm to solve system. With constant or changing parameters**
"""

# ╔═╡ 2680b820-8e32-11eb-28de-6f039360524b
begin
	function euler(f::Function,change_par!::Function,t0,tn,Δt,x₀,par)
		T = t0:Δt:tn
		F = Matrix{eltype(x₀)}(undef,(length(T),length(x₀)))
		F[1,:] = x₀
		for (n,t) ∈ enumerate(T[2:end])
			change_par!(par,t,F,n)
			F[n+1,:] = F[n,:] + f(t,F,n,par) .* Δt
		end
		return F
	end

	function no_change!(par,t,F,n) end

	euler(f::Function,t0,tn,Δt,x₀,par) = euler(f,no_change!,t0,tn,Δt,x₀,par)
end;

# ╔═╡ 58f42866-9146-11eb-29ea-87830bc8db66
md"""
**Defines the infinitesimal change of system**
"""

# ╔═╡ 3ab68bec-8e2f-11eb-203e-f9530bc39073
begin
	function SIR_step(t,x,β,δ,γ,ρ,N)
		S, I, R, D = x[1], x[2], x[3], x[4]
		nI = β*I*S/N
		return [
			-nI + δ*R,
			nI - γ*I,
			γ*(1-ρ)*I - δ*R,
			γ*ρ*I,
			nI / Δt
			]
	end
	SIR_step(t,x,par::Dict) = SIR_step(t,x,
				par["rate of infection"],par["rate of immuity loss"],
				par["rate of recovery"],par["rate of death of infected"],
				par["total population"]
			)
	SIR_step(t,H,n,par::Dict) = SIR_step(t,H[n,:],par)
end;

# ╔═╡ a7ac774a-913b-11eb-267b-8f9841c44cb4
begin
	function SIR_step_2D(t,H,n,par₁::Dict,par₂::Dict,par_mix::Dict)
		x₁,x₂ = view(H,n,1:4), view(H,n,5:8)
		x = vcat(SIR_step(t,x₁,par₁),SIR_step(t,x₂,par₂))
		return x + SIR_interaction(t,x₁,x₂,
					par_mix["rate of infection 1 -> 2"],
					par_mix["rate of infection 2 -> 1"],
					par₁["total population"],par₂["total population"])
	end

	function SIR_interaction(t,x₁,x₂,β₃,β₄,N₁,N₂)
		S₁, I₁, S₂, I₂ = x₁[1], x₁[2], x₂[1], x₂[2]
		return [
			- β₄*S₁*I₂/N₁, β₄*S₁*I₂/N₁, 0, 0,
			- β₃*S₂*I₁/N₁, β₃*S₂*I₁/N₁, 0, 0
		]
	end
end;

# ╔═╡ de101aba-9145-11eb-2457-7b1831145eeb
md"""

**Supporting function for calculating the Incidence**
"""

# ╔═╡ 16a50518-906d-11eb-1e2b-3faa1765f30c
begin
	#Calculates at one time within the simulation
	function incidence(nI,t,pop_size;ndays=7,scale=100_000)
		#set starting time in range
		t₀ = max(t-round(Int,ndays/Δt),1)
		#return Incidence
		return (nI[t] - nI[t₀]) .* (scale/pop_size)
	end

	#Calculates for all times after simulation
	function incidence(nI,pop_size;ndays=7,scale=100_000)
		#rescale days
		ndays = round(Int,ndays/Δt)
		#sum new cases in first days
		I = nI[1:ndays]
		#iterate and add incidence
		for t₀ ∈ 1:length(nI)-ndays
			append!(I,nI[t₀+ndays] - nI[t₀])
		end
		#scale and return
		return I .* scale/pop_size
	end
end;

# ╔═╡ ce54bf76-906f-11eb-2480-75667d715789
begin
	#Additional parameter
	parameter_incidence_dependent_recovery_rate = merge(parameter,Dict(
			"high factor" => 2.8,
			"low incidence" => 17,
			"high incidence" => 200,
			"delay" => round(Int,5/Δt),
			"original recovery rate" => parameter["rate of recovery"]
			))

	#Dependency of rate of recovery on incidence
	function incidence_dependent_recovery_rate!(par,t,F,n)
		if n > par["delay"]
			i = incidence(view(F,:,5),n-par["delay"],par["total population"])
			if i > par["high incidence"]
				par["rate of recovery"] = par["original recovery rate"]*par["high factor"]
			end
			if i < par["low incidence"]
				par["rate of recovery"] = par["original recovery rate"]
			end
		end
	end
end;

# ╔═╡ c6bbb678-9080-11eb-2047-a3dcd34e8c77
begin
	N = parameter_incidence_dependent_recovery_rate["total population"]
	I₀ = parameter_incidence_dependent_recovery_rate["initial number of infected"]

	x₀ = [
		N-I₀,
		I₀,
		0.0,
		0.0,
		0.0
	]

	result = euler(
		SIR_step,
		incidence_dependent_recovery_rate!,
		0.0,T,Δt,
		x₀,
		parameter_incidence_dependent_recovery_rate
		)
end;

# ╔═╡ 3d9185a6-908c-11eb-301c-358ed4a3eb37
let
	p = plot(framestlye=:zerolines,size=(680,450))

	plot!(p, 0.0:Δt:T,result[:,1:4],
		label = ["S" "I" "R" "D"],
		color = [:blue :red :green :yellow],
		xlabel = "Days",
		ylabel = "Number of Individuals"
	)
end

# ╔═╡ b2e6d55c-9142-11eb-3008-89e56b2af8db
let
	p = plot(framestlye=:zerolines,size=(680,450))

	plot!(p, 0.0:Δt:T,result[:,5],
		label = "I",
		color = :blue ,
		xlabel = "Days",
		ylabel = "Number of Cases"
	)

	p
end

# ╔═╡ d492c47e-908c-11eb-36bd-45fbb77d545e
let
	p = plot(framestlye=:zerolines,size=(680,450))

	I = incidence(result[:,5],parameter["total population"])

	plot!(p, 0.0:Δt:T,I,
		label = "I",
		color = :blue ,
		xlabel = "Days",
		ylabel = "7 day incidence per 100 000"
	)

	p
end

# ╔═╡ ceedca8e-9148-11eb-189c-a3a088520c2c
md"""
---
Package Management
"""

# ╔═╡ Cell order:
# ╟─10a9374c-8e1a-11eb-10b4-fb8272dd46d1
# ╠═43970fa0-8e1b-11eb-3052-1701966a4478
# ╠═742d1fd0-907d-11eb-27d1-59ca2af59f97
# ╟─d59335ac-913b-11eb-228f-273be8d7a86a
# ╠═56b9df4c-913b-11eb-180a-eb784c2b5181
# ╠═a7ac774a-913b-11eb-267b-8f9841c44cb4
# ╟─6618e836-9080-11eb-39f4-21ac8e8af61a
# ╠═ce54bf76-906f-11eb-2480-75667d715789
# ╟─f8e23e02-9142-11eb-1700-6974da1a09c0
# ╠═c6bbb678-9080-11eb-2047-a3dcd34e8c77
# ╟─da008842-9142-11eb-388c-d7e1f662d9db
# ╟─3d9185a6-908c-11eb-301c-358ed4a3eb37
# ╟─d492c47e-908c-11eb-36bd-45fbb77d545e
# ╟─b2e6d55c-9142-11eb-3008-89e56b2af8db
# ╟─93b1fcc4-9145-11eb-3d5e-4b9bd577a545
# ╠═2680b820-8e32-11eb-28de-6f039360524b
# ╟─58f42866-9146-11eb-29ea-87830bc8db66
# ╠═3ab68bec-8e2f-11eb-203e-f9530bc39073
# ╟─de101aba-9145-11eb-2457-7b1831145eeb
# ╠═16a50518-906d-11eb-1e2b-3faa1765f30c
# ╟─ceedca8e-9148-11eb-189c-a3a088520c2c
# ╠═9f18ccd0-8e25-11eb-20e6-ab941f663550

max_stars_repo_name: roccminton/PlutoNotebooks

max_issues_repo_name: roccminton/PlutoNotebooks

max_forks_repo_name: roccminton/PlutoNotebooks


text:
using YaoPastaQ
using Documenter

DocMeta.setdocmeta!(YaoPastaQ, :DocTestSetup, :(using YaoPastaQ); recursive=true)

makedocs(;
    modules=[YaoPastaQ],
    authors="Roger-Luo <rogerluo.rl18@gmail.com> and contributors",
    repo="https://github.com/QuantumBFS/YaoPastaQ.jl/blob/{commit}{path}#{line}",
    sitename="YaoPastaQ.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://QuantumBFS.github.io/YaoPastaQ.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/QuantumBFS/YaoPastaQ.jl",
)

ext: jl


text:
# Adapted from:
# SOSDEMO7 --- Chebyshev polynomials
# Section 3.7 of SOSTOOLS User's Manual

using Test
using SumOfSquares
using DynamicPolynomials

function chebyshev_test(optimizer, config::MOIT.TestConfig)
    ndeg = 8   # Degree of Chebyshev polynomial

    @polyvar x

    Z = monomials((x,), 0:ndeg-1)

    model = _model(optimizer)

    @variable(model, γ)
    @variable(model, p1, Poly(Z))

    p = p1 + γ * x^ndeg # the leading coeff of p is γ

    dom = @set x >= -1 && x <= 1
    @constraint(model, 1 - p in SOSCone(), domain = dom)
    @constraint(model, p + 1 in SOSCone(), domain = dom)

    @objective(model, Max, γ)

    JuMP.optimize!(model)

    @test JuMP.termination_status(model) == MOI.OPTIMAL
    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT

    @test isapprox(JuMP.value(p), 128x^8 - 256x^6 + 160x^4 - 32x^2 + 1, ztol=4config.atol, atol=4config.atol, rtol=4config.rtol)
    @test isapprox(JuMP.value(γ), 128, atol=config.atol, rtol=config.rtol)
end
sd_tests["chebyshev"] = chebyshev_test


text:
#
# Protein - TMAO (compare new and old implementations)
#

using ComplexMixtures
using PDBTools
using Plots
using DelimitedFiles
ENV["GKSwstype"] = "nul"

# Here we use the PDBTools package to read the pdb file (from http://github.com/m3g/PDBTools)
atoms = PDBTools.readPDB("../../test/data/NAMD/structure.pdb")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
protein = PDBTools.select(atoms,"protein")
solute = ComplexMixtures.Selection(protein,nmols=1)

# The solvent is TMAO, which has 14 atoms. Use the natomspermol to indicate how many
# atoms each molecule has, such that there is no ambiguity on how to split the coordinates
# of the selection into individual molecules.
tmao = PDBTools.select(atoms,"resname TMAO")
solvent = ComplexMixtures.Selection(tmao,natomspermol=14)

# Initialize trajectroy data structure and open input stream
trajectory = ComplexMixtures.Trajectory("../../test/data/NAMD/trajectory.dcd",solute,solvent)

# Input options for the calculation
options = ComplexMixtures.Options(binstep=0.2)

# Run ComplexMixtures calculation, and get the results in the R structure
R = ComplexMixtures.mddf(trajectory,options)

# Save data for future loeading
ComplexMixtures.save(R,"example1.json")

# Save data in human-readable format
ComplexMixtures.write(R,"example1.dat",solute,solvent)

old = readdlm("./old_reference.dat",comments=true,comment_char='#')

plot(layout=(6,1))

sp=1
plot!(ylabel="mddf",subplot=sp)
plot!(old[:,1],old[:,2],subplot=sp,label="old")
plot!(R.d,R.mddf,subplot=sp,label="new - mddf")
plot!(R.d,R.rdf,subplot=sp,label="new - rdf")
plot!(legend=:topright,subplot=sp)

sp=2
plot!(ylabel="KB",subplot=sp)
plot!(old[:,1],old[:,3],subplot=sp,label="old")
plot!(R.d,R.kb,subplot=sp,label="new - mddf")
plot!(R.d,R.kb_rdf,subplot=sp,label="new - rdf")
plot!(legend=:topright,subplot=sp)

sp=3
plot!(ylabel="Count",subplot=sp)
plot!(old[:,1],old[:,4],subplot=sp,label="old")
scatter!(R.d,R.md_count,subplot=sp,label="new")
plot!(old[:,1],old[:,5],subplot=sp,label="old - rand")
scatter!(R.d,R.md_count_random,subplot=sp,label="new -rand")

sp=4
plot!(ylabel="Shell vol", subplot=sp)
plot!(old[:,1],old[:,8],subplot=sp,label="old")
plot!(R.d,R.volume.shell,subplot=sp,label="new")

sp=5
plot!(ylabel="Sum MD", subplot=sp)
plot!(old[:,1],old[:,6],subplot=sp,label="old - md")
scatter!(R.d,R.sum_md_count,subplot=sp,label="new - md")

sp=6
plot!(ylabel="Sum RAND", subplot=sp)
plot!(old[:,1],old[:,7],subplot=sp,label="old - rand")
scatter!(R.d,R.sum_md_count_random,subplot=sp,label="new - rand")
scatter!(R.d,R.sum_rdf_count,subplot=sp,label="new - rdf")
plot!(legend=:topleft,subplot=sp)

#sp=7
#plot!(ylabel="Count",subplot=sp)
#z = zeros(R.nbins)
#@. z = old[:,5]/R.md_count_random
#plot!(R.d,z,subplot=sp,label="old/rand")
#println(z[R.nbins])

plot!(size=(800,1300))
savefig("./example1.pdf")


text:
using Test
using Liga

include("testobjects.jl")

include("testoperations.jl")


text:
module Quad_Types

using Aerospace

export
    MotorCmd,
    Motors,
    SoftwareTelem,
    SoftwareMsg,
    Quadrotor,
    QuadForces

type MotorCmd
    motor1_cmd::Float64
    motor2_cmd::Float64
    motor3_cmd::Float64
    motor4_cmd::Float64
end

type Motors
    motor1::Actuator
    motor2::Actuator
    motor3::Actuator
    motor4::Actuator
end

type SoftwareTelem
    time::Float64
    navstates::NavStates
    axcmd::Float64
    aycmd::Float64
    azcmd::Float64
    pcmd::Float64
    qcmd::Float64
    rcmd::Float64
end

type SoftwareMsg
    time::Float64
    dt::Float64
    eta::Array{Float64,1}
    dz::Float64
    count::UInt;
    mode::UInt
    softwareTM::SoftwareTelem
    motorcmd::MotorCmd
end

type Quadrotor
    time::Float64
    dt::Float64
    airframe::Aero6DOF
    motors::Motors
    motorCmd::MotorCmd
    software::SoftwareMsg
end

type QuadForces
    f1body::Array{Float64,1}
    f2body::Array{Float64,1}
    f3body::Array{Float64,1}
    f4body::Array{Float64,1}
    m1body::Array{Float64,1}
    m2body::Array{Float64,1}
    m3body::Array{Float64,1}
    m4body::Array{Float64,1}
end

end ## End of Module


text:
using NullableArrays
using Lift

srand(1)
N = 5_000_000
f(x, y) = x * y
g(x, y) = x + y

X = NullableArray(rand(N))
Y = NullableArray(rand(N))

Xn = NullableArray(rand(N), rand(Bool, N))
Yn = NullableArray(rand(N), rand(Bool, N))

function tracedot(X, Y)
    res = Nullable(0.0)
    for i in 1:length(X)
         res += @^ f(X[i], Y[i]) Float64
    end
    return res
end

function profile_no_nulls()
    tracedot(X, Y)
    # tracedot(X, Y)
    @time(tracedot(X, Y))
end

function profile_half_nulls()
    tracedot(Xn, Yn)
    # tracedot(Xn, Yn)
    @time(tracedot(Xn, Yn))
end


text:
include("EmbeddingTest.jl")
include("../src/SequenceTools.jl")

@time output = SequencePad(output, TimeSteps)


text:
"""
get_width_height_35(image)
--> (w, h)::Quantity{Length}

'image' can be of types
- SVGimage (output from 'readsvg'),
- CairoSurfaceBase{UInt32} (output from 'readpng')
- matrix  - it is assumed that matrices follow the image convention: A column represent horizontal pixels.

The width and height of image based on scale_to_pt(m).
Note that the width can be inaccurate, as this is ambiguous for svg input.
"""
get_width_height_35(image) =  (scale_pt_to_unit(m)) .* (image.width, image.height)
get_width_height_35(image::Matrix) =  (scale_pt_to_unit(m)) .* reverse(size(matrix))

"""
    place_image_35(image::SVGimage, pos::Point;
                  width = missing, height = missing,
                  centered = false)
    -> (upper left point, lower right point)

Return the upper left and bottom right points of the placed and scaled image.
"""
function place_image_35(pos::Point, image::SVGimage;
                      width = missing, height = missing,
                      centered = false)
    @assert !ismissing(width) + !ismissing(height) < 2 "Width and height can not be specified simultaneously."
    original_width, original_height = get_width_height_35(image)
    # Find scaling from input to output
    scalefac = if ismissing(width) && ismissing(height)
        1.0
    elseif ismissing(height)
        scale_to_pt(width) / scale_to_pt(original_width)
    elseif ismissing(width)
        scale_to_pt(height) / scale_to_pt(original_height)
    end
    # Destination size
    dest_width_pix, dest_height_pix = scale_to_pt.(scalefac .* (original_width, original_height))
    # Destination upper left corner
    ptupleft = centered ?  pos - 0.5 .* (dest_width_pix, dest_height_pix) : pos
    @layer begin
        translate(ptupleft)
        scale(scalefac)
        placeimage(image, Point(0, 0); centered = false)
    end
    ptupleft, ptupleft + (dest_width_pix, dest_height_pix)
end


text:
# Examples presented in class - Lecture 1
# LP Resource allocation

using JuMP, Cbc
# JuMP is for implementing math. programming models;
# Cbc is for solving them.

# Example 1 - resource allocation

# Problem data
i = 1:2 # i=1: Seattle; i=2: San Diego
j = 1:3 # j=1: New York; j=2: Chicago; j=3: Miami

C = [350 600] # Capacities of the factories
D = [325 300 275] # Demand of clients
T = [2.5 1.7 1.8
     3.5 1.9 1.4] # Transportation costs

# Model implementation
# Creates a model and informs the solver to be used
m = Model(solver = CbcSolver())

# Decision variable for the total transported
@variable(m, x[i,j] >= 0)

# Capacity constraint
@constraint(m, cap[i = 1:2], sum(x[i,j] for j = 1:3) <= C[i])
# Demand constraint
@constraint(m, dem[j = 1:3], sum(x[i,j] for i = 1:2) >= D[j])

# Total distribution cost that we want to minimise
@objective(m, Min, sum(T[i,j]*x[i,j] for i=1:2, j=1:3))

println(m) # Prints the mathematical model for debugging
solve(m) # Solve the model

# Prints the optimal solution
println("\nDistribution plan: \n", getvalue(x))


text:
"Generate the non-linear system for a lower principal representation for odd n."
function LowerPrincipalOdd(s::Dictionary, moments = compute_moments(s))
    n = length(s) - 1
    @assert isodd(n)
    l = (n+1) >> 1
    QuadRuleFreePoints(QuadRuleData(s, l, moments))
end

"Generate the non-linear system for an upper principal representation for odd n."
function UpperPrincipalOdd(s::Dictionary, moments = compute_moments(s))
    n = length(s) - 1
    @assert isodd(n)
    l = ((n+1) >> 1) + 1
    QuadRuleFixedPoints(QuadRuleData(s, l, moments), [1, l], [supportleft(s), supportright(s)])
end

"Generate the non-linear system for a lower principal representation for even n."
function LowerPrincipalEven(s::Dictionary, moments = compute_moments(s))
    n = length(s) - 1
    @assert iseven(n)
    l = (n >> 1) + 1
    QuadRuleFixedPoints(QuadRuleData(s, l, moments), [1], [supportleft(s)])
end

"Generate the non-linear system for an upper principal representation for even n."
function UpperPrincipalEven(s::Dictionary, moments = compute_moments(s))
    n = length(s) - 1
    @assert iseven(n)
    l = (n >> 1) + 1
    QuadRuleFixedPoints(QuadRuleData(s, l, moments), [l], [supportright(s)])
end

# The fixed root is in K1 = [a,t1].
function CanonicalRepresentationOdd_K1(s::Dictionary, xstar, moments = compute_moments(s))
    n = length(s) - 1
    @assert isodd(n)
    l = ((n+1) >> 1) + 1
    QuadRuleFixedPoints(QuadRuleData(s, l, moments), [1, l], [xstar, supportright(s)])
end

# The fixed root is in J1 = [t1,s2].
function CanonicalRepresentationOdd_J1(s::Dictionary, xstar, moments = compute_moments(s))
    n = length(s) - 1
    @assert isodd(n)
    l = ((n+1) >> 1) + 1
    QuadRuleFixedPoints(QuadRuleData(s, l, moments), [1, 2], [supportleft(s), xstar])
end

# The fixed root is in J1 = [a,s1].
function CanonicalRepresentationEven_J1(s::Dictionary, xstar, moments = compute_moments(s))
    n = length(s) - 1
    @assert iseven(n)
    l = (n >> 1) + 1
    QuadRuleFixedPoints(QuadRuleData(s, l, moments), [1], [xstar])
end

# The fixed root is in K1 = [s1,t2].
function CanonicalRepresentationEven_K1(s::Dictionary, xstar, moments = compute_moments(s))
    n = length(s) - 1
    @assert iseven(n)
    l = (n >> 1) + 2
    QuadRuleFixedPoints(QuadRuleData(s, l, moments), [1, 2, l], [supportleft(s), xstar, supportright(s)])
end


text:
# ------------ Gradients and parameters ------------ #

const Gradients = DynamicMap{Value}

@inline function accumulate!(gs::Gradients, addr, v)
    set_sub!(gs, addr, Value(v + get(gs, addr, zero(v))))
end

@inline function accumulate!(gs::Gradients, addr, v::Value)
    val = get_value(v)
    set_sub!(gs, addr, Value(val + get(gs, addr, zero(val))))
end

@inline function accumulate!(gs1::Gradients, gs2::Gradients)
    for k in setdiff(keys(gs2.tree), keys(gs1.tree))
        accumulate!(gs1, k, get_sub(gs2, k))
    end
    inter = intersect(keys(gs1.tree), keys(gs2.tree))
    for k in inter
        accumulate!(gs1, k, get_sub(gs2, k))
    end
end

const LearnableByAddress = DynamicMap{Value}
Parameters() = LearnableByAddress()

function learnables(arr::Array{Pair{T, K}}) where {T <: Tuple, K}
    top = LearnableByAddress()
    map(arr) do (k, v)
        set_sub!(top, k, Value(v))
    end
    return top
end

function learnables(p::Pair{T, K}) where {T <: Symbol, K <: AddressMap}
    top = LearnableByAddress()
    set_sub!(top, p[1], p[2])
    return top
end

# ------------ update_learnables links into Flux optimiser APIs ------------ #

function update_learnables(opt, a::LearnableByAddress, b::Gradients)
    p_arr = array(a, Float64)
    gs_arr = array(b, Float64)
    update!(opt, p_arr, -gs_arr)
    return target(a, p_arr)
end


text:
using OpenGene
using OpenGene.Algorithm
using OpenGene.Reference
using Base.Test

include("Base/Base.jl")
include("IO/IO.jl")
include("Algorithm/Algorithm.jl")
include("Reference/Reference.jl")


text:
# custom wrappers for BLAS and LAPACK routines, together with some custom definitions
using LinearAlgebra: BlasFloat, Char, BlasInt, LAPACK, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare,
    triu!
using LinearAlgebra.BLAS: @blasfunc, libblas, BlasReal, BlasComplex
using LinearAlgebra.LAPACK: liblapack, chklapackerror

set_num_blas_threads(n::Integer) = LinearAlgebra.BLAS.set_num_threads(n)

# TODO: Remove once we drop support for Julia 1.5 or less
function get_num_blas_threads()
	blas = LinearAlgebra.BLAS.vendor()
	# Wrap in a try to catch unsupported blas versions
	if blas == :openblas
		return ccall((:openblas_get_num_threads, Base.libblas_name), Cint, ())
	elseif blas == :openblas64
		return ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
	elseif blas == :mkl
		return ccall((:mkl_get_max_threads, Base.libblas_name), Cint, ())
	end

	# OSX BLAS looks at an environment variable
	@static if Sys.isapple()
		return tryparse(Cint, get(ENV, "VECLIB_MAXIMUM_THREADS", "1"))
	end
    error("unknown BLAS")
end

# TODO: define for CuMatrix if we support this
function _one!(A::DenseMatrix)
    Threads.@threads for j = 1:size(A, 2)
        @simd for i = 1:size(A, 1)
            @inbounds A[i, j] = i == j
        end
    end
    return A
end

# MATRIX factorizations
#-----------------------
abstract type FactorizationAlgorithm end
abstract type OrthogonalFactorizationAlgorithm <: FactorizationAlgorithm end

struct QRpos <: OrthogonalFactorizationAlgorithm
end
struct QR <: OrthogonalFactorizationAlgorithm
end
struct QL <: OrthogonalFactorizationAlgorithm
end
struct QLpos <: OrthogonalFactorizationAlgorithm
end
struct LQ <: OrthogonalFactorizationAlgorithm
end
struct LQpos <: OrthogonalFactorizationAlgorithm
end
struct RQ <: OrthogonalFactorizationAlgorithm
end
struct RQpos <: OrthogonalFactorizationAlgorithm
end
struct SVD <: OrthogonalFactorizationAlgorithm
end
struct Polar <: OrthogonalFactorizationAlgorithm
end
struct SDD <: OrthogonalFactorizationAlgorithm # lapack's default divide and conquer algorithm
end

Base.adjoint(::QRpos) = LQpos()
Base.adjoint(::QR) = LQ()
Base.adjoint(::LQpos) = QRpos()
Base.adjoint(::LQ) = QR()

Base.adjoint(::QLpos) = RQpos()
Base.adjoint(::QL) = RQ()
Base.adjoint(::RQpos) = QLpos()
Base.adjoint(::RQ) = QL()

Base.adjoint(alg::Union{SVD, SDD, Polar}) = alg

_safesign(s::Real) = ifelse(s<zero(s), -one(s), +one(s))
_safesign(s::Complex) = ifelse(iszero(s), one(s), s/abs(s))

function _leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{QR, QRpos}, atol::Real)
    iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
    m, n = size(A)
    k = min(m, n)
    A, T = LAPACK.geqrt!(A, min(minimum(size(A)), 36))
    Q = similar(A, m, k)
    for j = 1:k
        for i = 1:m
            Q[i, j] = i == j
        end
    end
    Q = LAPACK.gemqrt!('L', 'N', A, T, Q)
    R = triu!(A[1:k, :])

    if isa(alg, QRpos)
        @inbounds for j = 1:k
            s = _safesign(R[j,j])
            @simd for i = 1:m
                Q[i,j] *= s
            end
        end
        @inbounds for j = size(R, 2):-1:1
            for i = 1:min(k, j)
                R[i,j] = R[i,j]*conj(_safesign(R[i,i]))
            end
        end
    end
    return Q, R
end
# TODO: reconsider the following implementation
# Unfortunately, geqrfp seems a bit slower than geqrt in the intermediate region
# around matrix size 100, which is the interesting region. => Investigate and maybe fix
# function _leftorth!(A::StridedMatrix{<:BlasFloat})
#     m, n = size(A)
#     A, τ = geqrfp!(A)
#     Q = LAPACK.ormqr!('L', 'N', A, τ, eye(eltype(A), m, min(m, n)))
#     R = triu!(A[1:min(m, n), :])
#     return Q, R
# end

function _leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{QL, QLpos}, atol::Real)
    iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
    m, n = size(A)
    @assert m >= n

    nhalf = div(n, 2)
    #swap columns in A
    @inbounds for j = 1:nhalf, i = 1:m
        A[i,j], A[i,n+1-j] = A[i,n+1-j], A[i,j]
    end
    Q, R = _leftorth!(A, isa(alg, QL) ? QR() : QRpos() , atol)

    #swap columns in Q
    @inbounds for j = 1:nhalf, i = 1:m
        Q[i,j], Q[i,n+1-j] = Q[i,n+1-j], Q[i,j]
    end
    #swap rows and columns in R
    @inbounds for j = 1:nhalf, i = 1:n
        R[i,j], R[n+1-i,n+1-j] = R[n+1-i,n+1-j], R[i,j]
    end
    if isodd(n)
        j = nhalf+1
        @inbounds for i = 1:nhalf
            R[i,j], R[n+1-i,j] = R[n+1-i,j], R[i,j]
        end
    end
    return Q, R
end

function _leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD, SDD, Polar}, atol::Real)
    U, S, V = alg isa SVD ? LAPACK.gesvd!('S', 'S', A) : LAPACK.gesdd!('S', A)
    if isa(alg, Union{SVD, SDD})
        n = count(s-> s .> atol, S)
        if n != length(S)
            return U[:,1:n], lmul!(Diagonal(S[1:n]), V[1:n, :])
        else
            return U, lmul!(Diagonal(S), V)
        end
    else
        iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
        # TODO: check Lapack to see if we can recycle memory of A
        Q = mul!(A, U, V)
        Sq = map!(sqrt, S, S)
        SqV = lmul!(Diagonal(Sq), V)
        R = SqV'*SqV
        return Q, R
    end
end

function _leftnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{QR, QRpos}, atol::Real)
    iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
    m, n = size(A)
    m >= n || throw(ArgumentError("no null space if less rows than columns"))

    A, T = LAPACK.geqrt!(A, min(minimum(size(A)), 36))
    N = similar(A, m, max(0, m-n));
    fill!(N, 0)
    for k = 1:m-n
        N[n+k,k] = 1
    end
    N = LAPACK.gemqrt!('L', 'N', A, T, N)
end

function _leftnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD, SDD}, atol::Real)
    size(A, 2) == 0 && return _one!(similar(A, (size(A, 1), size(A, 1))))
    U, S, V = alg isa SVD ? LAPACK.gesvd!('A', 'N', A) : LAPACK.gesdd!('A', A)
    indstart = count(>(atol), S) + 1
    return U[:, indstart:end]
end

function _rightorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{LQ, LQpos, RQ, RQpos},
                        atol::Real)
    iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
    # TODO: geqrfp seems a bit slower than geqrt in the intermediate region around
    # matrix size 100, which is the interesting region. => Investigate and fix
    m, n = size(A)
    k = min(m, n)
    At = transpose!(similar(A, n, m), A)

    if isa(alg, RQ) || isa(alg, RQpos)
        @assert m <= n

        mhalf = div(m, 2)
        # swap columns in At
        @inbounds for j = 1:mhalf, i = 1:n
            At[i,j], At[i,m+1-j] = At[i,m+1-j], At[i,j]
        end
        Qt, Rt = _leftorth!(At, isa(alg, RQ) ? QR() : QRpos(), atol)

        @inbounds for j = 1:mhalf, i = 1:n
            Qt[i,j], Qt[i,m+1-j] = Qt[i,m+1-j], Qt[i,j]
        end
        @inbounds for j = 1:mhalf, i = 1:m
            Rt[i,j], Rt[m+1-i,m+1-j] = Rt[m+1-i,m+1-j], Rt[i,j]
        end
        if isodd(m)
            j = mhalf+1
            @inbounds for i = 1:mhalf
                Rt[i,j], Rt[m+1-i,j] = Rt[m+1-i,j], Rt[i,j]
            end
        end
        Q = transpose!(A, Qt)
        R = transpose!(similar(A, (m, m)), Rt) # TODO: efficient in place
        return R, Q
    else
        Qt, Lt = _leftorth!(At, alg', atol)
        if m > n
            L = transpose!(A, Lt)
            Q = transpose!(similar(A, (n, n)), Qt) # TODO: efficient in place
        else
            Q = transpose!(A, Qt)
            L = transpose!(similar(A, (m, m)), Lt) # TODO: efficient in place
        end
        return L, Q
    end
end

function _rightorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD, SDD, Polar}, atol::Real)
    U, S, V = alg isa SVD ? LAPACK.gesvd!('S', 'S', A) : LAPACK.gesdd!('S', A)
    if isa(alg, Union{SVD, SDD})
        n = count(s-> s .> atol, S)
        if n != length(S)
            return rmul!(U[:,1:n], Diagonal(S[1:n])), V[1:n,:]
        else
            return rmul!(U, Diagonal(S)), V
        end
    else
        iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
        Q = mul!(A, U, V)
        Sq = map!(sqrt, S, S)
        USq = rmul!(U, Diagonal(Sq))
        L = USq*USq'
        return L, Q
    end
end

function _rightnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{LQ, LQpos}, atol::Real)
    iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
    m, n = size(A)
    k = min(m, n)
    At = adjoint!(similar(A, n, m), A)
    At, T = LAPACK.geqrt!(At, min(k, 36))
    N = similar(A, max(n-m, 0), n);
    fill!(N, 0)
    for k = 1:n-m
        N[k,m+k] = 1
    end
    N = LAPACK.gemqrt!('R', eltype(At) <: Real ? 'T' : 'C', At, T, N)
end

function _rightnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD, SDD}, atol::Real)
    size(A, 1) == 0 && return _one!(similar(A, (size(A, 2), size(A, 2))))
    U, S, V = alg isa SVD ? LAPACK.gesvd!('N', 'A', A) : LAPACK.gesdd!('A', A)
    indstart = count(>(atol), S) + 1
    return V[indstart:end, :]
end

function _svd!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD, SDD})
    U, S, V = alg isa SVD ? LAPACK.gesvd!('S', 'S', A) : LAPACK.gesdd!('S', A)
    return U, S, V
end

# TODO: override Julia's eig interface
# eig!(A::StridedMatrix{<:BlasFloat}) = LinAlg.LAPACK.gees!('V', A)
# function eig!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasReal
#     n = size(A, 2)
#     n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
#     issymmetric(A) && return eigfact!(Symmetric(A))
#     A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
#     evec = zeros(Complex{T}, n, n)
#     j = 1
#     while j <= n
#         if WI[j] == 0
#             evec[:,j] = view(VR, :, j)
#         else
#             for i = 1:n
#                 evec[i,j]   = VR[i,j] + im*VR[i,j+1]
#                 evec[i,j+1] = VR[i,j] - im*VR[i,j+1]
#             end
#             j += 1
#         end
#         j += 1
#     end
#     return Eigen(complex.(WR, WI), evec)
# end
#
# function eigfact!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasComplex
#     n = size(A, 2)
#     n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
#     ishermitian(A) && return eigfact!(Hermitian(A))
#     return Eigen(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]...)
# end
#
#
# eigfact!(A::RealHermSymComplexHerm{<:BlasReal, <:StridedMatrix}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...)
#
#
#
# Modified / missing Lapack wrappers
#------------------------------------
# geqrfp!: computes qrpos factorization, missing in Base
geqrfp!(A::StridedMatrix{<:BlasFloat}) =
    ((m, n) = size(A); geqrfp!(A, similar(A, min(m, n))))

for (geqrfp, elty, relty) in
    ((:dgeqrfp_, :Float64, :Float64), (:sgeqrfp_, :Float32, :Float32),
        (:zgeqrfp_, :ComplexF64, :Float64), (:cgeqrfp_, :ComplexF32, :Float32))
    @eval begin
        function geqrfp!(A::StridedMatrix{$elty}, tau::StridedVector{$elty})
            chkstride1(A, tau)
            m, n  = size(A)
            if length(tau) != min(m, n)
                throw(DimensionMismatch("tau has length $(length(tau)), but needs length $(min(m, n))"))
            end
            work  = Vector{$elty}(1)
            lwork = BlasInt(-1)
            info  = Ref{BlasInt}()
            for i = 1:2                # first call returns lwork as work[1]
                ccall((@blasfunc($geqrfp), liblapack), Nothing,
                      (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
                       Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
                      Ref(m), Ref(n), A, Ref(max(1, stride(A, 2))),
                      tau, work, Ref(lwork), info)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                end
            end
            A, tau
        end
    end
end


text:
# Automatically generated using Clang.jl
#! format: off

const FILE = Cvoid

# Skipping MacroDefinition: CEED_QFUNCTION ( name ) static const char name ## _loc [ ] = __FILE__ ":" # name ; static int name

# Skipping MacroDefinition: CeedError ( ceed , ecode , ... ) ( CeedErrorImpl ( ( ceed ) , __FILE__ , __LINE__ , __func__ , ( ecode ) , __VA_ARGS__ ) ? : ( ecode ) )

const CeedInt = Int32
const CeedScalar = Cdouble
const Ceed_private = Cvoid
const Ceed = Ptr{Ceed_private}
const CeedRequest_private = Cvoid
const CeedRequest = Ptr{CeedRequest_private}
const CeedVector_private = Cvoid
const CeedVector = Ptr{CeedVector_private}
const CeedElemRestriction_private = Cvoid
const CeedElemRestriction = Ptr{CeedElemRestriction_private}
const CeedBasis_private = Cvoid
const CeedBasis = Ptr{CeedBasis_private}
const CeedQFunction_private = Cvoid
const CeedQFunction = Ptr{CeedQFunction_private}
const CeedQFunctionContext_private = Cvoid
const CeedQFunctionContext = Ptr{CeedQFunctionContext_private}
const CeedOperator_private = Cvoid
const CeedOperator = Ptr{CeedOperator_private}

@cenum CeedMemType::UInt32 begin
    CEED_MEM_HOST = 0
    CEED_MEM_DEVICE = 1
end

@cenum CeedCopyMode::UInt32 begin
    CEED_COPY_VALUES = 0
    CEED_USE_POINTER = 1
    CEED_OWN_POINTER = 2
end

@cenum CeedNormType::UInt32 begin
    CEED_NORM_1 = 0
    CEED_NORM_2 = 1
    CEED_NORM_MAX = 2
end

@cenum CeedTransposeMode::UInt32 begin
    CEED_NOTRANSPOSE = 0
    CEED_TRANSPOSE = 1
end

@cenum CeedEvalMode::UInt32 begin
    CEED_EVAL_NONE = 0
    CEED_EVAL_INTERP = 1
    CEED_EVAL_GRAD = 2
    CEED_EVAL_DIV = 4
    CEED_EVAL_CURL = 8
    CEED_EVAL_WEIGHT = 16
end

@cenum CeedQuadMode::UInt32 begin
    CEED_GAUSS = 0
    CEED_GAUSS_LOBATTO = 1
end

@cenum CeedElemTopology::UInt32 begin
    CEED_LINE = 65536
    CEED_TRIANGLE = 131073
    CEED_QUAD = 131074
    CEED_TET = 196611
    CEED_PYRAMID = 196612
    CEED_PRISM = 196613
    CEED_HEX = 196614
end

const CeedQFunctionUser = Ptr{Cvoid}

# Skipping MacroDefinition: CEED_INTERN CEED_EXTERN __attribute__ ( ( visibility ( "hidden" ) ) )

const CEED_MAX_RESOURCE_LEN = 1024
const CEED_ALIGN = 64
const CEED_COMPOSITE_MAX = 16
const CEED_EPSILON = 1.0e-16
const CEED_DEBUG_COLOR = 0

# Skipping MacroDefinition: CeedDebug1 ( ceed , format , ... ) CeedDebugImpl ( ceed , format , ## __VA_ARGS__ )
# Skipping MacroDefinition: CeedDebug256 ( ceed , color , ... ) CeedDebugImpl256 ( ceed , color , ## __VA_ARGS__ )
# Skipping MacroDefinition: CeedDebug ( ... ) CeedDebug256 ( ceed , ( unsigned char ) CEED_DEBUG_COLOR , ## __VA_ARGS__ )
# Skipping MacroDefinition: CeedChk ( ierr ) do { if ( ierr ) return ierr ; } while ( 0 )
# Skipping MacroDefinition: CeedMalloc ( n , p ) CeedMallocArray ( ( n ) , sizeof ( * * ( p ) ) , p )
# Skipping MacroDefinition: CeedCalloc ( n , p ) CeedCallocArray ( ( n ) , sizeof ( * * ( p ) ) , p )
# Skipping MacroDefinition: CeedRealloc ( n , p ) CeedReallocArray ( ( n ) , sizeof ( * * ( p ) ) , p )

const CeedTensorContract_private = Cvoid
const CeedTensorContract = Ptr{CeedTensorContract_private}
const CeedQFunctionField_private = Cvoid
const CeedQFunctionField = Ptr{CeedQFunctionField_private}
const CeedOperatorField_private = Cvoid
const CeedOperatorField = Ptr{CeedOperatorField_private}


text:
# This file was generated, do not modify it. # hide
using PyPlot
ioff() # hide

figure(figsize=(8,6))
plot(mpg)

savefig(joinpath(@OUTPUT, "ISL-lab-2-mpg.svg")) # hide


text:
# Julia wrapper for header: /home/shindo/local-lemon/cuda-8.0/include/curand.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function curandCreateGenerator(generator,rng_type)
    check_curandstatus(ccall((:curandCreateGenerator,libcurand),curandStatus_t,(Ptr{curandGenerator_t},curandRngType_t),generator,rng_type))
end

function curandCreateGeneratorHost(generator,rng_type)
    check_curandstatus(ccall((:curandCreateGeneratorHost,libcurand),curandStatus_t,(Ptr{curandGenerator_t},curandRngType_t),generator,rng_type))
end

function curandDestroyGenerator(generator)
    check_curandstatus(ccall((:curandDestroyGenerator,libcurand),curandStatus_t,(curandGenerator_t,),generator))
end

function curandGetVersion(version)
    check_curandstatus(ccall((:curandGetVersion,libcurand),curandStatus_t,(Ptr{Cint},),version))
end

function curandGetProperty(_type,value)
    check_curandstatus(ccall((:curandGetProperty,libcurand),curandStatus_t,(libraryPropertyType,Ptr{Cint}),_type,value))
end

function curandSetStream(generator,stream)
    check_curandstatus(ccall((:curandSetStream,libcurand),curandStatus_t,(curandGenerator_t,cudaStream_t),generator,stream))
end

function curandSetPseudoRandomGeneratorSeed(generator,seed)
    check_curandstatus(ccall((:curandSetPseudoRandomGeneratorSeed,libcurand),curandStatus_t,(curandGenerator_t,Culonglong),generator,seed))
end

function curandSetGeneratorOffset(generator,offset)
    check_curandstatus(ccall((:curandSetGeneratorOffset,libcurand),curandStatus_t,(curandGenerator_t,Culonglong),generator,offset))
end

function curandSetGeneratorOrdering(generator,order)
    check_curandstatus(ccall((:curandSetGeneratorOrdering,libcurand),curandStatus_t,(curandGenerator_t,curandOrdering_t),generator,order))
end

function curandSetQuasiRandomGeneratorDimensions(generator,num_dimensions)
    check_curandstatus(ccall((:curandSetQuasiRandomGeneratorDimensions,libcurand),curandStatus_t,(curandGenerator_t,UInt32),generator,num_dimensions))
end

function curandGenerate(generator,outputPtr,num)
    check_curandstatus(ccall((:curandGenerate,libcurand),curandStatus_t,(curandGenerator_t,Ptr{UInt32},Csize_t),generator,outputPtr,num))
end

function curandGenerateLongLong(generator,outputPtr,num)
    check_curandstatus(ccall((:curandGenerateLongLong,libcurand),curandStatus_t,(curandGenerator_t,Ptr{Culonglong},Csize_t),generator,outputPtr,num))
end

function curandGenerateUniform(generator,outputPtr,num)
    check_curandstatus(ccall((:curandGenerateUniform,libcurand),curandStatus_t,(curandGenerator_t,Ptr{Cfloat},Csize_t),generator,outputPtr,num))
end

function curandGenerateUniformDouble(generator,outputPtr,num)
    check_curandstatus(ccall((:curandGenerateUniformDouble,libcurand),curandStatus_t,(curandGenerator_t,Ptr{Cdouble},Csize_t),generator,outputPtr,num))
end

function curandGenerateNormal(generator,outputPtr,n,mean,stddev)
    check_curandstatus(ccall((:curandGenerateNormal,libcurand),curandStatus_t,(curandGenerator_t,Ptr{Cfloat},Csize_t,Cfloat,Cfloat),generator,outputPtr,n,mean,stddev))
end

function curandGenerateNormalDouble(generator,outputPtr,n,mean,stddev)
    check_curandstatus(ccall((:curandGenerateNormalDouble,libcurand),curandStatus_t,(curandGenerator_t,Ptr{Cdouble},Csize_t,Cdouble,Cdouble),generator,outputPtr,n,mean,stddev))
end

function curandGenerateLogNormal(generator,outputPtr,n,mean,stddev)
    check_curandstatus(ccall((:curandGenerateLogNormal,libcurand),curandStatus_t,(curandGenerator_t,Ptr{Cfloat},Csize_t,Cfloat,Cfloat),generator,outputPtr,n,mean,stddev))
end

function curandGenerateLogNormalDouble(generator,outputPtr,n,mean,stddev)
    check_curandstatus(ccall((:curandGenerateLogNormalDouble,libcurand),curandStatus_t,(curandGenerator_t,Ptr{Cdouble},Csize_t,Cdouble,Cdouble),generator,outputPtr,n,mean,stddev))
end

function curandCreatePoissonDistribution(lambda,discrete_distribution)
    check_curandstatus(ccall((:curandCreatePoissonDistribution,libcurand),curandStatus_t,(Cdouble,Ptr{curandDiscreteDistribution_t}),lambda,discrete_distribution))
end

function curandDestroyDistribution(discrete_distribution)
    check_curandstatus(ccall((:curandDestroyDistribution,libcurand),curandStatus_t,(curandDiscreteDistribution_t,),discrete_distribution))
end

function curandGeneratePoisson(generator,outputPtr,n,lambda)
    check_curandstatus(ccall((:curandGeneratePoisson,libcurand),curandStatus_t,(curandGenerator_t,Ptr{UInt32},Csize_t,Cdouble),generator,outputPtr,n,lambda))
end

function curandGeneratePoissonMethod(generator,outputPtr,n,lambda,method)
    check_curandstatus(ccall((:curandGeneratePoissonMethod,libcurand),curandStatus_t,(curandGenerator_t,Ptr{UInt32},Csize_t,Cdouble,curandMethod_t),generator,outputPtr,n,lambda,method))
end

function curandGenerateBinomial(generator,outputPtr,num,n,p)
    check_curandstatus(ccall((:curandGenerateBinomial,libcurand),curandStatus_t,(curandGenerator_t,Ptr{UInt32},Csize_t,UInt32,Cdouble),generator,outputPtr,num,n,p))
end

function curandGenerateBinomialMethod(generator,outputPtr,num,n,p,method)
    check_curandstatus(ccall((:curandGenerateBinomialMethod,libcurand),curandStatus_t,(curandGenerator_t,Ptr{UInt32},Csize_t,UInt32,Cdouble,curandMethod_t),generator,outputPtr,num,n,p,method))
end

function curandGenerateSeeds(generator)
    check_curandstatus(ccall((:curandGenerateSeeds,libcurand),curandStatus_t,(curandGenerator_t,),generator))
end

function curandGetDirectionVectors32(vectors,set)
    check_curandstatus(ccall((:curandGetDirectionVectors32,libcurand),curandStatus_t,(Ptr{Ptr{curandDirectionVectors32_t}},curandDirectionVectorSet_t),vectors,set))
end

function curandGetScrambleConstants32(constants)
    check_curandstatus(ccall((:curandGetScrambleConstants32,libcurand),curandStatus_t,(Ptr{Ptr{UInt32}},),constants))
end

function curandGetDirectionVectors64(vectors,set)
    check_curandstatus(ccall((:curandGetDirectionVectors64,libcurand),curandStatus_t,(Ptr{Ptr{curandDirectionVectors64_t}},curandDirectionVectorSet_t),vectors,set))
end

function curandGetScrambleConstants64(constants)
    check_curandstatus(ccall((:curandGetScrambleConstants64,libcurand),curandStatus_t,(Ptr{Ptr{Culonglong}},),constants))
end


text:
struct DDEProblem{uType,tType,lType,lType2,isinplace,P,F,H,K} <:
                          AbstractDDEProblem{uType,tType,lType,isinplace}
  f::F
  u0::uType
  h::H
  tspan::tType
  p::P
  constant_lags::lType
  dependent_lags::lType2
  kwargs::K
  neutral::Bool
  order_discontinuity_t0::Int

  @add_kwonly function DDEProblem{iip}(f::AbstractDDEFunction{iip}, u0, h, tspan, p=NullParameters();
                                       constant_lags = (),
                                       dependent_lags = (),
                                       neutral = f.mass_matrix !== I && det(f.mass_matrix) != 1,
                                       order_discontinuity_t0 = 0,
                                       kwargs...) where {iip}
    _tspan = promote_tspan(tspan)
    new{typeof(u0),typeof(_tspan),typeof(constant_lags),typeof(dependent_lags),isinplace(f),
        typeof(p),typeof(f),typeof(h),typeof(kwargs)}(
          f, u0, h, _tspan, p, constant_lags, dependent_lags, kwargs, neutral,
          order_discontinuity_t0)
  end

  function DDEProblem{iip}(f::AbstractDDEFunction{iip}, h, tspan::Tuple, p=NullParameters();
                           order_discontinuity_t0 = 1, kwargs...) where iip
    DDEProblem{iip}(f, h(p, first(tspan)), h, tspan, p;
                    order_discontinuity_t0 = max(1, order_discontinuity_t0), kwargs...)
  end

  function DDEProblem{iip}(f, args...; kwargs...) where iip
    DDEProblem{iip}(convert(DDEFunction{iip}, f), args...; kwargs...)
  end
end

DDEProblem(f, args...; kwargs...) =
  DDEProblem(convert(DDEFunction, f), args...; kwargs...)

DDEProblem(f::AbstractDDEFunction, args...; kwargs...) =
  DDEProblem{isinplace(f)}(f, args...; kwargs...)


text: # AlphaShapes.jl


text:
# Set up the environment to run this example. Make sure you're within the folder that this
# file lives in.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using AbstractGPs, Plots, Random, Stheno
gr();

# Define model.
σ², ω, T = 1e0, 1.0, 25.0;
f = @gppp let
    f1 = periodic(GP(SEKernel()), ω)
    f2 = GP(0.1 * SEKernel())
    f3 = f1 + f2
end;

# Sample from marginal process to generate toy data.
rng = MersenneTwister(123456);
S = 25;
x = GPPPInput(:f3, range(0.0, T; length=300));
fx = f(x, σ²);
y = rand(rng, fx);

xp_ = range(-2.5, T + 2.5; length=500);

# Plots for stuff.
f₁′_plot, f₂′_plot, f₃′_plot = plot(), plot(), plot();

posterior_plot = plot();
approx_in_marginal_posterior_plot = plot();
approx_in_latents_posterior_plot = plot();

##################################### Exact Inference ######################################

# Compute the posterior processes, sample from them, and compute marginals.
@show logpdf(fx, y)
f′ = posterior(fx, y);
xp_f1 = GPPPInput(:f1, xp_);
xp_f2 = GPPPInput(:f2, xp_);
xp_f3 = GPPPInput(:f3, xp_);
xp = BlockData(xp_f1, xp_f2, xp_f3);
f′xp = rand(f′(xp, 1e-9), S);
f₁′xp, f₂′xp, f₃′xp = split(xp, f′xp);

items = [
    (xp_f1, f₁′xp, :red, "exact", f₁′_plot),
    (xp_f2, f₂′xp, :red, "", f₂′_plot),
    (xp_f3, f₃′xp, :red, "", f₃′_plot),
];

# Plot posterior marginals and samples.
for (x, samples, colour, name, plt) in items
    plot!(plt, xp_, f′(x); fillalpha=0.2, color=colour, label=name);
    plot!(plt, xp_, samples; alpha=0.3, color=colour, label="", linewidth=1);
end

# Plot observations.
scatter!(posterior_plot, x.x, y;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);

##################### Approximate inference with pseudo-data in f3 #########################

# Compute approximate posterior process.
M3 = 25;
z = GPPPInput(:f3, collect(range(0, T; length=M3)));
u = f(z, 1e-9);
f′_approx = approx_posterior(VFE(), fx, y, u);

@show elbo(fx, y, u);

items = [
    (xp_f1, :green, "Z in f3", f₁′_plot),
    (xp_f2, :green, "", f₂′_plot),
    (xp_f3, :green, "", f₃′_plot),
];

# Posterior marginal variance.
for (x, colour, name, plt) in items
    plot!(plt, xp_, f′_approx(x); fillalpha=0.3, fillcolor=colour, label=name);
end

# Plot observations and pseudo-input locations.
scatter!(approx_in_marginal_posterior_plot, x.x, y;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);
scatter!(approx_in_marginal_posterior_plot, z.x, zeros(M3);
    markercolor=:black,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z",
);

############# Perform approximate inference by placing pseudo-data in f1 and f2 ############

# Compute approximate posterior process.
M1, M2 = 15, 10;
z1 = GPPPInput(:f1, collect(range(0.0; step=1 / ω, length=M1)));
z2 = GPPPInput(:f2, collect(range(0.0; step=T, length=M2)));
z12 = BlockData(z1, z2);
u12 = f(z12, 1e-9);
f′_approx_12 = approx_posterior(VFE(), fx, y, u12);

@show elbo(fx, y, u12);

items = [
    (xp_f1, :blue, "Z in f1 and f2", f₁′_plot),
    (xp_f2, :blue, "", f₂′_plot),
    (xp_f3, :blue, "", f₃′_plot),
];

# Plot posterior marginals.
for (x, colour, name, plt) in items
    plot!(plt, xp_, f′_approx_12(x); fillalpha=0.3, fillcolor=colour, label=name)
end

# Plot observations and pseudo-input locations.
scatter!(approx_in_latents_posterior_plot, x.x, y;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);
scatter!(approx_in_latents_posterior_plot, z1.x, zeros(M1);
    markercolor=:black,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z₁",
);
scatter!(approx_in_latents_posterior_plot, z2.x, zeros(M2);
    markercolor=:magenta,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z₂",
);

# display(posterior_plot);
# display(approx_in_marginal_posterior_plot);
# display(approx_in_latents_posterior_plot);

plot(f₁′_plot, f₂′_plot, f₃′_plot; layout=(3, 1))


text:
function DfpDec(arg0::Dfp)
    return DfpDec((Dfp,), arg0)
end

function new_instance(obj::DfpDec)
    return jcall(obj, "newInstance", Dfp, ())
end

function new_instance(obj::DfpDec, arg0::Dfp)
    return jcall(obj, "newInstance", Dfp, (Dfp,), arg0)
end

function new_instance(obj::DfpDec, arg0::JString)
    return jcall(obj, "newInstance", Dfp, (JString,), arg0)
end

function new_instance(obj::DfpDec, arg0::jbyte)
    return jcall(obj, "newInstance", Dfp, (jbyte,), arg0)
end

function new_instance(obj::DfpDec, arg0::jbyte, arg1::jbyte)
    return jcall(obj, "newInstance", Dfp, (jbyte, jbyte), arg0, arg1)
end

function new_instance(obj::DfpDec, arg0::jdouble)
    return jcall(obj, "newInstance", Dfp, (jdouble,), arg0)
end

function new_instance(obj::DfpDec, arg0::jint)
    return jcall(obj, "newInstance", Dfp, (jint,), arg0)
end

function new_instance(obj::DfpDec, arg0::jlong)
    return jcall(obj, "newInstance", Dfp, (jlong,), arg0)
end

function next_after(obj::DfpDec, arg0::Dfp)
    return jcall(obj, "nextAfter", Dfp, (Dfp,), arg0)
end


text:
# run this with: julia1.5 --project="." examples/jmlr/single_pend_feas.jl "small" |& tee examples/jmlr/single_pend_feas_log.txt

using OVERTVerify
using LazySets
using Dates

controller_type = ARGS[1] # pass from command line, e.g. "small"
controller = "nnet_files/jmlr/single_pendulum_$(controller_type)_controller.nnet"
println("Controller: ", controller)
query = OvertQuery(
	SinglePendulum,    # problem
	controller,        # network file
	Id(),              # last layer activation layer Id()=linear, or ReLU()=relu
	"MIP",             # query solver, "MIP" or "ReluPlex"
	25,                # ntime
	0.1,               # dt
	-1,                # N_overt
	)

input_set = Hyperrectangle(low=[1., 0.], high=[1.2, 0.2])
target_set = HalfSpace([1., 0.], -0.2617) # 1*x_1 + 0*x_2 <= -.2167  -->  x_1 <= -.2167
t1 = Dates.time()
SATus, vals, stats = symbolic_satisfiability(query, input_set, target_set)
t2 = Dates.time()
dt = (t2-t1)

using JLD2
JLD2.@save "src/examples/jmlr/data/single_pendulum_satisfiability_$(controller_type)_controller_data.jld2" query input_set target_set SATus vals stats dt controller


text:
"""
    FDIModel <: AbstractFDDObject

Type for synthesis models employed to solve fault detection and isolation problems.

If `sysf::FDIModel` is the synthesis model object, the underlying descriptor system model
can be obtained via `sysf.sys` and the indices of control, disturbance, fault, noise
and auxiliary inputs can be accessed as the integer vectors
contained in `sysf.controls`, `sysf.disturbances`, `sysf.faults`,
`sysf.noise` and `sysf.aux`, respectively.
"""
struct FDIModel{T} <: AbstractFDDObject where T
    sys::DescriptorStateSpace{T}
    controls::Vector{Int}
    disturbances::Vector{Int}
    faults::Vector{Int}
    noise::Vector{Int}
    aux::Vector{Int}
    function FDIModel{T}(sys::DescriptorStateSpace{T}; controls::Vector{Int} = Int[], disturbances::Vector{Int} = Int[],
                         faults::Vector{Int} = Int[], noise::Vector{Int} = Int[], aux::Vector{Int} = Int[]) where T
        fdimodel_validation(sys, controls, disturbances, faults, noise, aux)
        new{T}(sys, controls, disturbances, faults, noise, aux)
    end
end
function fdimodel_validation(sys::DescriptorStateSpace{T}, controls::Vector{Int}, disturbances::Vector{Int},
                             faults::Vector{Int}, noise::Vector{Int}, aux::Vector{Int}) where T
    m = size(sys,2)
    isempty(intersect(controls,disturbances)) || error("control and disturbance inputs must be distinct")
    isempty(intersect(disturbances,faults)) || error("disturbance and fault inputs must be distinct")
    isempty(intersect(noise,faults)) || error("noise and fault inputs must be distinct")

    any(controls .> m) && error("control inputs indices larger than the number of system inputs $m")
    any(disturbances .> m) && error("disturbance inputs indices larger than the number of system inputs $m")
    any(faults .> m) && error("fault inputs indices larger than the number of system inputs $m")
    any(noise .> m) && error("noise inputs indices larger than the number of system inputs $m")
    any(aux .> m) && error("auxiliary inputs indices larger than the number of system inputs $m")
end
FDIModel(sys::DescriptorStateSpace, controls::Vector{Int}, disturbances::Vector{Int}, faults::Vector{Int}, noise::Vector{Int}, aux::Vector{Int})  =
           FDIModel{eltype(sys)}(sys; controls, disturbances, faults, noise, aux)

"""
    FDFilter <: AbstractFDDObject

Type for fault detection filters resulted as solutions of fault detection problems.

If `filter::FDFilter` is the fault detection filter object, the underlying descriptor system model
can be obtained via `filter.sys` and the indices of output and control inputs
can be accessed as the integer vectors
contained in `filter.outputs` and `filter.controls`, respectively.
"""
struct FDFilter{T} <: AbstractFDDObject where T
    sys::DescriptorStateSpace{T}
    outputs::Vector{Int}
    controls::Vector{Int}
    function FDFilter{T}(sys::DescriptorStateSpace{T}, outputs::Vector{Int}, controls::Vector{Int}) where T
        inpout, inpu = fdfilter_validation(sys, outputs, controls)
        mout = length(inpout)
        mu = length(inpu)
        new{T}(sys[:,[inpout; inpu]], Vector(1:mout), mout .+ Vector(1:mu))
    end
end
function fdfilter_validation(sys::DescriptorStateSpace{T}, outputs::Vector{Int}, controls::Vector{Int}) where T
    m = size(sys,2)
    isempty(intersect(controls,outputs)) || error("output and control inputs must be distinct")

    inpout = unique(outputs)
    inpu = unique(controls)

    any(inpout .> m) && error("output inputs indices larger than the number of system inputs $m")
    any(inpu .> m) && error("control inputs indices larger than the number of system inputs $m")

    return inpout, inpu
end
FDFilter(sys::DescriptorStateSpace, outputs::Vector{Int}, controls::Vector{Int})  = FDFilter{eltype(sys)}(sys, outputs, controls)
"""
    FDFilter(sys, p, mu) -> Q::FDFilter

Build for a given linear time-invariant descriptor system model `sys = (A-λE,B,C,D)`,
a fault detection filter object `Q`, as determined with the synthesis functions of FDI filters.
`p` and `mu` are the number of measured outputs and the number of control inputs, respectively.
It is assumed that `B = [By Bu Bv]` and `D = [Dy Du Dv]` are partitioned matrices such that
`By` and `Dy` have `p` columns, and `Bu` and `Du` have `mu` columns,
where `By` and `Bu` are the input matrices from the measured outputs `y` and
control inputs `u`, `Dy` and `Du` are the feedthrough matrices from the measured outputs `y` and
control inputs `u`.

The resulting `Q` contains the partitioned system
`Q.sys = (A-λE,[By Bd],C,[Dy Du])` and the indices of inputs corresponding
to the measured outputs and control inputs are contained in the associated
integer vectors `Q.outputs` and `Q.controls`.
"""
function FDFilter(sys::DescriptorStateSpace, p::Int, mu::Int)
    m = size(sys,2)
    p < 0 && error("number of measured outputs must be non-negative")
    mu < 0 && error("number of control inputs must be non-negative")
    p+mu > m && error("number of measured outputs and control inputs exceeds the number of filter inputs $m")
    return FDFilter{eltype(sys)}(sys[:,1:p+mu], Vector(1:p), Vector(p+1:p+mu))
end
function FDFilter(sys::DescriptorStateSpace{T}; outputs::VRS = Int[], controls::VRS = Int[]) where T
    return FDFilter{T}(sys, vec([outputs; Int[]]), vec([controls; Int[]]))
end
function *(filter1::FDFilter{T1}, filter2::FDFilter{T2}) where {T1,T2}
    return FDFilter(filter1.sys*filter2.sys, length(filter2.outputs), length(filter2.controls))
end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, filter::FDFilter)
    summary(io, filter); println(io)
    display(filter.sys)
    p = length(filter.outputs)
    mu = length(filter.controls)
    if p+mu > 0
       println(io, "Input groups:")
       println(io,          "Name      Channels")
       p > 0 &&  println(io,"outputs   $(filter.outputs)")
       mu > 0 && println(io,"controls  $(filter.controls)")
   end
end
"""
    FDIFilter <: AbstractFDDObject

Type for fault detection and isolation filters resulted as solutions of fault detection and isolation problems.

If `filter::FDIFilter` is the fault detection and isolation filter object,
the underlying `i`-th descriptor system model
can be obtained via `filter.sys[i]` and the indices of output and control inputs
can be accessed as the integer vectors
contained in `filter.outputs` and `filter.controls`, respectively.
"""
struct FDIFilter{T} <: AbstractFDDObject where T
    sys::Vector{DescriptorStateSpace{T}}
    outputs::Vector{Int}
    controls::Vector{Int}
    function FDIFilter{T}(sys::Vector{DescriptorStateSpace{T}}, outputs::Vector{Int}, controls::Vector{Int}) where T
        N = length(sys)
        sysn = similar(sys,N)
        inpout, inpu = fdfilter_validation(sys[1], outputs, controls)
        sysn[1] = sys[1][:,[inpout; inpu]]
        #sys[1] = sys[1][:,[inpout; inpu]]
        for i = 2:N
            inpouti, inpui = fdfilter_validation(sys[i], outputs, controls)
            (inpout == inpouti && inpu == inpui) || error("all component filters must have the same inputs from output and control inputs ")
            sysn[i] = sys[i][:,[inpout; inpu]]
            #sys[i] = sys[i][:,[inpout; inpu]]
        end
        mout = length(inpout)
        mu = length(inpu)
        new{T}(sysn, Vector(1:mout), mout .+ Vector(1:mu))
        #new{T}(sys, Vector(1:mout), mout .+ Vector(1:mu))
    end
end
"""
    FDIFilter(sys, p, mu) -> Q::FDIFilter

Build for a vector of linear time-invariant descriptor system models `sys[i] = (Ai-λEi,Bi,Ci,Di)`
with the same number of inputs, a fault detection and isolation filter object `Q`,
as determined with the synthesis functions of FDI filters.
`p` and `mu` are the number of measured outputs and the number of control inputs, respectively.
It is assumed that each `Bi = [Byi Bui Bvi]` and `Di = [Dyi Dui Dvi]` are partitioned matrices such that
`Byi` and `Dyi` have `p` columns, and `Bui` and `Dui` have `mu` columns,
where `Byi` and `Bui` are the input matrices from the measured outputs `y` and
control inputs `u`, `Dyi` and `Dui` are the feedthrough matrices from the measured outputs `y` and
control inputs `u`.

The resulting `Q` contains the vector of partitioned systems
`Q.sys[i] = (Ai-λEi,[Byi Bdi],Ci,[Dyi Dui])` and the indices of inputs corresponding
to the measured outputs and control inputs are contained in the associated
integer vectors `Q.outputs` and `Q.controls`.
"""
function FDIFilter(sys::Vector{DescriptorStateSpace{T}}, p::Int, mu::Int) where T
    m = size(sys[1],2)
    p < 0 && error("number of measured outputs must be non-negative")
    mu < 0 && error("number of control inputs must be non-negative")
    p+mu > m && error("number of measured outputs and control inputs exceeds the number of filter inputs $m")
    return FDIFilter{T}(sys, Vector(1:p), Vector(p+1:p+mu))
end
function FDIFilter(sys::Vector{DescriptorStateSpace}, p::Int, mu::Int)
    m = size(sys[1],2)
    p < 0 && error("number of measured outputs must be non-negative")
    mu < 0 && error("number of control inputs must be non-negative")
    p+mu > m && error("number of measured outputs and control inputs exceeds the number of filter inputs $m")
    return FDIFilter{eltype(sys[1])}(sys, Vector(1:p), Vector(p+1:p+mu))
end

function FDIFilter(sys::Vector{DescriptorStateSpace{T}}; outputs::VRS = Int[], controls::VRS = Int[]) where T
    return FDIFilter{T}(sys, vec([outputs; Int[]]), vec([controls; Int[]]))
end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, filter::FDIFilter)
    summary(io, filter); println(io)
    for i = 1:length(filter.sys)
        println(io, "Filter number $i:")
        display(filter.sys[i])
    end
    p = length(filter.outputs)
    mu = length(filter.controls)
    if p+mu > 0
       println(io, "Input groups:")
       println(io,          "Name      Channels")
       p > 0 &&  println(io,"outputs   $(filter.outputs)")
       mu > 0 && println(io,"controls  $(filter.controls)")
   end
end

"""
    FDFilterIF <: AbstractFDDObject

Type for the internal form of fault detection filters resulted as solutions of fault detection problems.

If `filter::FDFilterIF` is the fault detection filter internal form object,
the underlying descriptor system model
can be obtained via `filter.sys` and the indices of control, disturbance, fault, noise
and auxiliary inputs can be accessed as the integer vectors
contained in `filter.controls`, `filter.disturbances`, `filter.faults`,
`filter.noise` and `filter.aux`, respectively.
"""
struct FDFilterIF{T} <: AbstractFDDObject where T
    sys::DescriptorStateSpace{T}
    controls::Vector{Int}
    disturbances::Vector{Int}
    faults::Vector{Int}
    noise::Vector{Int}
    aux::Vector{Int}
    function FDFilterIF{T}(sys::DescriptorStateSpace{T}, controls::Vector{Int}, disturbances::Vector{Int},
                           faults::Vector{Int}, noise::Vector{Int}, aux::Vector{Int}) where T
        inpu, inpd, inpf, inpw, inpaux = fdfilterIF_validation(sys, controls, disturbances, faults, noise, aux)
        mu = length(inpu)
        md = length(inpd)
        mf = length(inpf)
        mw = length(inpw)
        maux = length(inpaux)
        new{T}(sys[:,[inpu; inpd; inpf; inpw; inpaux]], Vector(1:mu), mu .+ Vector(1:md) , (mu+md) .+ Vector(1:mf),
               (mu+md+mf) .+ Vector(1:mw), (mu+md+mf+mw) .+ Vector(1:maux))
    end
end
function fdfilterIF_validation(sys::DescriptorStateSpace{T}, controls::Vector{Int}, disturbances::Vector{Int},
                             faults::Vector{Int}, noise::Vector{Int}, aux::Vector{Int}) where T
    m = size(sys,2)
    isempty(intersect(controls,disturbances)) || error("control and disturbance inputs must be distinct")
    isempty(intersect(disturbances,faults)) || error("disturbance and fault inputs must be distinct")
    isempty(intersect(faults,noise)) || error("fault and noise inputs must be distinct")
    isempty(intersect(noise,aux)) || error("noise and aux inputs must be distinct")

    inpu = unique(controls);
    inpd = unique(disturbances);
    inpf = unique(faults);
    inpw = unique(noise);
    inpaux = unique(aux);

    isempty(intersect(inpu,inpd)) || error("control and disturbance inputs must be distinct")
    isempty(intersect(inpd,inpf)) || error("disturbance and fault inputs must be distinct")
    isempty(intersect(inpw,inpf)) || error("noise and fault inputs must be distinct")
    isempty(intersect(inpw,inpaux)) || error("noise and aux inputs must be distinct")

    any(inpu .> m) && error("control inputs indices larger than the number of system inputs $m")
    any(inpd .> m) && error("disturbance inputs indices larger than the number of system inputs $m")
    any(inpf .> m) && error("fault inputs indices larger than the number of system inputs $m")
    any(inpw .> m) && error("noise inputs indices larger than the number of system inputs $m")
    any(inpaux .> m) && error("auxiliary inputs indices larger than the number of system inputs $m")

    return inpu, inpd, inpf, inpw, inpaux
end
FDFilterIF(sys::DescriptorStateSpace, controls::Vector{Int}, disturbances::Vector{Int}, faults::Vector{Int}, noise::Vector{Int}, aux::Vector{Int})  =
           FDFilterIF{eltype(sys)}(sys, controls, disturbances, faults, noise, aux)
"""
    FDFilterIF(sys, mu, md, mf, mw = 0, maux = 0; moff = 0 ) -> R::FDFilterIF

Build for a given linear time-invariant descriptor system model `sys = (A-λE,B,C,D)`,
a fault detection filter internal form object `R`, as determined with the synthesis functions of FDI filters.
`mu`, `md`, `mf`, `mw` and `maux` are the dimensions of control, disturbance, fault, noise and auxiliary input vectors, respectively.
It is assumed that `B = [Boff Bu Bd Bf Bw Bv]` and `D = [Doff Du Dd Df Dw Dv]` are partitioned matrices such that
`Boff` and `Doff` have `moff` columns, `Bu` and `Du` have `mu` columns, `Bd` and `Dd` have `md` columns,
`Bf` and `Df` have `mf` columns,  `Bw` and `Dw` have `mw` columns, and `Bv` and `Dv` have `maux` columns.

The resulting `R` contains the partitioned system
`R.sys = (A-λE,[Bu Bd Bf Bw Bv],C,[Du Dd Df Dw Dv])` and the indices of inputs corresponding
to the control, disturbance, fault, noise and auxiliary inputs are contained in the associated
integer vectors `R.controls`, `R.disturbances`, `R.faults`, `R.noise` and `R.aux`.
"""
function FDFilterIF(sys::DescriptorStateSpace{T}, mu::Int, md::Int, mf::Int, mw::Int = 0, maux::Int = 0; moff::Int = 0) where T
    m = moff+mu+md+mf+mw+maux
    m > size(sys,2) && error("number of selected inputs exceeds the number of system inputs")
    return FDFilterIF{T}(sys, moff .+ Vector(1:mu), (moff+mu) .+ Vector(1:md), (moff+mu+md) .+ Vector(1:mf),
                         (moff+mu+md+mf) .+  Vector(1:mw), (moff+mu+md+mf+mw) .+ Vector(1:maux))
end
function FDFilterIF(sys::DescriptorStateSpace{T}; controls::VRS = Int[], disturbances::VRS = Int[],
                faults::VRS = Int[], noise::VRS = Int[], aux::VRS = Int[]) where T
    return FDFilterIF{T}(sys, vec([controls; Int[]]), vec([disturbances; Int[]]),
                         vec([faults; Int[]]), vec([noise; Int[]]), vec([aux; Int[]]))
end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sysf::Union{FDIModel,FDFilterIF})
    summary(io, sysf); println(io)
    display(sysf.sys)
    mu = length(sysf.controls)
    md = length(sysf.disturbances)
    mf = length(sysf.faults)
    mw = length(sysf.noise)
    maux = length(sysf.aux)
    if mu+md+mf+mw+maux > 0
       println(io, "Input groups:")
       println(io, "Name         Channels")
       mu > 0 && println(io,"controls     $(sysf.controls)")
       md > 0 && println(io,"disturbances $(sysf.disturbances)")
       mf > 0 && println(io,"faults       $(sysf.faults)")
       mw > 0 && println(io,"noise        $(sysf.noise)")
       maux > 0 && println(io,"aux          $(sysf.aux)")
    end
end
"""
    FDIFilterIF <: AbstractFDDObject

Type for the internal form of fault detection and isolation filters resulted as solutions of fault detection and isolation problems.

If `filter::FDIFilterIF` is the fault detection and isolation filter internal form object,
the underlying `i`-th descriptor system model
can be obtained via `filter.sys[i]` and the indices of control, disturbance, fault, noise
and auxiliary inputs can be accessed as the integer vectors
contained in `filter.controls`, `filter.disturbances`, `filter.faults`,
`filter.noise` and `filter.aux`, respectively.
"""
struct FDIFilterIF{T} <: AbstractFDDObject where T
    sys::Vector{DescriptorStateSpace{T}}
    controls::Vector{Int}
    disturbances::Vector{Int}
    faults::Vector{Int}
    noise::Vector{Int}
    aux::Vector{Int}
    function FDIFilterIF{T}(sys::Vector{DescriptorStateSpace{T}}, controls::Vector{Int}, disturbances::Vector{Int},
                           faults::Vector{Int}, noise::Vector{Int}, aux::Vector{Int}) where T
        N = length(sys)
        sysn = similar(sys,N)
        inpu, inpd, inpf, inpw, inpaux = fdfilterIF_validation(sys[1], controls, disturbances, faults, noise, aux)
        sysn[1] = sys[1][:,[inpu; inpd; inpf; inpw; inpaux]]
        #sys[1] = sys[1][:,[inpu; inpd; inpf; inpw; inpaux]]
        for i = 2:N
            inpui, inpdi, inpfi, inpwi, inpauxi = fdfilterIF_validation(sys[i], controls, disturbances, faults, noise, aux)
            (inpu == inpui && inpd == inpdi && inpf == inpfi && inpw == inpwi && inpaux == inpauxi) ||
                 error("all component filters must have the same inputs groups")
            sysn[i] = sys[i][:,[inpu; inpd; inpf; inpw; inpaux]]
            #sysn[i] = sys[i][:,[inpu; inpd; inpf; inpw; inpaux]]
            #sys[i] = sys[i][:,[inpu; inpd; inpf; inpw; inpaux]]
        end
        mu = length(inpu)
        md = length(inpd)
        mf = length(inpf)
        mw = length(inpw)
        maux = length(inpaux)
        new{T}(sysn, Vector(1:mu), mu .+ Vector(1:md) , (mu+md) .+ Vector(1:mf),
               (mu+md+mf) .+ Vector(1:mw), (mu+md+mf+mw) .+ Vector(1:maux))
        # new{T}(sys, Vector(1:mu), mu .+ Vector(1:md) , (mu+md) .+ Vector(1:mf),
        #        (mu+md+mf) .+ Vector(1:mw), (mu+md+mf+mw) .+ Vector(1:maux))
    end
end
function FDIFilterIF(sys::Vector{DescriptorStateSpace{T}}; controls::VRS = Int[], disturbances::VRS = Int[],
                     faults::VRS = Int[], noise::VRS = Int[], aux::VRS = Int[]) where T
    return FDIFilterIF{T}(sys, vec([controls; Int[]]), vec([disturbances; Int[]]),
                          vec([faults; Int[]]), vec([noise; Int[]]), vec([aux; Int[]]))
end
FDIFilterIF(sys::Vector{DescriptorStateSpace{T}}, controls::Vector{Int}, disturbances::Vector{Int}, faults::Vector{Int}, noise::Vector{Int}, aux::Vector{Int}) where T  =
           FDIFilterIF{T}(sys, controls, disturbances, faults, noise, aux)
"""
    FDIFilterIF(sys, mu, md, mf, mw = 0, maux = 0; moff = 0 ) -> R::FDIFilterIF

Build for a a vector of linear time-invariant descriptor system models `sys[i] = (Ai-λEi,Bi,Ci,Di)`
with the same number of inputs, ,
a fault detection and isolation filter internal form object `R`, as determined with the synthesis functions of FDI filters.
`mu`, `md`, `mf`, `mw` and `maux` are the dimensions of control, disturbance, fault, noise and auxiliary input vectors, respectively.
It is assumed that each `Bi = [Boffi Bui Bdi Bfi Bwi Bvi]` and `Di = [Doffi Dui Ddi Dfi Dwi Dvi]` are partitioned matrices such that
`Boffi` and `Doffi` have `moff` columns, `Bui` and `Dui` have `mu` columns, `Bdi` and `Ddi` have `md` columns,
`Bfi` and `Dfi` have `mf` columns,  `Bwi` and `Dwi` have `mw` columns, and `Bvi` and `Dvi` have `maux` columns.

The resulting `R` contains the vector of partitioned systems
`R.sys[i] = (A-λE,[Bui Bdi Bfi Bwi Bvi],C,[Dui Ddi Dfi Dwi Dvi])` and the indices of inputs corresponding
to the control, disturbance, fault, noise and auxiliary inputs are contained in the associated
integer vectors `R.controls`, `R.disturbances`, `R.faults`, `R.noise` and `R.aux`.
"""
function FDIFilterIF(sys::Vector{DescriptorStateSpace{T}}, mu::Int, md::Int, mf::Int, mw::Int = 0, maux::Int = 0; moff::Int = 0) where T
    m = moff+mu+md+mf+mw+maux
    m > size(sys[1],2) && error("number of selected inputs exceeds the number of system inputs")
    return FDIFilterIF{T}(sys, moff .+ Vector(1:mu), (moff+mu) .+ Vector(1:md), (moff+mu+md) .+ Vector(1:mf),
                         (moff+mu+md+mf) .+  Vector(1:mw), (moff+mu+md+mf+mw) .+ Vector(1:maux))
end

function FDIFilterIF(sys::Vector{DescriptorStateSpace}, mu::Int, md::Int, mf::Int, mw::Int = 0, maux::Int = 0; moff::Int = 0)
    m = moff+mu+md+mf+mw+maux
    m > size(sys[1],2) && error("number of selected inputs exceeds the number of system inputs")
    return FDIFilterIF{eltype(sys[1])}(sys, moff .+ Vector(1:mu), (moff+mu) .+ Vector(1:md), (moff+mu+md) .+ Vector(1:mf),
                         (moff+mu+md+mf) .+  Vector(1:mw), (moff+mu+md+mf+mw) .+ Vector(1:maux))
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, filter::FDIFilterIF)
    summary(io, filter); println(io)
    for i = 1:length(filter.sys)
        println(io, "Filter number $i:")
        display(filter.sys[i])
    end
    mu = length(filter.controls)
    md = length(filter.disturbances)
    mf = length(filter.faults)
    mw = length(filter.noise)
    maux = length(filter.aux)
    if mu+md+mf+mw+maux > 0
       println(io, "Input groups:")
       println(io, "Name         Channels")
       mu > 0 && println(io,"controls     $(filter.controls)")
       md > 0 && println(io,"disturbances $(filter.disturbances)")
       mf > 0 && println(io,"faults       $(filter.faults)")
       mw > 0 && println(io,"noise        $(filter.noise)")
       maux > 0 && println(io,"aux          $(filter.aux)")
    end
end

"""
    fdimodset(sys; controls, c, disturbances, d, faults, f, fa, faults_sen, fs, noise, n, aux) -> sysf::FDIModel

Build for a given linear time-invariant system model `sys = (A-λE,B,C,D)`, a synthesis model object `sysf::FDIModel`
to be used in conjunction with the analysis and synthesis functions of FDI filters.
If `sys` is a vector of system models, then `sysf` results as a vector of synthesis model objects.

The information on the partition of the input components in control, disturbance, fault, noise and auxiliary
inputs can be specified using the following keyword arguments:

`controls = inpu` or `c = inpu` specifies the indices `inpu` of the control inputs (default: void)

`disturbances = inpd` or `d = inpd` specifies the indices `inpd` of the disturbance inputs (default: void)

`faults = inpf` or `f  = inpf` specifies the indices `inpf` of the fault inputs (default: void)

`fa = inpfa ` specifies the indices `inpfa` of control inputs subject to actuator fault (default: void)

`faults_sen = inpfs` or `fs = inpfs` specifies the indices `inpfs` of the system outputs subject to sensor fault inputs (default: void)

`noise = inpn` or  `noise = inpn` specifies the indices `inpn` of the noise inputs (default: void)

`aux = inpa` specifies the indices `inpa` of the auxiliary inputs (default: void)

The indices of inputs or outputs can be specified as integer vectors, integer scalars or integer `UnitRange`s.

The resulting `sysf` contains the partitioned system
`sysf.sys = (A-λE,[Bu Bd Bf Bw Bv],C,[Du Dd Df Dw Dv])`, where
`Bu`, `Bd`, `Bf`, `Bw` and `Bv` are the input matrices from the control inputs `u`, disturbance inputs `d`, fault inputs `f`,
noise inputs `w` and auxiliary inputs `v`, respectively, and `Du`, `Dd`, `Df`, `Dw` and `Dv` are the feedthrough matrices from those inputs.
The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors
`sysf.controls`, `sysf.disturbances`, `sysf.faults`, `sysf.noise` and `sysf.aux`, respectively.

_Method:_ If `G(λ)` is the `p x m` transfer function matrix of `sys`, then the resulting system `sysf` has an
equivalent input output form `[Gu(λ) Gd(λ) Gf(λ) Gw(λ) Gv(λ)]`, where the following relations define the component matrices:
`Gu(λ) = G(λ)*Su`, `Gd(λ) = G(λ)*Sd`,  `Gf(λ) = [G(λ)*Sf Ss]`, `Gw(λ) = G(λ)*Sw`,  `Gv(λ) = G(λ)*Sv`, with
the selection matrices `Su`, `Sd`, `Sf`, `Sw` and `Sv` formed from the columns of the `m`-th order identity matrix and
`Ss` is formed  from the columns of the `p`-th order identity matrix.
"""
function fdimodset(sys::DescriptorStateSpace{T};
                   controls::VRS = Int[], c::VRS = Int[],
                   disturbances::VRS = Int[], d::VRS = Int[],
                   faults::VRS = Int[], f::VRS = Int[], fa::VRS = Int[],
                   faults_sen::VRS = Int[], fs::VRS = Int[],
                   noise::VRS = Int[], n::VRS = Int[],
                   aux::VRS = Int[]) where T
    p, m = size(sys)
    nx = order(sys)
    inpu = unique([controls;c])
    inpd = unique([disturbances;d])
    isempty(intersect(inpu,inpd)) || error("control and disturbance inputs must be distinct")
    inpf1 = unique([faults;f;fa])
    isempty(intersect(inpd,inpf1)) || error("disturbance and fault inputs must be distinct")
    inpfs = unique([faults_sen;fs])
    inpw = unique([noise;n])
    isempty(intersect(inpw,inpf1)) || error("noise and fault inputs must be distinct")
    inpaux = unique([aux;Int[]])
    mu = length(inpu)
    md = length(inpd)
    mf1 = length(inpf1)
    mf2 = length(inpfs)
    mf = mf1+mf2
    mw = length(inpw)
    maux = length(inpaux)
    mu == 0 || maximum(inpu) <= m || error("selected index/indices of control inputs larger than the number of system inputs $m")
    md == 0 || maximum(inpd) <= m || error("selected index/indices of disturbance inputs larger than the number of system inputs $m")
    mf1 == 0 || maximum(inpf1) <= m || error("selected index/indices of fault inputs larger than the number of system inputs $m")
    mf2 == 0 || maximum(inpfs) <= p || error("selected index/indices of sensor fault inputs larger than the number of system outputs $p")
    mw == 0 || maximum(inpw) <= m || error("selected index/indices of noise inputs larger than the number of system inputs $m")
    maux == 0 || maximum(inpaux) <= m || error("selected index/indices of auxiliary inputs larger than the number of system inputs $m")
    Dsf = mf2 > 0 ? eye(T, p, p)[:,inpfs] : zeros(T, p, 0)

    Be = [sys.B[:,inpu] sys.B[:,inpd] sys.B[:,inpf1] zeros(T, nx, mf2) sys.B[:,inpw] sys.B[:,inpaux]];
    De = [sys.D[:,inpu] sys.D[:,inpd] sys.D[:,inpf1] Dsf sys.D[:,inpw] sys.D[:,inpaux]];
    return FDIModel{T}(dss(sys.A, sys.E, Be, sys.C, De, Ts = sys.Ts),
            controls = Vector(1:mu), disturbances = mu .+ Vector(1:md) , faults = (mu+md) .+ Vector(1:mf),
            noise = (mu+md+mf) .+ Vector(1:mw), aux = (mu+md+mf+mw) .+ Vector(1:maux))
end
function fdimodset(sys::Vector{DescriptorStateSpace}; kwargs...)
    N = length(sys)
    sysf = similar(Vector{FDIModel},N)
    for i = 1:N
        sysf[i] = fdimodset(sys[i]; kwargs...)
    end
    return sysf
end
"""
    fdIFeval(sysQ::FDFilter, sysf::FDIModel; minimal = false, atol, atol1 = atol, atol2 = atol, rtol, fast = true) -> sysR::FDFilterIF

Compute the internal form `sysR` of the fault detection filter `sysQ` applied to the synthesis model `sysf`.
If `sysf` has the partitioned transfer function matrix `G(λ) = [ Gu(λ)  Gd(λ) Gf(λ) Gw(λ) Gv(λ) ]` in accordance with
the partitioned system inputs as `controls`, `disturbances`, `faults`, `noise` and `auxiliary` inputs, respectively,
and `Q(λ) = [ Qy(λ) Qu(λ) ]` is the partitioned transfer function matrix of the fault detection filter `sysQ`
in accordance with the partitioned filter inputs as `outputs` and `controls`, then
the transfer function matrix `R(λ)` of the resulting internal form `sysR` is given by

     R(λ) = | Qy(λ)  Qu(λ) | * | Gu(λ)  Gd(λ) Gf(λ) Gw(λ) Gv(λ) |
                               |  I       0     0     0     0   |

A minimal descriptor realization is computed if `minimal = true` and a possibly non-minimal
realization is determined if `minimal = false` (default).

The minimal realization computation relies on pencil manipulation algorithms which
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

If `(Ar-λEr,Br,Cr,Dr)` is the full order descriptor realization of `sysR.sys`, then
the keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `Ar`, `Br`, `Cr`, `Dr`, the absolute tolerance for the nonzero elements of `Er`,
and the relative tolerance for the nonzero elements of `Ar`, `Br`, `Cr`, `Dr` and `Er`.
The default relative tolerance is `nϵ`, where `ϵ` is the working _machine epsilon_
and `n` is the order of the system `sysR`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

"""
function fdIFeval(Q::FDFilter, sysf::FDIModel; minimal::Bool = false,
                  atol::Real = 0., atol1::Real = atol, atol2::Real = atol, rtol::Real = 0., fast::Bool = true)
   # compute the internal form as
   #   R = Q * [ Gu Gd Gf Gw Gaux]
   #           [ I  0  0  0  0   ]
   p, m = size(sysf.sys)
   mu = length(sysf.controls)
   (length(Q.outputs) == p && length(Q.controls) == mu) || error("filter Q is incompatible with the given system")
   # return fdRset(gminreal(Q.sys*[sysf.sys;  I zeros(eltype(sysf.sys),mu,m-mu)]; atol1, atol2, rtol);
   #               sysf.controls, sysf.disturbances, sysf.faults, sysf.noise, sysf.aux)
   sysR = Q.sys*[sysf.sys;  I zeros(eltype(sysf.sys),mu,m-mu)]
   return FDFilterIF(minimal ? gminreal(sysR; atol1, atol2, rtol, fast) : sysR,
                 sysf.controls, sysf.disturbances, sysf.faults, sysf.noise, sysf.aux)

end
"""
    fdIFeval(sysQ::FDIFilter, sysf::FDIModel; minimal = false, atol, atol1 = atol, atol2 = atol, rtol, fast = true) -> sysR::FDIFilterIF

Compute the internal form `sysR` of the fault detection and isolation filter `sysQ` applied to the synthesis model `sysf`.
If `sysf` has the partitioned transfer function matrix `G(λ) = [ Gu(λ)  Gd(λ) Gf(λ) Gw(λ) Gv(λ) ]` in accordance with
the partitioned system inputs as `controls`, `disturbances`, `faults`, `noise` and `auxiliary` inputs, respectively,
and `Qi(λ) = [ Qyi(λ) Qui(λ) ]` is the partitioned transfer function matrix of the `i`-th filter `sysQ.sys[i]`
in accordance with the partitioned filter inputs as `outputs` and `controls`, then
the transfer function matrix `Ri(λ)` of the `i`-th filter in the resulting internal form `sysR.sys[i]` is given by

     Ri(λ) = | Qyi(λ)  Qui(λ) | * | Gu(λ)  Gd(λ) Gf(λ) Gw(λ) Gv(λ) |
                                  |  I       0     0     0     0   |

Minimal descriptor realizations are computed if `minimal = true` and a possibly non-minimal
realization is determined if `minimal = false` (default).

The minimal realization computation relies on pencil manipulation algorithms which
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

If `(Ari-λEri,Bri,Cri,Dri)` is the full order descriptor realization of `sysR.sys[i]`, then
the keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `Ari`, `Bri`, `Cri`, `Dri`, the absolute tolerance for the nonzero elements of `Eri`,
and the relative tolerance for the nonzero elements of `Ari`, `Bri`, `Cri`, `Dri` and `Eir`.
The default relative tolerance is `ni*ϵ`, where `ϵ` is the working _machine epsilon_
and `ni` is the order of the system `sysR.sys[i]`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

"""
function fdIFeval(Q::FDIFilter, sysf::FDIModel; minimal::Bool = false,
                  atol::Real = 0., atol1::Real = atol, atol2::Real = atol, rtol::Real = 0., fast::Bool = true)
   # compute the internal form as
   #   R[i] = Q[i] * [ Gu Gd Gf Gw Gaux]
   #                 [ I  0  0  0  0   ]
   N = length(Q.sys)
   sysR = similar(Q.sys,N)
   p, m = size(sysf.sys)
   mu = length(sysf.controls)
   (length(Q.outputs) == p && length(Q.controls) == mu) || error("filter Q is incompatible with the given system")
   for i = 1:N
       sysR[i] = minimal ? gminreal(Q.sys[i]*[sysf.sys;  I zeros(eltype(sysf.sys),mu,m-mu)]; atol1, atol2, rtol, fast) :
                           Q.sys[i]*[sysf.sys;  I zeros(eltype(sysf.sys),mu,m-mu)]
   end
   return FDIFilterIF{eltype(Q.sys[1])}(sysR, sysf.controls, sysf.disturbances, sysf.faults, sysf.noise, sysf.aux)
end
gbalmr(Q::FDFilter{T}; kwargs...) where T = FDFilter(gbalmr(Q.sys; kwargs...)[1], Q.outputs, Q.controls)
gminreal(Q::FDFilter{T}; kwargs...) where T = FDFilter(gminreal(Q.sys; kwargs...), Q.outputs, Q.controls)
function gbalmr(Q::FDIFilter{T}; kwargs...) where T
    for i = 1:length(Q.sys)
        Q.sys[i] = gbalmr(Q.sys[i]; kwargs...)[1]
    end
    return Q
end
function gminreal(Q::FDIFilter{T}; kwargs...) where T
    for i = 1:length(Q.sys)
        Q.sys[i] = gminreal(Q.sys[i]; kwargs...)
    end
    return Q
end
gpole(Q::FDFilter{T}; kwargs...) where T = gpole(Q.sys;  kwargs...)
gpole(Q::FDIFilter{T}; kwargs...) where T = gpole.(Q.sys;  kwargs...)

# left multiplication of FD objects with descriptor system object
function *(sys::DescriptorStateSpace{T1}, sysr::FDFilterIF{T2}) where {T1,T2}
    return FDFilterIF(sys*sysr.sys, sysr.controls, sysr.disturbances, sysr.faults, sysr.noise, sysr.aux)
end
function *(sys::DescriptorStateSpace{T1}, sysr::FDIModel{T2}) where {T1,T2}
    return FDIModel(sys*sysr.sys, sysr.controls, sysr.disturbances, sysr.faults, sysr.noise, sysr.aux)
end
function *(sys::Vector{DescriptorStateSpace{T1}}, sysr::FDIFilterIF{T2}) where {T1,T2}
    return FDIFilterIF(sys .* sysr.sys, sysr.controls, sysr.disturbances, sysr.faults, sysr.noise, sysr.aux)
end

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
using EchoStateNetworks
using LinearAlgebra
using CSV
using Random

T = Float64
filepath = joinpath(dirname(@__FILE__),"MackeyGlass_t17.txt")
# if !isfile(filepath)
#     println("download data...")
#     run(`curl -o $(filepath) http://minds.jacobs-university.de/sites/default/files/uploads/mantas/code/MackeyGlass_t17.txt`)
#     println("finish.")
# end
data = convert(Array{T},CSV.read(filepath))'

train_range = 1:2000
test_range = 2001:4000
data_range = 1:4000

esn = EchoStateNetwork{T}(Ni=1, No=1, Nr=500,
                          spectral_radius=1.2,
                          sparsity=0.95,
                          leaking_rate=0.99,
                          teacher_forcing=true,
                          rng=MersenneTwister(0x59))

y_train = data[:,train_range]
x_train = ones(T, 1,length(train_range))
y_test = data[:,test_range]
x_test = ones(T, 1,length(test_range))

y_pred = train!(esn, x_train, y_train, reg=1.0e-8, discard=100)
println("train RMSE: $(norm(y_train-y_pred)/sqrt(length(train_range)))")
y_pred = predict!(esn, x_test)
println("test RMSE:  $(norm(y_test-y_pred)/sqrt(length(test_range)))")


text:
using FastGaussQuadrature

function nearlyequal{T,U}(x::T,y::T,τ::U)
    X = norm(x)
    Y = norm(y)
    D = norm(x-y)

    x == y && return true

    (X <= τ && Y <= τ) && return true
    2 * D / (X+Y) < τ
end

function legendre(n,a,b)
    x,w = FastGaussQuadrature.gausslegendre(n)
    w *= (b-a)/2
    x = a+(x+1)*(b-a)/2
    return x, w
end

function dblquadints1{N}(v1,v2,v3,x,::Type{Val{N}},ri=-1.0,ro=1.0e15)
    G = 10000
    s, w = legendre(G, 0.0, 1.0)
    t1 = v1-v3
    t2 = v2-v3
    n = cross(t1,t2)
    I = zeros(eltype(x),N+3)
    K = zeros(typeof(x),N+3)
    a2 = norm(n)
    a2 < eps(eltype(x)) && return I, K
    n = normalize(n)
    d = dot(x-v1,n)
    ξ = x - d*n
    for g in 1:G
        u = s[g]
        for h in 1:G
            v = (1-u)*s[h]
            y = v3 + u*t1 + v*t2
            R = norm(x-y)
            (ri <= R <= ro) || continue
            j = a2 * w[g] * w[h] * (1-u)
            I[1] += j * d / R^3
            K[1] += j * (y-ξ) / R^3
            for n = -1 : N
                I[n+3] += j * R^n
                K[n+3] += j * (y-ξ) * R^n
            end
        end
    end
    return I, K
end


text:
module Tuxedo

export @tux

include("tux.jl")
include("src.jl")
include("test.jl")

end

max_stars_repo_name: milktrader/Tuxedo

max_issues_repo_name: milktrader/Tuxedo

max_forks_repo_name: milktrader/Tuxedo


text:
using Pkg
Pkg.update()
ENV["PYTHON_JL_RUNTIME_PYTHON"] = Sys.which("python")
Pkg.build("PyCall")
module HHNSGA
    using Pkg
    using Plots
    using UnicodePlots
    using OrderedCollections
    using LinearAlgebra
    using UnicodePlots
    using PyCall
    using NSGAII
    export NSGAII

    export SNN
    Pkg.add(url="https://github.com/Datseris/SpikeSynchrony.jl")
    using SpikeSynchrony
    #import Pkg;
    Pkg.add("OrderedCollections")
    Pkg.add("NSGAII")#` to install the NSGAII package.

    include("../src/SpikingNeuralNetworks.jl")
    include("../src/units.jl")
    include("../src/plot.jl")

    SNN = SpikingNeuralNetworks.SNN
    using Random
    using PyCall
    py"""
    import copy
    import matplotlib
    import matplotlib.pyplot as plt
    import numpy as np
    matplotlib.get_backend()
    import sys
    import os
    sys.path.append(os.getcwd())
    """
    Main.include("../src/SpikingNeuralNetworks.jl")
    Main.eval("SNN = SpikingNeuralNetworks.SNN")

    """
    ##model = py"IZModel()"
    simple = py"IZModel()"

    py"""
    import pickle
    import numpy as np
    import random
    from sciunit.scores.collections import ScoreArray
    from sciunit import TestSuite
    import sciunit
    import random

    from neuronunit.tests.fi import RheobaseTestP
    from collections import OrderedDict
    #from neuronunit import get_neab
    #import pdb
    #pdb.set_trace()
    cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
    for test in cell_tests.values():
        if "Rheobase test" in test.keys():
            temp_test = {k:v for k,v in test.items()}
            break
    rt = temp_test["Rheobase test"]
    #rtp = RheobaseTestP(rt.observation)
    JUIZI = {
        'a': 0.02,
        'b': 0.2,
        'c': -65,
        'd': 8,
    }

    #ranges = OrderedDict({k:[v-0.5*np.abs(v),v+0.5*np.abs(v)] for k,v in copy.copy(JUIZI).items()})
    temp = {'a':[0.02,0.1],'b':[0.2,0.26],'c':[-65,-50],'d':[0.05,8]}
    ranges = OrderedDict(temp)

    N = len(JUIZI)
    """
    N = py"N"
    ranges = py"ranges"
    #H1=[values(ranges)]
    #current_params = py"rt.params"
    #simple.attrs = py"JUIZI"
    py"""
    from izhi import IZModel

    def evaluate(test,god):
        model = IZModel(attrs=god)
        model.attrs = god
        rt = test['Rheobase test']

        rheo = rt.generate_prediction(model)

        scores_ = []
        test = {k:v for k,v in test.items() if v.observation is not None}
        for temp_test in test.values():
            if 'injected_square_current' not in temp_test.params.keys():
               temp_test.params['injected_square_current'] = {}
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "InjectedCurrentAPWidthTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "Injected CurrentAPThresholdTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "InjectedCurrentAPAmplitudeTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "RheobaseTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

        to_map = [(tt,model) for tt in test.values() ]
        def map_score(content):
            temp_test = content[0]
            model = content[1]
            temp_test.score_type = sciunit.scores.ZScore
            try:
                score = temp_test.judge(model)
                score = np.abs(score.log_norm_score)
            except:
                try:
                    score = np.abs(float(score.raw_score))
                except:
                    score = 100
            if isinstance(score, sciunit.scores.incomplete.InsufficientDataScore):
                score = 100
            return score

        scores_ = list(map(map_score,to_map))
        tt = TestSuite(list(test.values()))
        SA = ScoreArray(tt, scores_)
        errors = tuple(SA.values,)
        return errors
    """

    py"""
    def z_py(evaluate,god):

        cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
        tests = cell_tests['Neocortex pyramidal cell layer 5-6']

        model_type="RAW"
        #from neuronunit.optimisation import make_sim_tests
        fps = ['a','b','c','d']
        #sim_tests, OM, target = make_sim_tests.test_all_objective_test(fps,model_type=model_type)
        SA = evaluate(tests,god)
        return SA
    """

    export transdict
    function transdict(x)
        bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

        py"""
        def map_dict(x):
            genes_out = x
            genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
            return genes_out, genes_out_dic
        """
        decoded = x

        try
           decoded = NSGAII.decode(x,bc)
        catch
           decoded = x
        end
        _,god_dict = py"map_dict"(decoded)
        #bincoded = NSGAII.encode(genes_out, bc)
        return god_dict#, genes_out_dic
    end
    py"""
    lower_list = [v[0] for k,v in ranges.items()]
    upper_list = [v[1] for k,v in ranges.items()]
    """

    export init_function
    function init_function()
        #const
        bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

        py"""
        def pheno_map():

            lower_list = [v[0] for k,v in ranges.items()]
            upper_list = [v[1] for k,v in ranges.items()]
            gene_out = [random.uniform(lower, upper) for lower, upper in zip(lower_list, upper_list)]
            #gene_out_dic = OrderedDict({k:gene_out[i] for i,k in enumerate(ranges.keys()) })
            return gene_out#, gene_out_dic
        """
        gene_out = py"pheno_map"()
        bincoded = NSGAII.encode(gene_out, bc)
        return bincoded
    end

    function references_for_H()
        bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

        py"""
        from sklearn.model_selection import ParameterGrid
        grid = ParameterGrid(ranges)
        genes = []
        for g in grid:
            genes.append(list(g.values()))
        genes_out_dic = OrderedDict({k:genes[i] for i,k in enumerate(ranges.keys()) })

        """
        genes_out = py"genes"
        bcd = []
        for g in genes_out
             bincoded = NSGAII.encode(g, bc)
             append!(bcd,bincoded)
        end
        return bcd
    end

    export z
    function z(x)
        god = transdict(x)
        contents = py"z_py"(py"evaluate",god)
        return contents
    end

end
#create_sysimage(:HHNSGA, sysimage_path="sys_HHNSGA.so", precompile_execution_file="precompile_HHNSGA.jl")


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using UniqueVectors: UniqueVector
using AdmittanceModels: lossless_modes_dense, apply_transform

geo_code = GeoCode(gdspy_open(joinpath(@__DIR__, "boxes.gds"), "main"),
    layer_translation=Dict(12 => -50, 14 => -20),
    layer_extrusion=Dict(12 => 100, 14 => 40),
    layer_volume_difference=[(12, 14, false)])
file_name = joinpath(@__DIR__, "boxes.geo")
geo_write!(file_name, geo_code, characteristic_length_factor=.9)

initialize!()
gmsh_open!(file_name)
N, K = 3, 4
mesh!(K)
node_tags, points, tcomp = get_triangulated_complex(N, K)
group_dict = get_physical_groups(node_tags, points)
comp = tcomp.complex
orient!(comp)
m = Metric(N)
mesh = Mesh(tcomp, circumcenter(m))

μ⁻, Λ⁻, σ, ϵ = 1/DiscretePDEs.μ₀, 0, 0, DiscretePDEs.ϵ₀
μ⁻_form = get_material(comp, μ⁻, 3)
Λ⁻_form = get_material(comp, Λ⁻, 2)
σ_form = get_material(comp, σ, 2)
ϵ_form = get_material(comp, ϵ, 2)

pso, null_basis = electrodynamics_pso(m, mesh, Vector{Cell{N}}[], boundary(comp),
    μ⁻_form, Λ⁻_form, σ_form, ϵ_form)
constrained_pso = apply_transform(pso, null_basis)
λs, vs = lossless_modes_dense(constrained_pso, min_freq=1e6)
freq = imag.(λs[1])/(2π)

hfss_answer = 1.97194e6
@testset "modes of box with superconductors" begin
    @test isapprox(freq, hfss_answer, rtol=6e-2)
end

if false # set to true to display plots with Gmsh
    v = null_basis * vs[:,1]
    v /= maximum(abs.(v))

    vec_A = sharp(m, comp, v)
    comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
    ordering = [findfirst(isequal(p), comp_points) for p in points]
    add_field!("Vector potential", node_tags, vec_A[ordering])
    gui!()
end
rm(file_name)


text:
using BinDeps

@BinDeps.setup

using Compat

group = library_group("groonga")

deps = [
    groonga = library_dependency("groonga", aliases = ["libgroonga0", "libgroonga","libgroonga.so.0"], depends = [], group = group)
]

provides(AptGet,
    @compat Dict(
        "libgroonga0" => groonga
    ))

const version = "5.0.4"

provides(Sources,
    @compat Dict(
         URI("http://packages.groonga.org/source/groonga/groonga-$(version).tar.gz") => groonga
    ))

prefix = joinpath(BinDeps.depsdir(groonga), "usr")
srcdir = joinpath(BinDeps.depsdir(groonga), "src", "groonga-$(version)")

provides(SimpleBuild,
    (@build_steps begin
        GetSources(groonga)
        @build_steps begin
            ChangeDirectory(srcdir)
            `./configure --prefix=$prefix`
            `make`
            `make install`
        end
     end), groonga, os = :Unix)

@BinDeps.install [:groonga => :libgroonga]


text:
# ---
# title: 1448. Count Good Nodes in Binary Tree
# id: problem1448
# author: Jimmy Shen
# date: 2020-08-18
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/count-good-nodes-in-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree `root`, a node _X_ in the tree is named  **good** if in
# the path from root to _X_ there are no nodes with a value _greater than_ X.
#
# Return the number of **good** nodes in the binary tree.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/02/test_sample_1.png)**
#
#
#
#     Input: root = [3,1,4,3,null,1,5]
#     Output: 4
#     Explanation: Nodes in blue are **good**.
#     Root Node (3) is always a good node.
#     Node 4 -> (3,4) is the maximum value in the path starting from the root.
#     Node 5 -> (3,4,5) is the maximum value in the path
#     Node 3 -> (3,1,3) is the maximum value in the path.
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/02/test_sample_2.png)**
#
#
#
#     Input: root = [3,3,null,4,2]
#     Output: 3
#     Explanation: Node 2 -> (3, 3, 2) is not good, because "3" is higher than it.
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: 1
#     Explanation: Root is considered as **good**.
#
#
#
# **Constraints:**
#
#   * The number of nodes in the binary tree is in the range `[1, 10^5]`.
#   * Each node's value is between `[-10^4, 10^4]`.
#
#
## @lc code=start
using LeetCode

function good_nodes(root::TreeNode)::Int
    dfs(::Nothing, args...) = 0
    function dfs(n::TreeNode, mval::Int)::Int
        res = 0
        (n.val >= mval) && (res += 1)
        n_max = max(mval, n.val)
        res += dfs(n.left, n_max)
        res += dfs(n.right, n_max)
        return res
    end
    return dfs(root, typemin(Int))
end
## @lc code=end