text:
"""
AddColPetro(args...)  -->DataFrames

It adds specified colums to a given DataFrames

<br> The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Log|Mandatory|--|DataFrame|This is the entire WellLog|
|DenPhi|Opcional|false|DenPhi=Bool|Adds a Porosity from density column to Log DataFramed called DenPhi|
|Phie|Opcional|false|Phie=Bool|Adds a Effective porosity column to Log DataFramed called Phie|
|Sw|Opcional|false|Sw=Bool|Adds a Water Saturarion column to Log DataFramed called Sw|
|Perm|Opcional|false|Perm=Bool|Adds a Permeability column to Log DataFramed called Perm|
|PayFlag|Opcional|false|PayFlag=Bool|Adds a Pay column to Log DataFramed called PayFlag|
|Kh|Opcional|false|Kh=Bool|Adds a Flow Capacity column to Log DataFramed called Sw|
|All|Opcional|false|All=Bool|Adds all of column listed above|
"""
function AddColPetro(Log;Vsh=false,
                     DenPhi=false,
                     Phie=false,
                     Sw=false,
                     Perm=false,
                     PayFlag=false,
                     Kh=false,
                     All=false)

    if All==true
       Log.Vsh=zeros(size(Log,1)).+1
       Log.Phie=zeros(size(Log,1))
       Log.Sw=zeros(size(Log,1))
       Log.Perm=zeros(size(Log,1))
       Log.PayFlag=zeros(size(Log,1))
       Log.Kh=zeros(size(Log,1))
       Log.KhCum=zeros(size(Log,1))
       Log.DenPhi=zeros(size(Log,1))
    else

        if Vsh==true
        Log.Vsh=zeros(size(Log,1)).+1
        end

        if Phie==true
        Log.Phie=zeros(size(Log,1))
        end

        if Sw==true
        Log.Sw=zeros(size(Log,1))
        end

        if Perm==true
        Log.Perm=zeros(size(Log,1))
        end

        if PayFlag==true
        Log.PayFlag=zeros(size(Log,1))
        end

        if Kh==true
        Log.Kh=zeros(size(Log,1))
        Log.KhCum=zeros(size(Log,1))
        end

        if DenPhi==true
        Log.DenPhi=zeros(size(Log,1))
        end
    end

    return Log
end


text:
using Documenter, Example, clima_demo
makedocs(
  sitename="clima_demo",
  doctest = false,
  strict = false,
  format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",),
  modules = [Documenter, clima_demo],
  clean = false,
  pages = Any[
    "Home" => "index.md",
  ]
)
deploydocs(
           repo = "github.com/costachris/clima_demo.git",
           target = "build",
          )

ext: jl


text:
using InvertibleNetworks, LinearAlgebra, Test, Flux

import InvertibleNetworks: reset!
import ChainRulesCore: rrule

# Initialize invertible/non-invertible layers
nx = 32
ny = 32
n_ch = 16
n_hidden = 64
batchsize = 2
logdet = false
N1 = CouplingLayerHINT(n_ch, n_hidden; logdet=logdet, permute="full")
N2 = CouplingLayerHINT(n_ch, n_hidden; logdet=logdet, permute="full")
N3 = CouplingLayerHINT(n_ch, n_hidden; logdet=logdet, permute="full")
N4 = Chain(Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1), x -> relu.(x), Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1))
N5 = Chain(Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1), x -> relu.(x), Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1))
N6 = CouplingLayerHINT(n_ch, n_hidden; logdet=logdet, permute="full")
N7 = CouplingLayerHINT(n_ch, n_hidden; logdet=logdet, permute="full")
N8 = Chain(Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1), x -> relu.(x), Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1))
N9 = Chain(Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1), x -> relu.(x), Conv((3, 3), n_ch => n_ch; stride = 1, pad = 1))
N10 = CouplingLayerHINT(n_ch, n_hidden; logdet=logdet, permute="full")

# Forward pass + gathering pullbacks
function fw(X)
    X1, ∂1 = rrule(N1, X)
    X2, ∂2 = rrule(N2, X1)
    X3, ∂3 = rrule(N3, X2)
    X5, ∂5 = Flux.Zygote.pullback(Chain(N4, N5), X3)
    X6, ∂6 = rrule(N6, X5)
    X7, ∂7 = rrule(N7, X6)
    X9, ∂9 = Flux.Zygote.pullback(Chain(N8, N9), X7)
    X10, ∂10 = rrule(N10, X9)
    d1 = x -> ∂1(x)[2]
    d2 = x -> ∂2(x)[2]
    d3 = x -> ∂3(x)[2]
    d5 = x -> ∂5(x)[1]
    d6 = x -> ∂6(x)[2]
    d7 = x -> ∂7(x)[2]
    d9 = x -> ∂9(x)[1]
    d10 = x -> ∂10(x)[2]
    return X10, ΔY -> d1(d2(d3(d5(d6(d7(d9(d10(ΔY))))))))
end

# Gradient Test
X = randn(Float32, nx, ny, n_ch, batchsize)
Y0 = randn(Float32, nx, ny, n_ch, batchsize)

loss(X) = 0.5f0*norm(N(X) - Y0)^2

Y, ∂Y = fw(X)

g = ∂Y(Y-Y0)

ΔX_ = randn(Float32, nx, ny, n_ch, batchsize)

# # Integrated backward pass w/ Zygote
reset!(GLOBAL_STATE_INVOPS)
N = Chain(N1, N2, N3, N4, N5, N6, N7, N8, N9, N10);

g2 = gradient(X -> loss(X), X)

@test g ≈ g2[1] rtol=1f-6


text:
using HDF5

function save_results(file_name, key, data, overwrite=false)
    mode = overwrite ? "w" : "cw"
    h5open(file_name, mode) do file
        id = update_current_count(file, key)
        file["clients/$key/$id"] = data
    end
end

function update_current_count(file, key)
    key = "counters/$key"
    if haskey(file, key)
        tmp = read(file, key) + 1
        write(file[key], tmp)
        tmp
    else
        file[key] = 0
    end
end


text:
__precompile__()
module LibIAEAPhsp

include("constants.jl")
include("raw.jl")

end # module


text:
module Test_efdsyn

using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Random
using Test

# Testing examples for EFDSYN
println("Test_efdsyn")
@testset "efdsyn" begin
rand()

## Example without control and disturbance inputs
p = 3; mf = 2;
sysf = fdimodset(rss(1,p,mf),faults = 1:mf);

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = false, rdim = 2); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys[:,Rf.faults]-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[1 1; 1 1; 1 1] && info.HDesign == [0.0 1.0 0.0; 0.0 0.0 1.0]

## Example without disturbance inputs
p = 3; mf = 2; mu = 2
sysf = fdimodset(rss(1,p,mf+mu),controls = 1:mu, faults = mu+1:mu+mf);

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; FDtol = 0.00001, atol = 1.e-7, minimal = false, rdim = 2); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && iszero(R.sys[:,R.controls],atol=1.e-7) &&
      info.S == Bool[1 1; 1 1; 1 1] && info.HDesign == [0.0 1.0 0.0; 0.0 0.0 1.0]

# solve using observer based nullspace
@time Q, Rf, info = efdsyn(sysf; nullspace = false, atol3 = 1.e-7, minimal = false, rdim = 2); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && iszero(R.sys[:,R.controls],atol=1.e-7) &&
      info.S == Bool[1 1; 1 1; 1 1] && info.HDesign == [1.0 0.0 0.0; 0.0 1.0 0.0]

@time Q, Rf, info = efdsyn(sysf; nullspace = false, atol3 = 1.e-7, minimal = true, rdim = 1); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && iszero(R.sys[:,R.controls],atol=1.e-7) &&
      info.S == Bool[1 1; 1 1; 1 1] && info.HDesign == [1.0 0.0 0.0]

## Example with solvable strong synthesis (frequencies are also zeros)
#
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s/(s+1) 0; (s^2+4)/(s+3)^2 1]; # enter Gf(s)
# new information in smat and zero local H-minus index values
smat, gains = fdisspec_(dss(Gf), [0, 2]; FDGainTol = 1.e-3)
@test smat[:,:,1] == Bool[0 0; 1 1] && smat[:,:,2] == Bool[1 0; 0 1] &&
      gains[:,:,1] ≈ [0.0 0.0; 0.4444444444444444 1.0] && gains[:,:,2] ≈ [0.894427190999916 0.0; 0.0 1.0]
ssmat = fditspec_(dss(Gf), atol = 1.e-7,FDtol = 1.e-3, FDfreq = [0, 2])
@test smat == ssmat
swmat = fditspec_(dss(Gf), atol=1.e-7,FDtol = 1.e-3)
@test swmat == Bool[1 0; 1 1]

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1);

@time Q, Rf, info = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2]); info
# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2])[:] == Bool[1]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2], FDGainTol = 0.001, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1:2);
@time Q, Rf, info = efdsyn(sysf; sdeg = -3, rdim = 2, FDfreq = [0, 2]); info
# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2]; block = false) == Bool[0 1; 0 0] &&
      fdisspec(Rf, [0]; block = false) == Bool[1 1; 0 0] &&
      fdisspec(Rf, [2]; block = false) == Bool[0 1; 1 0] &&
      fdisspec(Rf, [0, 2]; block = true) == Bool[1 1]

## Example with solvable strong synthesis (frequencies are also zeros)
#  Using simple basis based synthesis
#
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s/(s+1) 0; (s^2+4)/(s+3)^2 1]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1);

@time Q, Rf, info = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2], simple = true); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2])[:] == Bool[1]

Q1, Rf1, info1 = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2], FDGainTol = 0.001, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

## Example with solvable strong synthesis (frequencies are also poles)
#
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [1/s; 1/(s^2+4)]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1);

@time Q, Rf, info = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2]); info
# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2])[:] == Bool[1]

Q1, Rf1, info1 = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2], FDGainTol = 0.001, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

## Example with solvable strong synthesis (frequencies are also poles)
#  Using simple basis based synthesis
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [1/s; 1/(s^2+4)]; # enter Gf(s)

@time Q, Rf, info = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2], simple = true); info
# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2])[:] == Bool[1]

Q1, Rf1, info1 = efdsyn(sysf; sdeg = -3, rdim = 1, FDfreq = [0, 2], FDGainTol = 0.001, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

## Example with solvable strong synthesis (frequencies are also zeros)
#
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s; s^2+4]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1);

@time Q, Rf, info = efdsyn(sysf; sdeg = -3, smarg = -3, rdim = 1, FDfreq = [0, 2]); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2])[:] == Bool[1]

Q1, Rf1, info1 = efdsyn(sysf; sdeg = -3, smarg = -3, rdim = 1, FDfreq = [0, 2], HDesign = info.HDesign); info1
@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

@time Q, Rf, info = efdsyn(sysf; sdeg = -3, smarg = -3, rdim = 1, FDfreq = [0, 2], simple = true); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2])[:] == Bool[1]

Q1, Rf1, info1 = efdsyn(sysf; sdeg = -3, smarg = -3, rdim = 1, FDfreq = [0, 2], HDesign = info.HDesign); info1
@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

## Example with solvable strong synthesis (frequencies are also zeros)
#  Using simple basis based synthesis
#
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s; s^2+4]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1);

@time Q, Rf, info = efdsyn(sysf; sdeg = -3, smarg = -3, rdim = 1, FDfreq = [0, 2], simple = true); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0, 2])[:] == Bool[1]

Q1, Rf1, info1 = efdsyn(sysf; sdeg = -3, smarg = -3, rdim = 1, FDfreq = [0, 2], simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

# Example 5.2 - EFDP has no solution

# define s as an improper transfer function
s = rtf('s');
# define Gu(s), Gd(s), Gf(s)
Gu = [(s+1)/(s+2); (s+2)/(s+3)];     # enter Gu(s)
Gd = [1/(s+2); 0];                   # enter Gd(s)
Gf = [(s+1)/(s+2) 0; 0 1];           # enter Gf(s)
p = 2; mu = 1; md = 1; mf = 2;       # set dimensions

# build model with additive faults
#sysf = fdimodset(ss([Gu Gd Gf]),struct('c',1:mu,'d',mu+(1:md),'f',mu+md+(1:mf)));
sysf = fdimodset(dss([Gu Gd Gf]),c = 1:mu, d = mu.+(1:md), f = (mu+md).+(1:mf))

# call of EFDSYN with default options
try
  Q, Rf, info = efdsyn(sysf)
  @test false
catch err
  display(err)
  # compute achievable structure matrix
  S_achievable = fdigenspec(sysf)
  @test true
end

# Example 5.3 - Solution of an EFDP using EFDSYN

# define s as an improper transfer function
s = rtf('s');
# define Gu(s) and Gd(s)
Gu = [(s+1)/(s+2); (s+2)/(s+3)];     # enter Gu(s)
Gd = [(s-1)/(s+2); 0];               # enter Gd(s)
p = 2; mu = 1; md = 1; mf = 2;       # set dimensions

# setup the synthesis model with additive faults with Gf(s) = [ Gu(s) [0;1]]
sysf = fdimodset(dss([Gu Gd]),c =1,d = 2,f = 1,fs = 2);

# call of EFDSYN with the options for stability degree -3 and the synthesis
# of a scalar output filter
@time Q, Rf, info = efdsyn(sysf, smarg =-3, sdeg = -3, rdim = 1);

# normalize Q and Rf to match example
scale = evalfr(Rf.sys[1,1],Inf)[1,1];
dss2rm(Q.sys/scale, atol = 1.e-7)
dss2rm(Rf.sys/scale, atol = 1.e-7)

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0]) == Bool[1 1]

# Example 5.4c - Solution of an EFDP using EFDSYN

# define s as an improper transfer function
s = rtf('s');
# define Gu(s) and Gd(s)
Gu = [(s+1)/(s-2); (s+2)/(s-3)];     # enter Gu(s)
Gd = [(s-1)/(s+2); 0];               # enter Gd(s)
Gf = [(s+1)/(s-2) 0; (s+2)/(s-3) 1]; # enter Gf(s)
p = 2; mu = 1; md = 1; mf = 2;       # set dimensions

# setup the synthesis model with faults with Gf(s) = [ Gu(s) [0;1]]
sysf = fdimodset(dss([Gu Gd]),c =1,d = 2,f = 1,fs = 2);

# call of EFDSYN with the options for stability degree -3 and the synthesis
# of a scalar output filter
@time Q, Rf, info = efdsyn(sysf, sdeg = -3, rdim = 1);

# normalize Q and Rf to match example
scale = evalfr(Rf.sys[1,1],Inf)[1,1];
dss2rm(Q.sys/scale, atol = 1.e-7)
dss2rm(Rf.sys/scale, atol = 1.e-7)

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(R.sys[:,[R.controls;R.disturbances]], atol = 1.e-7) &&
      iszero(R.sys[:,R.faults]-Rf.sys) &&
      fdisspec(Rf, [0]) == Bool[1 1]

@test iszero(gbalmr(Q; balance = true).sys-Q.sys, atol=1.e-7)

## Test covering design matrices
##
p = 3; mf = 2;
sysf = fdimodset(rss(1,p,mf),faults = 1:mf);

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = false, rdim = 2); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[1 1; 1 1; 1 1] && info.HDesign == [0.0 1.0 0.0; 0.0 0.0 1.0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

## solve an EFP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = false, rdim = 2, simple = true); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[1 1; 1 1; 1 1] && info.HDesign == [0.0 1.0 0.0; 0.0 0.0 1.0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

## solve an EFP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = false, rdim = 1, simple = true); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[1 1; 1 1; 1 1] && info.HDesign == [0.0  1.0  0.0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

## solve an EFD using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = true, simple = true); info
R = fdIFeval(Q, sysf; atol = 1.e-7);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && iszero(R.sys[:,R.controls],atol=1.e-7) &&
      info.S == Bool[1 1; 1 1; 1 1] && info.HDesign ==[0.0 1.0 0.0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

## solve an EFD using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = true, simple = false); info
R = fdIFeval(Q, sysf; atol = 1.e-7);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && iszero(R.sys[:,R.controls],atol=1.e-7) &&
      info.S == Bool[1 1; 1 1; 1 1] && info.HDesign ==[0.0 1.0 0.0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, simple = false, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

##
## Example with fault and noise inputs
p = 3; mf = 2; mw = 2; mu = 0
sysf = fdimodset(rss(1,p,mf+mw), noise = (mu+mf).+ (1:mw), faults = mu+1:mu+mf);

## solve an AMMP using the nullspace based approach
Q, Rfw, info = efdsyn(sysf; atol = 1.e-7, atol3 = 1.e-7); info
R = fdIFeval(Q, sysf; atol = 1.e-7);
@test iszero(Rfw.sys[:,Rfw.faults]-R.sys[:,R.faults],atol=1.e-7) && iszero(R.sys[:,R.controls],atol=1.e-7) &&
      iszero(Rfw.sys[:,Rfw.noise]-R.sys[:,R.noise],atol=1.e-7) &&
      info.S == Bool[1 1; 1 1; 1 1] && info.HDesign ==[0.0 1.0 0.0]

# the same results  must result using the resulting design matrix H
Q1, Rfw1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, simple = false, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rfw.sys-Rfw1.sys, atol = 1.e-7)

## special cases which require forming rdim combinations of nout vecotors, with nout > rdim

# Gf = [ I;0]
p = 4; mf = 2;
sys = dss([eye(2);zeros(2,2)]);
sysf = fdimodset(sys,faults = 1:mf);

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = true, rdim = 1); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[0 0; 0 0; 0 1; 1 0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = false, rdim = 1); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[0 0; 0 0; 0 1; 1 0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = true, rdim = 1, simple = true); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[0 0; 0 0; 0 1; 1 0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = true, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = efdsyn(sysf; atol3 = 1.e-7, minimal = false, rdim = 1, simple = true); info
R = fdIFeval(Q, sysf);
@test iszero(Rf.sys-R.sys[:,R.faults],atol=1.e-7) && info.S == Bool[0 0; 0 0; 0 1; 1 0]

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = efdsyn(sysf; atol3 = 1.e-7, minimal = false, simple = true, HDesign = info.HDesign); info1

@test iszero(Q.sys-Q1.sys, atol = 1.e-7) && iszero(Rf.sys-Rf1.sys, atol = 1.e-7)

## Model Niemann 1998, Optim. Appl. Meth.
n = 5; mu = 2; md = 1; mf = 2; m = mu+md+mf; p = 5;
a = [
-0.0782 0.2939 0.0220 0.0208 -0.0291
0.0077 -0.0278 0.0015 -0.0015 0.0026
-1.3251 5.32608 -0.5263 0.2214 -0.4777
1.0809 -4.4452 0.3770 -0.4631 0.4032
2.1532 -8.6386 0.7811 -0.5745 0.7816]*1000;
bd = [
    0.0443
   -0.0042
    0.7910
   -0.6598
   -1.2881]*1000;
bu = [
    -0.007197 0.003005
0.003473 0.000284
1.218502 -0.032866
1.322502 0.020147
-0.082265 0.024388 ];
bf = bu;
c = eye(n); dd = zeros(n,md); du = zeros(n,mu); df = zeros(n,mf);

sys = dss(a,[bu bd bf],c,[du dd df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

# opt = struct('tol',1.e-7,'FDGainTol',.01,...
#     'FDFreq',0,'sdeg',-0.05,'m1',mu+md);
# apply fdigenspec to [Gu Gd Gf; I 0 0]
#S_strong = fdigenspec([sys; eye(mu,mu+md+mf)],opt)
S_strong = Bool[1 0;0 1;1 1]
SFDI = S_strong

nb = size(SFDI,1)
orders = zeros(Int,nb)

Qtot = dss(zeros(0,p+mu)); Rftilde = dss(zeros(0,mf));
info = Any[]; HDesign = Any[]
for i = 1:nb
   # set input groups for i-th design
   indd = Vector(1:mf)[SFDI[i,:] .== false]
   indf = Vector(1:mf)[SFDI[i,:] .== true]
   sysc = fdimodset(sysf.sys, c = 1:mu, d = [Vector(mu .+ (1:md)); (mu+md) .+ indd],
                    f = (mu+md) .+ indf, aux = (mu+md) .+ (1:mf));
   #
   Qi, Rfi, infoi = efdsyn(sysc, atol = 1.e-7, sdeg =-5, smarg = -0.05, minimal = true,
                           FDfreq = 0, FDGainTol = 0.0001, rdim = 1, simple = false);  infoi
   info = [info; infoi]
   HDesign = [HDesign; [infoi.HDesign]]
   Qtot = [Qtot; Qi.sys];
   Rftilde = [Rftilde; Rfi.sys[:,Rfi.aux]];
   orders[i] = order(Qi.sys)
end

# checks: Q[Gu Gd;I 0] = 0 and Q[Gf; 0] = Rftilde
syse = [sysf.sys; eye(mu,mu+md+mf)];
@test iszero(Qtot*syse[:,[sysf.controls; sysf.disturbances]], atol = 1.e-7)
@test iszero(Qtot*syse[:,sysf.faults]-Rftilde,atol=1.e-7)

@test all(orders .== 1)
#step(Rftilde,10)

S, gains = fdisspec_(Rftilde, 0;  FDGainTol = 0.0001)
@test isequal(SFDI,S[:,:,1])
S1 = fditspec_(Rftilde; atol = 1.e-7,FDStol = 0.0001, FDfreq = 0)
@test isequal(SFDI,S1[:,:,1])

## Example 1, Xiong & Saif 2000, Int. J. Robust Nonlinear Control
n = 4; mu = 2; md = 2; mf = 2; m = mu+md+mf; p = 4;
a = [ -9.9477 -0.7476 0.2632 -5.0337
      52.1659 2.7452 5.5532 -24.4221
      26.0922 2.6361 -4.1975 -19.2774
      0.0 0.0 1.0 0.0];
bu = [ 0.4422 0.1761
       3.5446 -7.5922
        -5.52 4.49
        0.0 0.0];
bd = [ 0 0; 0 1 ; 1 0; 0 0 ];
bf = [bu[:,1] zeros(n,1)];
c = [ 1 0 0 0; 0 1 0 0; 0 0 1 0; 0 1 1 1];
dd = zeros(n,md); du = zeros(n,mu); df = [zeros(n,mf-1) [0 0 1 0]'];

sys = dss(a,[bu bd bf],c,[du dd df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

# opt = struct('tol',1.e-7,'FDGainTol',.01,...
#     'FDFreq',0,'sdeg',-0.05,'m1',mu+md);
# apply fdigenspec to [Gu Gd Gf; I 0 0]
#S_strong = fdigenspec([sys; eye(mu,mu+md+mf)],opt)
S_strong = Bool[1 0;0 1;1 1]
SFDI = S_strong

nb = size(SFDI,1)
orders = zeros(Int,nb)

Qtot = dss(zeros(0,p+mu)); Rftilde = dss(zeros(0,mf));
info = Any[]; HDesign = Any[]
for i = 1:nb
   # set input groups for i-th design
   indd = Vector(1:mf)[SFDI[i,:] .== false]
   indf = Vector(1:mf)[SFDI[i,:] .== true]
   sysc = fdimodset(sysf.sys, c = 1:mu, d = [Vector(mu .+ (1:md)); (mu+md) .+ indd],
                    f = (mu+md) .+ indf, aux = (mu+md) .+ (1:mf));
   #
   Qi, Rfi, infoi = efdsyn(sysc, atol = 1.e-7, sdeg =-5, smarg = -0.05, minimal = true,
                           FDfreq = 0, FDGainTol = 0.0001, rdim = 1, simple = false);  infoi
   info = [info; infoi]
   HDesign = [HDesign; [infoi.HDesign]]
   Qtot = [Qtot; Qi.sys];
   Rftilde = [Rftilde; Rfi.sys[:,Rfi.aux]];
   orders[i] = order(Qi.sys)
end

# checks: Q[Gu Gd;I 0] = 0 and Q[Gf; 0] = Rftilde
syse = [sysf.sys; eye(mu,mu+md+mf)];
@test iszero(Qtot*syse[:,[sysf.controls; sysf.disturbances]], atol = 1.e-7)
@test iszero(Qtot*syse[:,sysf.faults]-Rftilde,atol=1.e-7)

@test orders == [2, 1, 1]
#step(Rftilde,10)

S, gains = fdisspec_(Rftilde, 0;  FDGainTol = 0.0001)
@test isequal(SFDI,S[:,:,1])
S1 = fditspec_(Rftilde; atol = 1.e-7,FDStol = 0.0001, FDfreq = 0)
@test isequal(SFDI,S1[:,:,1])

## Example 2, Xiong & Saif 2000, Int. J. Robust Nonlinear Control
n = 5; mu = 3; md = 2; mf = 2; m = mu+md+mf; p = 3;
a = [
0 0 -0.0034 0 0
0 -0.041 0.0013 0 0
0 0 -1.1471 0 0
0 0 -0.0036 0 0
0 0.094 0.0057 0 -0.051];
bu = [
    -1 0 0
     0 0 0
     0 0 0.948
     0.916 -1 0
    -0.598 0 0];
bd = [
     0 1
     0.062 -0.132
     0 -7.189
     0 0
     0 0 ];
bf = bu[:,1:2];
c = [
1 0 0 0 0
0 0 1 0 0
0 0 0 1 0];
dd = zeros(p,md); du = zeros(p,mu); df = zeros(p,mu-1);

sys = dss(a,[bu bd bf],c,[du dd df]);

sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

# opt = struct('tol',1.e-7,'FDGainTol',.01,...
#     'FDFreq',0,'sdeg',-0.05,'m1',mu+md);
# apply fdigenspec to [Gu Gd Gf; I 0 0]
#S_strong = fdigenspec([sys; eye(mu,mu+md+mf)],opt)
S_strong = Bool[1 0;0 1;1 1]
SFDI = S_strong

nb = size(SFDI,1)
orders = zeros(Int,nb)

Qtot = dss(zeros(0,p+mu)); Rftilde = dss(zeros(0,mf));
info = Any[]; HDesign = Any[]
for i = 1:nb
   # set input groups for i-th design
   indd = Vector(1:mf)[SFDI[i,:] .== false]
   indf = Vector(1:mf)[SFDI[i,:] .== true]
   sysc = fdimodset(sysf.sys, c = 1:mu, d = [Vector(mu .+ (1:md)); (mu+md) .+ indd],
                    f = (mu+md) .+ indf, aux = (mu+md) .+ (1:mf));
   #
   Qi, Rfi, infoi = efdsyn(sysc, atol = 1.e-7, sdeg =-5, smarg = -0.05, minimal = true,
                           FDfreq = 0, FDGainTol = 0.0001, rdim = 1, simple = false);  infoi
   info = [info; infoi]
   HDesign = [HDesign; [infoi.HDesign]]
   Qtot = [Qtot; Qi.sys];
   Rftilde = [Rftilde; Rfi.sys[:,Rfi.aux]];
   orders[i] = order(Qi.sys)
end

# checks: Q[Gu Gd;I 0] = 0 and Q[Gf; 0] = Rftilde
syse = [sysf.sys; eye(mu,mu+md+mf)];
@test iszero(Qtot*syse[:,[sysf.controls; sysf.disturbances]], atol = 1.e-7)
@test iszero(Qtot*syse[:,sysf.faults]-Rftilde,atol=1.e-7)

@test orders == [1, 1, 1]
#step(Rftilde,10)

S, gains = fdisspec_(Rftilde, 0;  FDGainTol = 0.0001)
@test isequal(SFDI,S[:,:,1])
S1 = fditspec_(Rftilde; atol = 1.e-7,FDStol = 0.0001, FDfreq = 0)
@test isequal(SFDI,S1[:,:,1])

# use individual synthesis with specified structure vectors
nb = size(SFDI,1)
orders = zeros(Int,nb)

Qtot = dss(zeros(0,p+mu)); Rftilde = dss(zeros(0,mf));
info = Any[]; HDesign = Any[]
for i = 1:nb
   Qi, Rfi, infoi = efdsyn(sysf, SFDI[i,:]; atol = 1.e-7, sdeg =-5, smarg = -0.05, minimal = true,
                           FDfreq = 0, FDGainTol = 0.0001, rdim = 1, simple = false);  infoi
   info = [info; infoi]
   HDesign = [HDesign; [infoi.HDesign]]
   Qtot = [Qtot; Qi.sys];
   Rftilde = [Rftilde; Rfi.sys];
   orders[i] = order(Qi.sys)
end

# checks: Q[Gu Gd;I 0] = 0 and Q[Gf; 0] = Rftilde
syse = [sysf.sys; eye(mu,mu+md+mf)];
@test iszero(Qtot*syse[:,[sysf.controls; sysf.disturbances]], atol = 1.e-7)
@test iszero(Qtot*syse[:,sysf.faults]-Rftilde,atol=1.e-7)

@test orders == [1, 1, 1]
#step(Rftilde,10)

S, gains = fdisspec_(Rftilde, 0;  FDGainTol = 0.0001)
@test isequal(SFDI,S[:,:,1])
S1 = fditspec_(Rftilde; atol = 1.e-7,FDStol = 0.0001, FDfreq = 0)
@test isequal(SFDI,S1[:,:,1])

## Yuan et al. IJC (1997)  strong synthesis
p = 3; mu = 1; md = 0; mf = 8; n = 4;
A = [
-1 1 0 0
1 -2 1 0
0 1 -2 1
0 0 1 -2  ];
Bu = [1 0 0 0]';
Bf = [
1 0 0 0 1 0 0 0
0 1 0 0 -1 1 0 0
0 0 1 0 0 -1 1 0
0 0 0 1 0 0 -1 1];
C = [
 1 0 0 0
0 0 1 0
0 0 0 1];
Du = zeros(p,mu); Df = zeros(p,mf);

sys = dss(A,[Bu Bf],C,[Du Df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

# opt = struct('tol',1.e-7,'FDGainTol',.01,...
#     'FDFreq',0,'sdeg',-0.05,'m1',mu+md);
# # apply fdigenspec to [Gu Gd Gf; I 0 0]
# S_strong = fdigenspec([sys; eye(mu,mu+md+mf)],opt)
#

S_strong = Bool[   0   0   0   1   0   0   1   1
0   1   1   0   1   1   1   0
0   1   1   1   1   1   0   1
0   1   1   1   1   1   1   1
1   0   1   0   1   1   1   0
1   0   1   1   1   1   0   1
1   0   1   1   1   1   1   1
1   1   0   0   1   1   0   0
1   1   0   1   1   1   1   1
1   1   1   0   1   1   1   0
1   1   1   1   1   1   0   1
1   1   1   1   1   1   1   1]

SFDI = S_strong

nb = size(SFDI,1)
orders = zeros(Int,nb)

Qtot = dss(zeros(0,p+mu)); Rftilde = dss(zeros(0,mf));
info = Any[]; HDesign = Any[]
for i = 1:nb
   # set input groups for i-th design
   indd = Vector(1:mf)[SFDI[i,:] .== false]
   indf = Vector(1:mf)[SFDI[i,:] .== true]
   sysc = fdimodset(sysf.sys, c = 1:mu, d = [Vector(mu .+ (1:md)); (mu+md) .+ indd],
                    f = (mu+md) .+ indf, aux = (mu+md) .+ (1:mf));
   #
   Qi, Rfi, infoi = efdsyn(sysc, atol = 1.e-7, sdeg =-5, smarg = -0.05, minimal = true,
                           FDfreq = 0, FDGainTol = 0.0001, rdim = 1, simple = false);  infoi
   info = [info; infoi]
   HDesign = [HDesign; [infoi.HDesign]]
   Qtot = [Qtot; Qi.sys];
     Rftilde = [Rftilde; Rfi.sys[:,Rfi.aux]];
   orders[i] = order(Qi.sys)
end

# checks: Q[Gu Gd;I 0] = 0 and Q[Gf; 0] = Rftilde
syse = [sysf.sys; eye(mu,mu+md+mf)];
@test iszero(Qtot*syse[:,[sysf.controls; sysf.disturbances]], atol = 1.e-7)
@test iszero(Qtot*syse[:,sysf.faults]-Rftilde,atol=1.e-7)

@test orders == [1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2]
#step(Rftilde,10)

S, gains = fdisspec_(Rftilde, 0;  FDGainTol = 0.0001)
@test isequal(SFDI,S[:,:,1])
S1 = fditspec_(Rftilde; atol = 1.e-7,FDStol = 0.0001, FDfreq = 0)
@test isequal(SFDI,S1[:,:,1])

## Yuan et al. IJC (1997)  weak synthesis with stabillization
p = 3; mu = 1; md = 0; mf = 8; n = 4;
A = [
-1 1 0 0
1 -2 1 0
0 1 -2 1
0 0 1 -2  ];
Bu = [1 0 0 0]';
Bf = [
1 0 0 0 1 0 0 0
0 1 0 0 -1 1 0 0
0 0 1 0 0 -1 1 0
0 0 0 1 0 0 -1 1];
C = [
 1 0 0 0
0 0 1 0
0 0 0 1];
Du = zeros(p,mu); Df = zeros(p,mf);

sys = dss(A,[Bu Bf],C,[Du Df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

# opt = struct('tol',1.e-7,'FDGainTol',.01,...
#     'FDFreq',0,'sdeg',-0.05,'m1',mu+md);
# # apply fdigenspec to [Gu Gd Gf; I 0 0]
# S_strong = fdigenspec([sys; eye(mu,mu+md+mf)],opt)
#

S_weak = Bool[   0   0   0   1   0   0   1   1
0   1   1   0   1   1   1   0
0   1   1   1   1   1   0   1
0   1   1   1   1   1   1   1
1   0   1   0   1   1   1   0
1   0   1   1   1   1   0   1
1   0   1   1   1   1   1   1
1   1   0   0   1   1   0   0
1   1   0   1   1   1   1   1
1   1   1   0   0   1   1   0
1   1   1   0   1   0   1   0
1   1   1   0   1   1   1   0
1   1   1   1   0   1   0   1
1   1   1   1   0   1   1   1
1   1   1   1   1   0   0   1
1   1   1   1   1   0   1   1
1   1   1   1   1   1   0   1
1   1   1   1   1   1   1   1]

SFDI = S_weak;
nb = size(SFDI,1)
orders = zeros(Int,nb)

Qtot = dss(zeros(0,p+mu)); Rftilde = dss(zeros(0,mf));
info = Any[]; HDesign = Any[]
for i = 1:nb
   # set input groups for i-th design
   indd = Vector(1:mf)[SFDI[i,:] .== false]
   indf = Vector(1:mf)[SFDI[i,:] .== true]
   sysc = fdimodset(sysf.sys, c = 1:mu, d = [Vector(mu .+ (1:md)); (mu+md) .+ indd],
                    f = (mu+md) .+ indf, aux = (mu+md) .+ (1:mf));
   #
   Qi, Rfi, infoi = efdsyn(sysc, atol = 1.e-7, sdeg =-5, smarg = -0.05, minimal = true,
                           FDtol = 0.0001, rdim = 1, simple = false);  infoi
   info = [info; infoi]
   HDesign = [HDesign; [infoi.HDesign]]
   Qtot = [Qtot; Qi.sys];
   Rftilde = [Rftilde; Rfi.sys[:,Rfi.aux]];
   orders[i] = order(Qi.sys)
end

# checks: Q[Gu Gd;I 0] = 0 and Q[Gf; 0] = Rftilde
syse = [sysf.sys; eye(mu,mu+md+mf)];
@test iszero(Qtot*syse[:,[sysf.controls; sysf.disturbances]], atol = 1.e-7)
@test iszero(Qtot*syse[:,sysf.faults]-Rftilde,atol=1.e-7)

@test orders == [1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
#step(Rftilde,10)

S = fditspec_(Rftilde;  atol = 1.e-6, FDtol = 0.00001)
@test isequal(SFDI,S)

## Yuan et al. IJC (1997)  weak synthesis with pole assignment and
#  enforcing least global order of 6!
p = 3; mu = 1; md = 0; mf = 8; n = 4;
A = [
-1 1 0 0
1 -2 1 0
0 1 -2 1
0 0 1 -2  ];
Bu = [1 0 0 0]';
Bf = [
1 0 0 0 1 0 0 0
0 1 0 0 -1 1 0 0
0 0 1 0 0 -1 1 0
0 0 0 1 0 0 -1 1];
C = [
 1 0 0 0
0 0 1 0
0 0 0 1];
Du = zeros(p,mu); Df = zeros(p,mf);

sys = dss(A,[Bu Bf],C,[Du Df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

# opt = struct('tol',1.e-7,'FDGainTol',.01,...
#     'FDFreq',0,'sdeg',-0.05,'m1',mu+md);
# # apply fdigenspec to [Gu Gd Gf; I 0 0]
# S_strong = fdigenspec([sys; eye(mu,mu+md+mf)],opt)
#

S_weak = Bool[   0   0   0   1   0   0   1   1
0   1   1   0   1   1   1   0
0   1   1   1   1   1   0   1
0   1   1   1   1   1   1   1
1   0   1   0   1   1   1   0
1   0   1   1   1   1   0   1
1   0   1   1   1   1   1   1
1   1   0   0   1   1   0   0
1   1   0   1   1   1   1   1
1   1   1   0   0   1   1   0
1   1   1   0   1   0   1   0
1   1   1   0   1   1   1   0
1   1   1   1   0   1   0   1
1   1   1   1   0   1   1   1
1   1   1   1   1   0   0   1
1   1   1   1   1   0   1   1
1   1   1   1   1   1   0   1
1   1   1   1   1   1   1   1]

SFDI = S_weak;
nb = size(SFDI,1)
orders = zeros(Int,nb)

Qtot = dss(zeros(0,p+mu)); Rftilde = dss(zeros(0,mf));
info = Any[]; HDesign = Any[]
for i = 1:nb
   # set input groups for i-th design
   indd = Vector(1:mf)[SFDI[i,:] .== false]
   indf = Vector(1:mf)[SFDI[i,:] .== true]
   sysc = fdimodset(sysf.sys, c = 1:mu, d = [Vector(mu .+ (1:md)); (mu+md) .+ indd],
                    f = (mu+md) .+ indf, aux = (mu+md) .+ (1:mf));
   #
   Qi, Rfi, infoi = efdsyn(sysc, atol = 1.e-7, sdeg =-5, smarg = -5, poles = [-5, -6], minimal = true,
                            FDtol = 0.0001, rdim = 1, simple = false);  infoi
   info = [info; infoi]
   HDesign = [HDesign; [infoi.HDesign]]
   Qtot = [Qtot; Qi.sys];
   Rftilde = [Rftilde; Rfi.sys[:,Rfi.aux]];
   orders[i] = order(Qi.sys)
end

# checks: Q[Gu Gd;I 0] = 0 and Q[Gf; 0] = Rftilde
syse = [sysf.sys; eye(mu,mu+md+mf)];
@test iszero(Qtot*syse[:,[sysf.controls; sysf.disturbances]], atol = 1.e-7)
@test iszero(Qtot*syse[:,sysf.faults]-Rftilde,atol=1.e-7)

@test orders == [1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
#step(Rftilde,10)

S = fditspec_(Rftilde;  atol = 1.e-7, FDtol = 0.0001);
@test isequal(SFDI,S)

@test  count(ghanorm(Qtot)[2] .> 1.e-7) == 6

# ## DC Motor (Ding, 2013)
# s = rtf('s');
# Gu = 8.75/(1+1.225*s)/(1+0.03*s)/(1+0.005*s);
# Gd = -31.07/s/(1+0.005*s);
# Gc = 1.6*(1 + 1.225*s)/s;

# G = [Gu*Gc Gu 1 Gd ]./(1+Gu*Gc);
# # set input groups
# sysf = fdimodset(dss(G,minimal=true,atol=1.e-7),c = 1, f = 2:4);
# S_weak = fdigenspec(sysf,struct('tol',1.e-7))
# S_strong = fdigenspec(sysf,struct('tol',1.e-7,'FDFreq',0))

# # # check with expressions from Ding (2013)
# # gu = 14/(s*(1+0.03*s)*(1+0.005*s)+14);
# # gd = -31.07*(1+0.03*s)/(s*(1+0.03*s)*(1+0.005*s)+14);
# # gfa = 8.75*s/(1+1.225*s)/(s*(1+0.03*s)*(1+0.005*s)+14);
# # gfs = s*(1+0.03*s)*(1+0.005*s)/(s*(1+0.03*s)*(1+0.005*s)+14);
# # norm(gu-G(:,1),inf), norm(gd-G(:,4),inf), norm(gfa-G(:,2),inf), norm(gfs-G(:,3),inf),
# #

# ## case 0:
# f = [0 0 1]; T = [0:0.1:5]'; N = length(T); ucontrol = ones(N,1);
# Uf = [ucontrol, f(1)*ones(N,1) f(2)*ones(N,1) f(3)*ones(N,1)];
# lsim(G,Uf,T);

# ## case 1 (Echeveria, DC Motor Benchmark):
# f = [0.87 -0.12 0.53]; T = [0:0.1:5]'; N = length(T); ucontrol = rand(N,1);
# Uf = [ucontrol, f(1)*ones(N,1) f(2)*ones(N,1) f(3)*ones(N,1)];
# yf = lsim(G,Uf,T);
# yfp = yf.*(ones(N,1)+0.02*rand(N,1));
# yfp2 = yf.*(ones(N,1)+0.08*rand(N,1));

# x = rand(3,1);
# err = dcmotor( x, G, yf, T, ucontrol);
# lb = -ones(3,1); ub = ones(3,1);

# [x,err] = fmincon(@(x) dcmotor(x,G,yf,T,ucontrol),zeros(3,1),[],[],[],[],lb,ub)

# [xp,err1] = fmincon(@(x) dcmotor(x,G,yfp,T,ucontrol),zeros(3,1),[],[],[],[],lb,ub)
# [xp2,err2] = fmincon(@(x) dcmotor(x,G,yfp2,T,ucontrol),zeros(3,1),[],[],[],[],lb,ub)

end # test efdsyn

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
function instantaneous_population(
    sol::ODESolution,
    H::OpenQuantumBase.AbstractHamiltonian;
    lvl = 1,
    t_axis = sol.t,
)
    y = []
    if ndims(sol.u[1]) == 1
        for x in t_axis
            s = sol.prob.p(x)
            _, v = eigen_decomp(H, s, lvl = lvl)
            push!(y, abs2.(v' * sol(x)))
        end
    elseif ndims(sol.u[1]) == 2
        for x in t_axis
            s = sol.prob.p(x)
            _, v = eigen_decomp(H, s, lvl = lvl)
            push!(y, real.(diag(v' * sol(x) * v)))
        end
    else
        throw(ArgumentError("Solution needs to be state vector or density matrix."))
    end
    y = hcat(y...)'
end

pulse_on_density!(c::Vector, p) = c .= conj(p) ⊗ p * c
pulse_on_density!(c::Matrix, p) = c .= p * c * p'


text:
using Maker
using Base.Test

COUNT = 0

path = "tmpdir/tmp"

Maker.directory(path)

Maker.file("$path/x.csv", path) do t
    global COUNT += 1
    writecsv(t.name, rand(2,2))
end

@desc "Update y.csv"
function update_y_csv(t)
    global COUNT += 1
    y = readcsv(t.dependencies[1])
    writecsv(t.name, 2 * y)
end
Maker.file(update_y_csv, "$path/y.csv", "$path/x.csv")
show(tasks("$path/y.csv"))

Maker.task("default", "$path/y.csv")

Maker.task("clean") do
    rm("tmpdir", recursive = true)
end

make()

@test COUNT == 2
@test isdir("tmpdir")
@test isdir("tmpdir/tmp")
@test isfile("tmpdir/tmp/x.csv")
@test isfile("tmpdir/tmp/y.csv")

make()
make()
@test COUNT == 2

rm("tmpdir/tmp", recursive = true)
make()
@test COUNT == 4

make("clean")

@test !isdir("tmpdir")
@test !isdir("tmpdir/tmp")
@test !isfile("tmpdir/tmp/x.csv")
@test !isfile("tmpdir/tmp/y.csv")

Maker.fixcache()


text:
# Invertible network utilities
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August 2020

export get_pars, get_grads, set_pars!, tmapz2x_den, initialize_id!

function get_pars(H::InvertibleNetwork)
    θ = get_params(H)
    θvec = Array{Array{Float32}, 1}(undef, length(θ))
    for i = 1:length(θ)
        θvec[i] = θ[i].data
    end
    return θvec
end

function set_pars!(H::InvertibleNetwork, θ::Array)
    θH = get_params(H)
    for i = 1:length(θ)
        isa(θ[i], Parameter) ? θH[i].data = copy(θ[i].data) : θH[i].data = copy(θ[i])
    end
end

function get_grads(H::InvertibleNetwork)
    θ = get_params(H)
    θvec = Array{Array{Float32}, 1}(undef, length(θ))
    for i = 1:length(θ)
        θvec[i] = θ[i].grad
    end
    return θvec
end

function tmapz2x_den(T::InvertibleNetwork, X::Array{Float32, 4})
    Z = T.inverse(X)
    logdet = Array{Float32, 4}(undef, 1, 1, 1, size(X, 4))
    for i = 1:size(X, 4)
        _, logdet[i] = T.forward(Z[:, :, :, i:i])
    end
    return neglogden = 0.5f0*sum(abs.(Z).^2f0, dims = 3)+logdet
end

function initialize_id!(H::NetworkHINT1D; αw::Float32 = 0f0, αb::Float32 = 0f0)
    n = length(H.CL)
    m = length(H.CL[1].CL)
    for i = 1:n
        for j = 1:m
            # H.CL[i].CL[j].RB.W1.data .*= αw
            # H.CL[i].CL[j].RB.W2.data .*= αw
            H.CL[i].CL[j].RB.W3.data .*= 0f0
            # H.CL[i].CL[j].RB.b1.data .*= αb
            # H.CL[i].CL[j].RB.b2.data .*= αb
       end
   end
   # H.CL[n].CL[1].RB.W1.data .= 0f0
   # H.CL[n].CL[1].RB.W2.data .= 0f0
   # H.CL[n].CL[1].RB.W3.data .= 0f0
end

function initialize_id!(H::NetworkGlow1D; αw::Float32 = 0f0, αb::Float32 = 0f0)
    n = length(H.CL)
    for i = 1:n
        # H.CL[i].RB.W1.data .*= αw
        # H.CL[i].RB.W2.data .*= αw
        H.CL[i].RB.W3.data .*= 0f0
        # H.CL[i].RB.b1.data .*= αb
        # H.CL[i].RB.b2.data .*= αb
   end
   # H.CL[n].CL[1].RB.W1.data .= 0f0
   # H.CL[n].CL[1].RB.W2.data .= 0f0
   # H.CL[n].CL[1].RB.W3.data .= 0f0
end


text:
using HighlyAdaptiveLasso
using Test
using Random
using MLJBase
using RCall

@testset "Test HALRegressor" begin
    # Those are sanity checks as the original package is already well tested
    n = 100
    p = 3
    X = randn(n, p)
    y = X[:, 1] .* 3sin.(X[:, 2]) + randn(n)
    Xt = MLJBase.table(X)
    # With cross validation config
    model = HALRegressor(cv_select=true, n_folds=3)
    mach = machine(model, Xt, y)
    fit!(mach)
    @test mach.state == 1
    @test predict(mach, Xt) isa Vector{Float64}

    @rget n_folds
    @test n_folds == 3

    @test mach.report ==
        "fit_hal(X, y, family = 'gaussian', fit_control = list(cv_select = cv_select, n_folds = n_folds))"

    # Test specifying a lambda
    model = HALRegressor(cv_select=false, lambda=30, max_degree=2, num_knots=[100, 50, 20], smoothness_orders=2)
    mach = machine(model, Xt, y)
    fit!(mach)
    @test mach.state == 1
    @test predict(mach, Xt) isa Vector{Float64}

    @rget max_degree
    @rget num_knots
    @rget smoothness_orders
    @rget cv_select

    @test cv_select == false
    @test max_degree == 2
    @test smoothness_orders == 2
    @test num_knots == [100, 50, 20]

    @test mach.report ==
        "fit_hal(X, y, family = 'gaussian', max_degree = max_degree, smoothness_orders = smoothness_orders, "*
        "num_knots = num_knots, lambda = lambda, fit_control = list(cv_select = cv_select))"

    # Test specifying a formula
    formula_ = "~h(x1) + h(x2, x1) + h(x3)"
    model = HALRegressor(formula=formula_, lambda=30, cv_select=false)
    mach = machine(model, Xt, y)
    fit!(mach)
    @test mach.state == 1
    @test predict(mach, Xt) isa Vector{Float64}

    @rget formula
    @test formula == formula_

    @test mach.report ==
        "fit_hal(X, y, family = 'gaussian', formula = formula, lambda = lambda, fit_control = list(cv_select = cv_select))"

end

@testset "Test HALClassifier" begin
    # Those are sanity checks as the original package is already well tested
    n = 100
    p = 3
    X = randn(n, p)
    y = categorical(rand(["toto", "tata"], n))
    Xt = MLJBase.table(X)
    # With cross validation config
    model = HALClassifier()
    mach = machine(model, Xt, y)
    fit!(mach)
    @test mach.state == 1
    @test predict(mach, Xt) isa UnivariateFiniteVector

    @test mach.report == "fit_hal(X, y, family = 'binomial')"

end


text:
module constrsparsereg_test

using Base.Test, ConstrainedLasso, ECOS

info("Test lsq_constrsparsereg: sum-to-zero constraint")

# set up
srand(123)
n, p = 100, 20
# truth with sum constraint sum(β) = 0
β = zeros(p)
β[1:round(Int, p / 4)] = 0
β[(round(Int, p / 4) + 1):round(Int, p / 2)] = 1
β[(round(Int, p / 2) + 1):round(Int, 3p / 4)] = 0
β[(round(Int, 3p / 4) + 1):p] = -1
# generate data
X = randn(n, p)
y = X * β + randn(n)
Aeq = ones(1, p)
beq = [0.0]
penwt  = ones(p)
solver = ECOSSolver(maxit=10e8, verbose=0)

info("Optimize at a single tuning parameter value")
ρ = 10.0
β̂, = lsq_constrsparsereg(X, y, ρ; Aeq = Aeq, beq = beq,
    penwt = penwt, solver = solver)
@test sum(β̂)≈0.0 atol=1e-6

info("Optimize at multiple tuning parameter values")
ρlist = [0.0:10.0; Inf]
β̂, = lsq_constrsparsereg(X, y, ρlist; Aeq = Aeq, beq = beq,
    penwt = penwt, solver = solver)
@testset "zero-sum for multiple param values" begin
for si in sum(β̂, 1)
    @test si≈0.0 atol=1e-6
end
end

info("Optimize at multiple tuning parameter values (warm start)")
β̂ws, = lsq_constrsparsereg(X, y, ρlist; Aeq = Aeq, beq = beq,
    penwt = penwt, solver = solver, warmstart = true)
#@show sum(β̂ws, 1)
@testset "zero-sum for multiple param values" begin
for si in sum(β̂ws, 1)
    @test si≈0.0 atol=1e-6
end
end

end # end of module


text:
function map_color(color_type, bit_depth)
    if color_type == PNG_COLOR_TYPE_GRAY
        colors_type = Gray{bit_depth}
    elseif color_type == PNG_COLOR_TYPE_PALETTE
        colors_type = RGB{bit_depth}
    elseif color_type == PNG_COLOR_TYPE_RGB
        colors_type = RGB{bit_depth}
    elseif color_type == PNG_COLOR_TYPE_RGB_ALPHA
        colors_type = RGBA{bit_depth}
    elseif color_type == PNG_COLOR_TYPE_GRAY_ALPHA
        colors_type = GrayA{bit_depth}
    else
        error("Unknown color type: $color_type")
    end
    return colors_type
end

function readimage(filename::String, transforms::Int = 0)
    fp = open_png(filename)

    png_ptr = create_read_struct()
    info_ptr = create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)
    png_set_sig_bytes(png_ptr)

    png_read_png(png_ptr, info_ptr, transforms)

    width = png_get_image_width(png_ptr, info_ptr)
    height = png_get_image_height(png_ptr, info_ptr)
    color_type = png_get_color_type(png_ptr, info_ptr)
    bit_depth = png_get_bit_depth(png_ptr, info_ptr)
    num_channels = png_get_channels(png_ptr, info_ptr)

    if (color_type == PNG_COLOR_TYPE_PALETTE)
        error("The color type PNG_COLOR_TYPE_PALETTE is not currently supported.")
    end

    @debug "Read image info:
            width=$width,
            height=$height,
            color_type=$color_type,
            bit_depth=$bit_depth,
            num_channels=$num_channels"

    even_depth = ((bit_depth + 1) >> 1) << 1
    if bit_depth <= 8
        T = Normed{UInt8, 8}
    elseif bit_depth <= 16
        T = Normed{UInt16, even_depth}
    else
        # TODO UInt32?
        error("Unknown bit_depth: $bit_depth")
    end

    colors_type = map_color(color_type, T)

    buf = Array{colors_type}(height, width)

    get_image_pixels!(rawview(channelview(buf)), png_ptr, info_ptr)

    png_destroy_read_struct(png_ptr, info_ptr)

    close_png(fp)

    return buf
end

function get_image_pixels!{T<:Unsigned}(buf::AbstractArray{T, 2}, png_ptr::Ptr{Void}, info_ptr::Ptr{Void})
    height, width = size(buf)
    rows = ccall((:png_get_rows, libpng), Ptr{Ptr{T}}, (Ptr{Void}, Ptr{Void}), png_ptr, info_ptr)
    for i = 1:height
        row = unsafe_load(rows, i)
        for j = 1:width
            buf[i, j] = unsafe_load(row, j)
        end
    end
    buf
end

function get_image_pixels!{T<:Unsigned}(buf::AbstractArray{T, 3}, png_ptr::Ptr{Void}, info_ptr::Ptr{Void})
    num_channels, height, width = size(buf)
    rows = ccall((:png_get_rows, libpng), Ptr{Ptr{T}}, (Ptr{Void}, Ptr{Void}), png_ptr, info_ptr)
    for i = 1:height
        row = unsafe_load(rows, i)
        for j = 1:width
            for c = 1:num_channels
                buf[c, i, j] = unsafe_load(row, num_channels * (j - 1) + c)
            end
        end
    end
    buf
end

function get_image_pixels!{T, N}(buf::AbstractArray{T, N}, png_ptr::Ptr{Void}, info_ptr::Ptr{Void})
    error("Image array has invalid dimension $N")
end

to_raw{C<:Colorant}(A::Array{C})  = to_raw(channelview(A))
to_raw{T<:Normed}(A::Array{T})    = rawview(A)
to_raw{T<:Real}(A::Array{T})      = to_raw(convert(Array{N0f8}, A))
to_raw(A::ColorView) = channelview(A)

get_bit_depth{C<:Colorant}(img::AbstractArray{C}) = _get_bit_depth(eltype(C))
get_bit_depth{T<:Normed}(img::AbstractArray{T}) = _get_bit_depth(T)
_get_bit_depth{T, N}(::Type{Normed{T, N}}) = N
_get_bit_depth{T<:AbstractFloat}(img::Type{T}) = 8
_get_bit_depth(img::Type{Bool}) = 8

get_color_type{T}(::Type{Gray{T}})  = PNG_COLOR_TYPE_GRAY
get_color_type{T}(::Type{GrayA{T}}) = PNG_COLOR_TYPE_GRAY_ALPHA
get_color_type{T}(::Type{RGB{T}})   = PNG_COLOR_TYPE_RGB
get_color_type{T}(::Type{RGBA{T}})  = PNG_COLOR_TYPE_RGBA
get_color_type{T<:Normed}(::Type{T}) = PNG_COLOR_TYPE_RGB

map_image{T}(c::Gray{T}) = convert(Gray{N0f8}, c)
map_image{T<:Normed}(c::Gray{T}) = c
map_image{T}(c::GrayA{T}) = convert(GrayA{N0f8}, c)
map_image{T<:Normed}(c::GrayA{T}) = c
map_image{T}(c::RGB{T}) = convert(RGB{N0f8}, c)
map_image{T<:Normed}(c::RGB{T}) = c
map_image{T}(c::RGBA{T}) = convert(RGBA{N0f8}, c)
map_image{T<:Normed}(c::RGBA{T}) = c

map_image(x::Bool) = convert(Gray{N0f8}, x)
map_image(x::AbstractFloat) = convert(N0f8, x)
map_image(x::Normed) = x

get_image_size{T}(buffer::AbstractArray{T,2}) = size(buffer)
get_image_size{T,N}(buffer::AbstractArray{T,N}) = error("Number of dimensions in image of $ndims not supported.")
function get_image_size{T}(buffer::AbstractArray{T,3})
    n_channels, height, width = size(buffer)
    height, width
end

function writeimage{T}(filename::String, image::AbstractArray{T})

    fp = ccall((:fopen, "libc"), Ptr{Void}, (Cstring, Cstring), filename, "wb")
    fp == C_NULL && error("Could not open $(filename) for writing")

    png_ptr = png_create_write_struct(png_error_fn, png_warn_fn)
    info_ptr = png_create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)

    image = map(map_image, image)

    buffer = to_raw(image)

    height, width = get_image_size(buffer)
    bit_depth = get_bit_depth(image)

    color_type = get_color_type(eltype(image))
    interlace = 0        # Set to always off
    compression_type = 0 # Set to always off
    filter_type = 0      # Set to always off

    @debug "Write image info:
            width=$width,
            height=$height,
            bit_depth=$bit_depth,
            color_type=$color_type,
            interlace=$interlace,
            compression_type=$compression_type,
            filter_type=$filter_type"

    ccall((:png_set_IHDR, libpng), Void,
          (Ptr{Void}, Ptr{Void}, Cuint, Cuint, Cint, Cint, Cint, Cint, Cint),
          png_ptr, info_ptr, width, height, bit_depth, color_type, interlace,
          compression_type, filter_type)

    png_write_info(png_ptr, info_ptr)

    write_rows(buffer, png_ptr, info_ptr)

    png_destroy_write_struct(png_ptr, info_ptr)

    close_png(fp)

    return nothing
end

# 2 dim matrix
function write_rows{T}(buf::AbstractArray{T, 2}, png_ptr::Ptr{Void}, info_ptr::Ptr{Void})
    height, width = get_image_size(buf)
    for row = 1:height
        row_buf = buf[row, :]
        ccall((:png_write_row, libpng), Void, (Ptr{Void}, Ptr{T}), png_ptr, row_buf)
    end
    ccall((:png_write_end, libpng), Void, (Ptr{Void}, Ptr{Void}), png_ptr, info_ptr)
end

# 3-dim matrix
function write_rows{T}(buf::AbstractArray{T, 3}, png_ptr::Ptr{Void}, info_ptr::Ptr{Void})
    height, width = get_image_size(buf)
    for row = 1:height
        row_buf = buf[:, row, :]
        ccall((:png_write_row, libpng), Void, (Ptr{Void}, Ptr{T}), png_ptr, row_buf)
    end
    ccall((:png_write_end, libpng), Void, (Ptr{Void}, Ptr{Void}), png_ptr, info_ptr)
end

function write_rows{T, N}(buf::AbstractArray{T, N}, png_ptr::Ptr{Void}, info_ptr::Ptr{Void})
    error("Image has invalid dimension $N")
end


text:
# ---
# title: 593. Valid Square
# id: problem593
# author: Indigo
# date: 2021-06-26
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/valid-square/description/>
# hidden: true
# ---
#
# Given the coordinates of four points in 2D space `p1`, `p2`, `p3` and `p4`,
# return `true` _if  the four points construct a square_.
#
# The coordinate of a point `pi` is represented as `[xi, yi]`. The input is
# **not** given in any order.
#
# A **valid square** has four equal sides with positive length and four equal
# angles (90-degree angles).
#
#
#
# **Example 1:**
#
#
#
#     Input: p1 = [0,0], p2 = [1,1], p3 = [1,0], p4 = [0,1]
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: p1 = [0,0], p2 = [1,1], p3 = [1,0], p4 = [0,12]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: p1 = [1,0], p2 = [-1,0], p3 = [0,1], p4 = [0,-1]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `p1.length == p2.length == p3.length == p4.length == 2`
#   * `-104 <= xi, yi <= 104`
#
#
## @lc code=start
using LeetCode

function valid_square(ps::Vector{Vector{Int}})
    sort!(ps)
    ps[1] + ps[4] != ps[2] + ps[3] && return false
    sum((ps[1] .- ps[2]) .^ 2) != sum((ps[1] .- ps[3]) .^ 2) && return false
    return true
end
## @lc code=end


text:
using PrettyPrint

PrettyPrint.pprint_impl(io, sym :: Symbol, indent, newline) = print(io, repr(sym))

PrettyPrint.pprint_impl(io, expr :: Expr, indent, newline) = begin
    print(io, "Expr(\n")

    print(io, repeat(" ", indent + 2))
    print(io, repr(expr.head))
    print(io, ",\n")

    for each in expr.args
        print(io, repeat(" ", indent + 2))
        pprint(io, each, indent + 2, false)
        print(io, ",\n")
    end

    print(io, repeat(" ", indent))
    print(io, ")")
end

PrettyPrint.pprint_impl(io, expr :: GlobalRef, indent, newline) = begin
    print(io, "GlobalRef(\n")

    for each in [expr.mod, expr.name]
        print(io, repeat(" ", indent + 2))
        pprint(io, each, indent + 2, false)
        print(io, ",\n")
    end

    print(io, repeat(" ", indent))
    print(io, ")")
end

PrettyPrint.pprint_impl(io, lnode :: LineNumberNode, indent, newline) = begin
    print(io, "LineNumberNode(\n")

    for each in [lnode.line, lnode.file]
        print(io, repeat(" ", indent + 2))
        pprint(io, each, indent + 2, false)
        print(io, ",\n")
    end

    print(io, repeat(" ", indent))
    print(io, ")")
end

PrettyPrint.pprint_impl(io, qnode :: QuoteNode, indent, newline) = begin
    print(io, "QuoteNode(\n")

    print(io, repeat(" ", indent + 2))
    pprint(io, qnode.value, indent + 2, false)
    print(io, ",\n")

    print(io, repeat(" ", indent))
    print(io, ")")
end

PrettyPrint.pprint_impl(io, mod :: Union{Module, Type}, indent, newline) = begin
    print(io, mod)
end


text:
using Base.Test
using Distributions

n = 100
w = rand(100)

# auxiliary tools

function ccount(K, x)
	r = zeros(K)
	for i = 1:length(x)
		r[x[i]] += 1.0
	end
	r
end

function ccount(K, x, w)
	r = zeros(K)
	for i = 1:length(x)
		r[x[i]] += w[i]
	end
	r
end

# Beta - Bernoulli / Binomial

pri = Beta(1.0, 2.0)

x = rand(Bernoulli(0.3), n)
p = posterior(pri, Bernoulli, x)
@test isa(p, Beta)
@test_approx_eq p.alpha pri.alpha + sum(x)
@test_approx_eq p.beta  pri.beta + (n - sum(x))

p = posterior(pri, Bernoulli, x, w)
@test isa(p, Beta)
@test_approx_eq p.alpha pri.alpha + sum(x .* w)
@test_approx_eq p.beta  pri.beta + (sum(w) - sum(x .* w))

x = rand(Binomial(10, 0.3), n)
p = posterior(pri, Binomial, 10, x)
@test isa(p, Beta)
@test_approx_eq p.alpha pri.alpha + sum(x)
@test_approx_eq p.beta  pri.beta + (10n - sum(x))

p = posterior(pri, Binomial, 10, x, w)
@test isa(p, Beta)
@test_approx_eq p.alpha pri.alpha + sum(x .* w)
@test_approx_eq p.beta  pri.beta + (10 * sum(w) - sum(x .* w))

# Dirichlet - Categorical / Multinomial

pri = Dirichlet([1., 2., 3.])

x = rand(Categorical([0.2, 0.3, 0.5]), n)
p = posterior(pri, Categorical, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + ccount(3, x)

p = posterior(pri, Categorical, x, w)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + ccount(3, x, w)

x = rand(Multinomial(100, [0.2, 0.3, 0.5]), 1)
p = posterior(pri, Multinomial, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + x

x = rand(Multinomial(10, [0.2, 0.3, 0.5]), n)
p = posterior(pri, Multinomial, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + vec(sum(x, 2))

p = posterior(pri, Multinomial, x, w)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + vec(x * w)

# Gamma - Exponential

pri = Gamma(1.5, 2.0)

x = rand(Exponential(2.0), n)
p = posterior(pri, Exponential, x)
@test isa(p, Gamma)
@test_approx_eq p.shape pri.shape + n
@test_approx_eq rate(p) rate(pri) + sum(x)

p = posterior(pri, Exponential, x, w)
@test isa(p, Gamma)
@test_approx_eq p.shape pri.shape + sum(w)
@test_approx_eq rate(p) rate(pri) + sum(x .* w)

# Normal likelihood

pri = Normal(1.0, 5.0)

x = rand(Normal(2.0, 3.0), n)
p = posterior((pri, 3.0), Normal, x)
@test isa(p, Normal)
@test_approx_eq mean(p)  (mean(pri) / var(pri) + sum(x) / 9.0) / (1.0 / var(pri) + n / 9.0)
@test_approx_eq var(p) inv(1.0 / var(pri) + n / 9.0)

p = posterior((pri, 3.0), Normal, x, w)
@test isa(p, Normal)
@test_approx_eq mean(p)  (mean(pri) / var(pri) + dot(x, w) / 9.0) / (1.0 / var(pri) + sum(w) / 9.0)
@test_approx_eq var(p) inv(1.0 / var(pri) + sum(w) / 9.0)

pri = InvertedGamma(1.5, 0.5) # β = 2.0

x = rand(Normal(2.0, 3.0), n)
p = posterior((2.0, pri), Normal, x)
@test isa(p, InvertedGamma)
@test_approx_eq p.shape pri.shape + n / 2
@test_approx_eq rate(p) rate(pri) + sum(abs2(x - 2.0)) / 2

p = posterior((2.0, pri), Normal, x, w)
@test isa(p, InvertedGamma)
@test_approx_eq p.shape pri.shape + sum(w) / 2
@test_approx_eq rate(p) rate(pri) + dot(w, abs2(x - 2.0)) / 2


text:
# ---
# title: 78. Subsets
# id: problem78
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Backtracking, Bit Manipulation
# link: <https://leetcode.com/problems/subsets/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, return _all possible subsets (the power set)_.
#
# The solution set must not contain duplicate subsets.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3]
#     Output: [[],[1],[2],[1,2],[3],[1,3],[2,3],[1,2,3]]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [0]
#     Output: [[],[0]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10`
#   * `-10 <= nums[i] <= 10`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# crisis.jl
#######################

# Load packages.
using Plots, CSV, Dates, DataFrames

# Download data from the MHLW web site.
download("https://covid19.mhlw.go.jp/public/opendata/confirmed_cases_cumulative_daily.csv","./csv/mhlw_cases.csv");
download("https://covid19.mhlw.go.jp/public/opendata/deaths_cumulative_daily.csv","./csv/mhlw_deaths.csv");
download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv","./csv/jhu_cases.csv");
download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv","./csv/jhu_deaths.csv");

# Shape up the data.
Ccsv=DataFrame(CSV.File("./csv/mhlw_cases.csv", header=false, delim=','));
Dcsv=DataFrame(CSV.File("./csv/mhlw_deaths.csv", header=false, delim=','));
Wcsv=DataFrame(CSV.File("./csv/jhu_cases.csv", header=false, delim=','));
Xcsv=DataFrame(CSV.File("./csv/jhu_deaths.csv", header=false, delim=','));
(pa,qa)=size(Ccsv);
l=Int64(floor((pa-1)/48));
(pw,qw)=size(Wcsv);

# Plot the data
# d0: the initial date
# df: the final day
dd0=Date(2021,12,1)
DD=Int64(qw-683);
dd1=dd0+Day(floor((DD-1)/3));
dd2=dd0+Day(floor(2*(DD-1)/3));
dd3=dd0+Day(DD-1);
ll0=string(dd0);
ll1=string(dd1);
ll2=string(dd2);
ll3=string(dd3);

J=parse.(Float64,Ccsv[2:pa,2:qa]);
K=parse.(Float64,Dcsv[2:pa,2:qa]);

# Okinawa
POKNW=1.458870;
COKNW=J[:,48]/POKNW
NOKNW=zeros(DD);
for j=1:DD
    NOKNW[j]=(COKNW[j+571]-COKNW[j+564])/7
end
DOKNW=K[:,48]/POKNW
NDOKNW=zeros(DD);
for j=1:DD
    NDOKNW[j]=(DOKNW[j+571]-DOKNW[j+564])/7
end

# Tokyo
PTKY=14.049146;
CTKY=J[:,14]/PTKY
NTKY=zeros(DD);
for j=1:DD
    NTKY[j]=(CTKY[j+571]-CTKY[j+564])/7
end

# Argentina
PARG=45.670451;
BARG=parse.(Float64,Array(Xcsv[8,5:qw]))/PARG;
NDARG=zeros(DD);
for j=1:DD
    NDARG[j]=(BARG[j+679]-BARG[j+672])/7
end

# NSW
PNSW=8.196;
ANSW=parse.(Float64,Array(Wcsv[11,5:qw]))/PNSW;
NNSW=zeros(DD);
for j=1:DD
    NNSW[j]=(ANSW[679+j]-ANSW[672+j])/7
end

# Victoria
PVIC=6.7;
AVIC=parse.(Float64,Array(Wcsv[16,5:qw]))/PVIC;
NVIC=zeros(DD);
for j=1:DD
    NVIC[j]=(AVIC[679+j]-AVIC[672+j])/7
end

# Brazil
PBRA=214.289417;
BBRA=parse.(Float64,Array(Xcsv[32,5:qw]))/PBRA;
NDBRA=zeros(DD);
for j=1:DD
    NDBRA[j]=(BBRA[j+679]-BBRA[j+672])/7
end

# Brunei Darussalam
PBWN=0.442205;
ABWN=parse.(Float64,Array(Wcsv[33,5:qw]))/PBWN;
NBWN=zeros(DD);
for j=1:DD
    NBWN[j]=(ABWN[679+j]-ABWN[672+j])/7
end

# Colombia
PCOL=51.503463;
BCOL=parse.(Float64,Array(Xcsv[94,5:qw]))/PCOL;
NDCOL=zeros(DD);
for j=1:DD
    NDCOL[j]=(BCOL[j+679]-BCOL[j+672])/7
end

# Indonesia
PIDN=276.833206;
AIDN=parse.(Float64,Array(Wcsv[150,5:qw]))/PIDN;
NIDN=zeros(DD);
for j=1:DD
    NIDN[j]=(AIDN[679+j]-AIDN[672+j])/7
end
BIDN=parse.(Float64,Array(Xcsv[150,5:qw]))/PIDN;
NDIDN=zeros(DD);
for j=1:DD
    NDIDN[j]=(BIDN[j+679]-BIDN[j+672])/7
end

# Israel
PISR=8.789774;
AISR=parse.(Float64,Array(Wcsv[154,5:qw]))/PISR;
NISR=zeros(DD);
for j=1:DD
    NISR[j]=(AISR[679+j]-AISR[672+j])/7
end
BISR=parse.(Float64,Array(Xcsv[154,5:qw]))/PISR;
NDISR=zeros(DD);
for j=1:DD
    NDISR[j]=(BISR[j+679]-BISR[j+672])/7
end

# Japan
PJPN=125.36;
AJPN=parse.(Float64,Array(Wcsv[157,5:qw]))/PJPN;
NJPN=zeros(DD);
for j=1:DD
    NJPN[j]=(AJPN[679+j]-AJPN[672+j])/7
end

# South Korea
PKOR=51.318552;
AKOR=parse.(Float64,Array(Wcsv[162,5:qw]))/PKOR;
NKOR=zeros(DD);
for j=1:DD
    NKOR[j]=(AKOR[679+j]-AKOR[672+j])/7
end

# Malaysia
PMYS=32.66;
AMYS=parse.(Float64,Array(Wcsv[178,5:qw]))/PMYS;
NMYS=zeros(DD);
for j=1:DD
    NMYS[j]=(AMYS[679+j]-AMYS[672+j])/7
end

# Mexico
PMEX=130.482814;
BMEX=parse.(Float64,Array(Xcsv[185,5:qw]))/PMEX;
NDMEX=zeros(DD);
for j=1:DD
    NDMEX[j]=(BMEX[j+679]-BMEX[j+672])/7
end

# Peru
PPER=33.510887;
APER=parse.(Float64,Array(Wcsv[213,5:qw]))/PPER;
NPER=zeros(DD);
for j=1:DD
    NPER[j]=(APER[679+j]-APER[672+j])/7
end
BMYS=parse.(Float64,Array(Xcsv[178,5:qw]))/PMYS;
NDMYS=zeros(DD);
for j=1:DD
    NDMYS[j]=(BMYS[j+679]-BMYS[j+672])/7
end
BPER=parse.(Float64,Array(Xcsv[213,5:qw]))/PPER;
NDPER=zeros(DD);
for j=1:DD
    NDPER[j]=(BPER[j+679]-BPER[j+672])/7
end

# Philippines
PPHI=111.249116;
APHI=parse.(Float64,Array(Wcsv[214,5:qw]))/PPHI;
NPHI=zeros(DD);
for j=1:DD
    NPHI[j]=(APHI[679+j]-APHI[672+j])/7
end
BPHI=parse.(Float64,Array(Xcsv[214,5:qw]))/PPHI;
NDPHI=zeros(DD);
for j=1:DD
    NDPHI[j]=(BPHI[j+679]-BPHI[j+672])/7
end

# Russia
PRUS=146.013169;
ARUS=parse.(Float64,Array(Wcsv[219,5:qw]))/PRUS;
NRUS=zeros(DD);
for j=1:DD
    NRUS[j]=(ARUS[679+j]-ARUS[672+j])/7
end
BRUS=parse.(Float64,Array(Xcsv[219,5:qw]))/PRUS;
NDRUS=zeros(DD);
for j=1:DD
    NDRUS[j]=(BRUS[j+679]-BRUS[j+672])/7
end

# Singapore
PSIN=5.902011;
ASIN=parse.(Float64,Array(Wcsv[232,5:qw]))/PSIN;
NSIN=zeros(DD);
for j=1:DD
    NSIN[j]=(ASIN[679+j]-ASIN[672+j])/7
end
BSIN=parse.(Float64,Array(Xcsv[232,5:qw]))/PSIN;
NDSIN=zeros(DD);
for j=1:DD
    NDSIN[j]=(BSIN[j+679]-BSIN[j+672])/7
end

# Sri Lanka
PLKA=21.516097;
ALKA=parse.(Float64,Array(Wcsv[240,5:qw]))/PLKA;
NLKA=zeros(DD);
for j=1:DD
    NLKA[j]=(ALKA[679+j]-ALKA[672+j])/7
end
BLKA=parse.(Float64,Array(Xcsv[240,5:qw]))/PLKA;
NDLKA=zeros(DD);
for j=1:DD
    NDLKA[j]=(BLKA[j+679]-BLKA[j+672])/7
end

# Thailand
PTHA=70.000662;
ATHA=parse.(Float64,Array(Wcsv[250,5:qw]))/PTHA;
NTHA=zeros(DD);
for j=1:DD
    NTHA[j]=(ATHA[679+j]-ATHA[672+j])/7
end
BTHA=parse.(Float64,Array(Xcsv[250,5:qw]))/PTHA;
NDTHA=zeros(DD);
for j=1:DD
    NDTHA[j]=(BTHA[j+679]-BTHA[j+672])/7
end

# United Staes
PUSA=333.225477;
AUSA=parse.(Float64,Array(Wcsv[257,5:qw]))/PUSA;
NUSA=zeros(DD);
for j=1:DD
    NUSA[j]=(AUSA[679+j]-AUSA[672+j])/7
end
BUSA=parse.(Float64,Array(Xcsv[257,5:qw]))/PUSA;
NDUSA=zeros(DD);
for j=1:DD
    NDUSA[j]=(BUSA[j+679]-BUSA[j+672])/7
end

# United Kingdom
PGBR=68.294438;
AGBR=parse.(Float64,Array(Wcsv[272,5:qw]))/PGBR;
NGBR=zeros(DD);
for j=1:DD
    NGBR[j]=(max(AGBR[679+j]-AGBR[672+j],0))/7
end
BGBR=parse.(Float64,Array(Xcsv[272,5:qw]))/PGBR;
NDGBR=zeros(DD);
for j=1:DD
    NDGBR[j]=(BGBR[j+679]-BGBR[j+672])/7
end

# Vietnam
PVNM=98.341025;
AVNM=parse.(Float64,Array(Wcsv[277,5:qw]))/PVNM;
NVNM=zeros(DD);
for j=1:DD
    NVNM[j]=(AVNM[679+j]-AVNM[672+j])/7
end

p=plot([NKOR NMYS NOKNW NSIN NNSW NVIC NRUS NVNM NGBR NUSA],
    grid=false,
    linewidth=2,
    title="COVID-19 7-day average of daily new cases per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    xticks = ([1 floor(DD/3) floor(2*DD/3) DD;], [ll0 ll1 ll2 ll3]),
    xlabel="date",
    yaxis="cases/1M",
    legendfont=font(8),
    label=["South Korea" "Malaysia" "Okinawa" "Singapore" "New South Wales" "Victoria" "Russia" "Vietnam" "United Kingdom" "United States"],
    palette = :seaborn_bright,
    legend = :topleft)
savefig("./crisis/cases.png")

q=plot([NJPN NTKY NPHI NBWN NTHA NIDN],
    grid=false,
    linewidth=2,
    title="COVID-19 7-day average of daily new cases per 1M \n data sourced by JHU and MOH of Japan",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    xticks = ([1 floor(DD/3) floor(2*DD/3) DD;], [ll0 ll1 ll2 ll3]),
    xlabel="date",
    yaxis="cases/1M",
    legendfont=font(8),
    label=["Japan" "Tokyo" "Philippines" "Brunei Darussalam" "Thailand" "Indonesia"],
    palette = :seaborn_bright,
    legend = :topright)
savefig("./crisis/delta.png")

r=plot([NDMYS NDBRA NDOKNW NDSIN NDMEX NDLKA NDRUS NDTHA NDGBR NDUSA],
    grid=false,
    linewidth=2,
    title="COVID-19: 7-day average deaths per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(DD/3) floor(2*DD/3) DD;], [ll0 ll1 ll2 ll3]),
    xlabel="date",
    yaxis="deaths/1M",
    legendfont=font(8),
    label=["Malaysia" "Brazil" "Okinawa" "Singapore" "Mexico" "Sri Lanka" "Russia" "Thailand" "United Kingdam" "United States"],
    palette = :seaborn_bright,
    legend = :topleft)
savefig("./crisis/deaths.png")

plot(p, q, r,
     layout=3,
     size=(1260,840),
     left_margin=Plots.Measures.Length(:mm, 5.0),
     right_margin=Plots.Measures.Length(:mm, 15.0),
     top_margin=Plots.Measures.Length(:mm, 5.0),
     bottom_margin=Plots.Measures.Length(:mm, 5.0))
savefig("./crisis/crisis.png")


text:
function orderstatistic(Xs::IIDRandomSequence{Continuous,Uniform},k::Int64)
  if params(Xs.d) == (0.,1.)
    return Beta(k,length(Xs)+1-k)
  else
    error("location/scale aspect not yet implemented")
  end
end

function range(Xs::IIDRandomSequence{Continuous,Uniform})
  if params(Xs.d) == (0.,1.)
    if length(Xs)>1
      return Beta(length(Xs)-1,2)
    else
      error("sequence must have length greater than 1 for range to exist")
    end
  else
    error("location/scale aspect not yet implemented")
  end
end

function spacing(Xs::IIDRandomSequence{Continuous,Uniform},space::Integer)
  if params(Xs.d) == (0.,1.)
    if length(Xs)>1
      return Beta(1,length(Xs))
    else
      error("sequence must have length greater than 1 for range to exist")
    end
  else
    error("location/scale aspect not yet implemented")
  end
end

function jointspacing(Xs::IIDRandomSequence{Continuous,Uniform})
  if params(Xs.d) == (0.,1.)
    if length(Xs)>1
      return Dirichlet(length(Xs)+1,1)
    else
      error("sequence must have length greater than 1 for range to exist")
    end
  else
    error("location/scale aspect not yet implemented")
  end
end


text:
# test format spec parsing

using Formatting
using Test

# default spec
fs = FormatSpec("")
@test fs.typ == 's'
@test fs.fill == ' '
@test fs.align == '<'
@test fs.sign == '-'
@test fs.width == -1
@test fs.prec == -1
@test fs.ipre == false
@test fs.zpad == false
@test fs.tsep == false

# more cases

fs = FormatSpec("d")
@test fs == FormatSpec('d')
@test fs.align == '>'

@test FormatSpec("8x") == FormatSpec('x'; width=8)
@test FormatSpec("08b") == FormatSpec('b'; width=8, zpad=true)
@test FormatSpec("12f") == FormatSpec('f'; width=12, prec=6)
@test FormatSpec("12.7f") == FormatSpec('f'; width=12, prec=7)
@test FormatSpec("+08o") == FormatSpec('o'; width=8, zpad=true, sign='+')

@test FormatSpec("8") == FormatSpec('s'; width=8)
@test FormatSpec(".6f") == FormatSpec('f'; prec=6)
@test FormatSpec("<8d") == FormatSpec('d'; width=8, align='<')
@test FormatSpec("#<8d") == FormatSpec('d'; width=8, fill='#', align='<')
@test FormatSpec("⋆<8d") == FormatSpec('d'; width=8, fill='⋆', align='<')
@test FormatSpec("#8,d") == FormatSpec('d'; width=8, ipre=true, tsep=true)

# format string

@test fmt("", "abc") == "abc"
@test fmt("", "αβγ") == "αβγ"
@test fmt("s", "abc") == "abc"
@test fmt("s", "αβγ") == "αβγ"
@test fmt("2s", "abc") == "abc"
@test fmt("2s", "αβγ") == "αβγ"
@test fmt("5s", "abc") == "abc  "
@test fmt("5s", "αβγ") == "αβγ  "
@test fmt(">5s", "abc") == "  abc"
@test fmt(">5s", "αβγ") == "  αβγ"
@test fmt("*>5s", "abc") == "**abc"
@test fmt("⋆>5s", "αβγ") == "⋆⋆αβγ"
@test fmt("*<5s", "abc") == "abc**"
@test fmt("⋆<5s", "αβγ") == "αβγ⋆⋆"

# format char

@test fmt("", 'c') == "c"
@test fmt("", 'γ') == "γ"
@test fmt("c", 'c') == "c"
@test fmt("c", 'γ') == "γ"
@test fmt("3c", 'c') == "c  "
@test fmt("3c", 'γ') == "γ  "
@test fmt(">3c", 'c') == "  c"
@test fmt(">3c", 'γ') == "  γ"
@test fmt("*>3c", 'c') == "**c"
@test fmt("⋆>3c", 'γ') == "⋆⋆γ"
@test fmt("*<3c", 'c') == "c**"
@test fmt("⋆<3c", 'γ') == "γ⋆⋆"

# format integer

@test fmt("", 1234) == "1234"
@test fmt("d", 1234) == "1234"
@test fmt("n", 1234) == "1234"
@test fmt("x", 0x2ab) == "2ab"
@test fmt("X", 0x2ab) == "2AB"
@test fmt("o", 0o123) == "123"
@test fmt("b", 0b1101) == "1101"

@test fmt("d", 0) == "0"
@test fmt("d", 9) == "9"
@test fmt("d", 10) == "10"
@test fmt("d", 99) == "99"
@test fmt("d", 100) == "100"
@test fmt("d", 1000) == "1000"

@test fmt("06d", 123) == "000123"
@test fmt("+6d", 123) == "  +123"
@test fmt("+06d", 123) == "+00123"
@test fmt(" d", 123) == " 123"
@test fmt(" 6d", 123) == "   123"
@test fmt("<6d", 123) == "123   "
@test fmt(">6d", 123) == "   123"
@test fmt("*<6d", 123) == "123***"
@test fmt("⋆<6d", 123) == "123⋆⋆⋆"
@test fmt("*>6d", 123) == "***123"
@test fmt("⋆>6d", 123) == "⋆⋆⋆123"
@test fmt("< 6d", 123) == " 123  "
@test fmt("<+6d", 123) == "+123  "
@test fmt("> 6d", 123) == "   123"
@test fmt(">+6d", 123) == "  +123"

@test fmt("+d", -123) == "-123"
@test fmt("-d", -123) == "-123"
@test fmt(" d", -123) == "-123"
@test fmt("06d", -123) == "-00123"
@test fmt("<6d", -123) == "-123  "
@test fmt(">6d", -123) == "  -123"

# format floating point (f)

@test fmt("", 0.125) == "0.125"
@test fmt("f", 0.0) == "0.000000"
@test fmt("f", -0.0) == "-0.000000"
@test fmt("f", 0.001) == "0.001000"
@test fmt("f", 0.125) == "0.125000"
@test fmt("f", 1.0/3) == "0.333333"
@test fmt("f", 1.0/6) == "0.166667"
@test fmt("f", -0.125) == "-0.125000"
@test fmt("f", -1.0/3) == "-0.333333"
@test fmt("f", -1.0/6) == "-0.166667"
@test fmt("f", 1234.5678) == "1234.567800"
@test fmt("8f", 1234.5678) == "1234.567800"

@test fmt("8.2f", 8.376) == "    8.38"
@test fmt("<8.2f", 8.376) == "8.38    "
@test fmt(">8.2f", 8.376) == "    8.38"
@test fmt("8.2f", -8.376) == "   -8.38"
@test fmt("<8.2f", -8.376) == "-8.38   "
@test fmt(">8.2f", -8.376) == "   -8.38"

@test fmt("<08.2f", 8.376) == "00008.38"
@test fmt(">08.2f", 8.376) == "00008.38"
@test fmt("<08.2f", -8.376) == "-0008.38"
@test fmt(">08.2f", -8.376) == "-0008.38"
@test fmt("*<8.2f", 8.376) == "8.38****"
@test fmt("⋆<8.2f", 8.376) == "8.38⋆⋆⋆⋆"
@test fmt("*>8.2f", 8.376) == "****8.38"
@test fmt("⋆>8.2f", 8.376) == "⋆⋆⋆⋆8.38"
@test fmt("*<8.2f", -8.376) == "-8.38***"
@test fmt("⋆<8.2f", -8.376) == "-8.38⋆⋆⋆"
@test fmt("*>8.2f", -8.376) == "***-8.38"
@test fmt("⋆>8.2f", -8.376) == "⋆⋆⋆-8.38"

@test fmt(".2f", 0.999) == "1.00"
@test fmt(".2f", 0.996) == "1.00"
@test fmt("6.2f", 9.999) == " 10.00"
# Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
@test (fmt(".2f", 0.995) == "1.00" || fmt(".2f", 0.995) == "0.99")
@test fmt(".2f", 0.994) == "0.99"

# format floating point (e)

@test fmt("E", 0.0) == "0.000000E+00"
@test fmt("e", 0.0) == "0.000000e+00"
@test fmt("e", 0.001) == "1.000000e-03"
@test fmt("e", 0.125) == "1.250000e-01"
@test fmt("e", 100/3) == "3.333333e+01"
@test fmt("e", -0.125) == "-1.250000e-01"
@test fmt("e", -100/6) == "-1.666667e+01"
@test fmt("e", 1234.5678) == "1.234568e+03"
@test fmt("8e", 1234.5678) == "1.234568e+03"

@test fmt("<12.2e", 13.89) == "1.39e+01    "
@test fmt(">12.2e", 13.89) == "    1.39e+01"
@test fmt("*<12.2e", 13.89) == "1.39e+01****"
@test fmt("⋆<12.2e", 13.89) == "1.39e+01⋆⋆⋆⋆"
@test fmt("*>12.2e", 13.89) == "****1.39e+01"
@test fmt("⋆>12.2e", 13.89) == "⋆⋆⋆⋆1.39e+01"
@test fmt("012.2e", 13.89) == "00001.39e+01"
@test fmt("012.2e", -13.89) == "-0001.39e+01"
@test fmt("+012.2e", 13.89) == "+0001.39e+01"

@test fmt(".1e", 0.999) == "1.0e+00"
@test fmt(".1e", 0.996) == "1.0e+00"
# Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
@test (fmt(".1e", 0.995) == "1.0e+00" || fmt(".1e", 0.995) == "9.9e-01")
@test fmt(".1e", 0.994) == "9.9e-01"
@test fmt(".1e", 0.6) == "6.0e-01"
@test fmt(".1e", 0.9) == "9.0e-01"

# issue #61
@test fmt("1.0e", 1e-21) == "1.e-21"
@test fmt("1.1e", 1e-21) == "1.0e-21"

@test fmt("10.2e", 1.2e100) == " 1.20e+100"
@test fmt("11.2e", BigFloat("1.2e1000")) == " 1.20e+1000"
@test fmt("11.2e", BigFloat("1.2e-1000")) == " 1.20e-1000"
@test fmt("9.2e", 9.999e9) == " 1.00e+10"
@test fmt("10.2e", 9.999e99) == " 1.00e+100"
@test fmt("11.2e", BigFloat("9.999e999")) == " 1.00e+1000"
@test fmt("10.2e", -9.999e-100) == " -1.00e-99"

# format special floating point value

@test fmt("f", NaN) == "NaN"
@test fmt("e", NaN) == "NaN"
@test fmt("f", NaN32) == "NaN"
@test fmt("e", NaN32) == "NaN"

@test fmt("f", Inf) == "Inf"
@test fmt("e", Inf) == "Inf"
@test fmt("f", Inf32) == "Inf"
@test fmt("e", Inf32) == "Inf"

@test fmt("f", -Inf) == "-Inf"
@test fmt("e", -Inf) == "-Inf"
@test fmt("f", -Inf32) == "-Inf"
@test fmt("e", -Inf32) == "-Inf"

@test fmt("<5f", Inf) == "Inf  "
@test fmt(">5f", Inf) == "  Inf"
@test fmt("*<5f", Inf) == "Inf**"
@test fmt("⋆<5f", Inf) == "Inf⋆⋆"
@test fmt("*>5f", Inf) == "**Inf"
@test fmt("⋆>5f", Inf) == "⋆⋆Inf"


text:
"""
Tests for `estimate_h2`.
Tests all fail because `ld_score_regression` is only half-implemented.
"""

include("../src/LDScoreJulia.jl")

using Main
using Main.LDScoreJulia
using Test

approx(x, y, eps) = abs(x - y) <= eps

function vector_approx(x, y, eps)
    if size(x) != size(y)
        return false
    end
    for (i, j) in enumerate(x)
        if !approx(x[i], y[i], eps)
            return false
        end
    end
    return true
end

# Test case from original ldsc repo
x = LDScoreJulia.estimate_h2(
    Dict([
        ("ref_id", "LDScoreJulia/test/test_ldscore/oneld_onefile"),
        ("w_ld", "LDScoreJulia/test/test_ldscore/w"),
        ("h2", "LDScoreJulia/test/test_sumstats/1"),
        ("out", "LDScoreJulia/test/test_out/1"),

        # default parser arguments
        ("n_blocks", 200),
        ("χ²_max", nothing),
        ("two_step", nothing),
        ("intercept_h2", nothing),
        ("overlap_annot", false),
    ])
)

# TODO finish implementing `ld_score_regression`


text:
using LIBSVM
using ReservoirComputing
using Plots
using ParameterizedFunctions
using DifferentialEquations
using DynamicalSystems
using Statistics

function nrmse(y, yt, sigma)
    rmse = 0.0
    for i=1:size(y, 1)
        rmse += (y[i]-yt[i])^2.0
    end
    rmse = sqrt(rmse/(size(y, 1)*sigma))
    return rmse
end

function data_prep(data, shift, train_len, test_len, h)

    new_d = Matrix(embed(Dataset(data), 4, 6))'
    y_target = data[h:end]
    train_in, train_out = new_d[:, shift:shift+train_len-1], y_target[shift:shift+train_len-1]
    test_in, test_out = new_d[:, shift+train_len:shift+train_len+test_len-1], y_target[shift+train_len:shift+train_len+test_len-1]

    return train_in, train_out, test_in, test_out
end

function f_mackey_glass(u, h, p, t)
  z = h(p, t - 17)

  0.2 * z / (1 + z^10) - 0.1 * u
end

function h_mackey_glass(p, t)
  t ≤ 0 || error("history function is only implemented for t ≤ 0")

  0.5
end

prob = DDEProblem(f_mackey_glass, h_mackey_glass, (0.0, 700.0); constant_lags = [17])
sol = solve(prob, RK4(), adaptive=false,dt=0.1)
v = sol.u
data = Matrix(hcat(v...)')

#data parameters
const shift = 200
const train_len = 1000
const test_len = 500
const h = 84

#model parameters
const approx_res_size = 700
const sparsity = 0.02
const activation = tanh
const radius = 0.98
const sigma = 0.25

const alpha = 1.0
const nla_type = NLADefault()
const extended_states = true

#creating data
train_in, train_out, test_in, test_out = data_prep(data, shift, train_len, test_len, h)

println("Training and testing on noiseless data...")

#create echo state network
W = init_reservoir_givensp(approx_res_size, radius, sparsity)
W_in = init_dense_input_layer(approx_res_size, size(train_in, 1), sigma)
esn = ESN(W, train_in, W_in, activation, alpha, nla_type, extended_states)

#train and predict using svesm
m = SVESMtrain(EpsilonSVR(kernel = Kernel.Linear), esn, train_out)
output = SVESM_direct_predict(esn, test_in, m)
#computing nrmse normalized on the variance of the original time series
sig = var(data)
println("SVESM nmrse for the noiseless test: ", nrmse(test_out, output, sig))
svesm_nn = plot(output, label="SVESM")
svesm_nn = plot!(test_out, label="actual")
savefig(svesm_nn, "svesm_noiseless_comparison")

#svm training
m1 = fit!(EpsilonSVR(kernel = Kernel.Polynomial), train_in', train_out)
t1 = predict(m1, test_in')

m2 = fit!(EpsilonSVR(kernel = Kernel.RadialBasis), train_in', train_out)
t2 = predict(m2, test_in')

println("SVM Poly kernel nmrse for the noiseless test: ", nrmse(test_out, t1, sig))
println("SVM RadialBasis kernel nmrse for the noiseless test: ", nrmse(test_out, t2, sig))
svm_nn = plot(t1, label="SVM Polynomial kernel")
svm_nn = plot!(t2, label="SVM RadialBasis kernel")
svm_nn = plot!(test_out, label="actual")
savefig(svm_nn, "svm_noiseless_comparison")


text:
@testset "Testing MBTR Representation with dc Si" begin
@info("Testing MBTR Representation for dc Si.")
using DescriptorZoo, JuLIP, Test

at = bulk(:Si, cubic=true)
desc = mbtr(at)
mbtr_ref = [
       6.9496778788414781e-01,
       7.1986220284027114e-01,
       7.3354857103049786e-01,
       7.3666491362136577e-01,
       7.3044762991316925e-01,
       7.1649441177128370e-01,
       6.9649901670096226e-01,
       6.7200606559977982e-01,
       6.4422928693791326e-01,
       6.1396345774896632e-01,
       5.8160117982024029e-01,
       5.4724434055754667e-01,
       5.1088099007051713e-01,
       4.7258543971849859e-01,
       4.3269540091997261e-01,
       3.9192569929718940e-01,
       3.5139206883965507e-01,
       3.1253747322315001e-01,
       2.7697310097245148e-01,
       2.4626260205440162e-01,
       2.2168837199196242e-01,
       2.0404146122179259e-01,
       1.9347226369615067e-01,
       1.8942888993870177e-01,
       1.9069597290259666e-01,
       1.9553067571537033e-01,
       2.0187700957831695e-01,
       2.0762645100533930e-01,
       2.1088450453095892e-01,
       2.1020118590035677e-01,
       2.0472935788111615e-01,
       1.9428785263606121e-01,
       1.7932399339294011e-01,
       1.6078869816994643e-01,
       1.3995256578723689e-01,
       1.1819976052859571e-01,
       9.6836522223254493e-02,
       7.6943433911371095e-02,
       5.9287846045598072e-02,
       4.4298754123766174e-02,
       3.2094480531120281e-02,
       2.2546050707531993e-02,
       1.5356935355672306e-02,
       1.0142102152565563e-02,
       6.4943824952075910e-03,
       4.0321099048365632e-03,
       2.4272265052827208e-03,
       1.4166771385286556e-03,
       8.0170343488886913e-04,
       4.3988429175749678e-04,
       2.3401514207793466e-04,
       1.2070688477923311e-04,
       6.0367240455477009e-05,
       2.9271955306051653e-05,
       1.3762052062538477e-05,
       6.2732939120686559e-06,
       2.7726111754441095e-06,
       1.1881276847952497e-06,
       4.9364868828203077e-07,
       1.9886277542685130e-07,
       7.7672905494426746e-08,
       2.9414810799660235e-08,
       1.0800476431045606e-08,
       3.8450325381032670e-09,
       1.3272030624847988e-09,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       1.2398167472199218e-04,
       2.0466965823382139e-04,
       3.3061207241311385e-04,
       5.2173848485611489e-04,
       8.0094009760583779e-04,
       1.1899435089817198e-03,
       1.7035053236855724e-03,
       2.3441813321358261e-03,
       3.1009483013927528e-03,
       3.9545161634549318e-03,
       4.8899112139089628e-03,
       5.9133597846795141e-03,
       7.0668491984957685e-03,
       8.4314424400194384e-03,
       1.0111020521181869e-02,
       1.2193069865332386e-02,
       1.4692894495666737e-02,
       1.7499888137509644e-02,
       2.0353055031603211e-02,
       2.2869855093105729e-02,
       2.4633725130110831e-02,
       2.5317115164213068e-02,
       2.4793837213310951e-02,
       2.3192772531288913e-02,
       2.0868780277482388e-02,
       1.8303511496475900e-02,
       1.5977201757051088e-02,
       1.4257195930872905e-02,
       1.3331644622708305e-02,
       1.3192932812138475e-02,
       1.3660594425678267e-02,
       1.4431816164603270e-02,
       1.5151385256259876e-02,
       1.5492278430193351e-02,
       1.5231225266263989e-02,
       1.4297925443659779e-02,
       1.2780983788951015e-02,
       1.0889283804535254e-02,
       8.8861144078574011e-03,
       7.0233923526266361e-03,
       5.4988368908542662e-03,
       4.4434391768754886e-03,
       3.9297916162201411e-03,
       3.9832623543391305e-03,
       4.5816420042560837e-03,
       5.6419071831224138e-03,
       7.0071435486950030e-03,
       8.4531220907505273e-03,
       9.7269400243558211e-03,
       1.0611971687423156e-02,
       1.0994671163099179e-02,
       1.0902061669702406e-02,
       1.0489646076980672e-02,
       9.9821047964225289e-03,
       9.5903851746582964e-03,
       9.4377295304821206e-03,
       9.5215457809027897e-03,
       9.7228581543024953e-03,
       9.8573486302087846e-03,
       9.7472439130392233e-03,
       9.2859442661513743e-03,
       8.4708227519770026e-03,
       7.3939955443089889e-03,
       6.1999404901693931e-03,
       5.0325661627771412e-03,
       3.9953700625636581e-03,
       3.1374375752149155e-03,
       2.4631659427485192e-03,
       1.9536325377179468e-03,
       1.5863710138811391e-03,
       1.3457544477740399e-03,
       1.2231373935468047e-03,
       1.2104996217773907e-03,
       1.2926302382092784e-03,
       1.4420095041531046e-03,
       1.6188292358654966e-03,
       1.7766316270724227e-03,
       1.8719822636020650e-03,
       1.8748519980087389e-03,
       1.7757892555370336e-03,
       1.5870889798633835e-03,
       1.3376348109504447e-03,
       1.0636595744390671e-03,
       7.9902381402190408e-04,
       5.6822246901199369e-04,
       3.8366975535995037e-04,
       2.4690734482912038e-04,
       1.5212987305029596e-04,
       9.0171398219280604e-05,
       5.1629948778467158e-05,
       2.8628346186686757e-05,
       1.5372481389457191e-05,
       7.9716621968390962e-06,
       3.9724906391045570e-06,
       1.8907325837018450e-06,
       8.5416133369256520e-07,
       3.6420537297391663e-07,
       1.4588838652304635e-07,
       5.4698430242977022e-08,
       1.9143156223214137e-08,
]
mbtr_now = vcat(desc[:]...)
tst = [isapprox(mbtr_now[x], mbtr_ref[x]; atol=1e-8) for x=1:length(mbtr_ref)]
println(@test all(tst))
end


text:
module Hashids
	export Hashid
    export encode, decode

	# Base code
	include("base.jl")
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6e3ef066-d115-11eb-2338-013a707dfe8a
begin
	using PlutoUI
	using DataFrames
	using Statistics
	using RollingFunctions
	using HTTP
	using JSON
	using Chain
	using Tables
	using PerceptualColourMaps
	using Gadfly
	using Match
	Gadfly.set_default_plot_size(16cm, 14cm)
end

# ╔═╡ 0be25163-f9c1-41bc-82cf-9f63ff8e609a
md"# Stabilizing the GPT architecture"

# ╔═╡ 654ad80a-a802-49d5-8373-b0b4056fd8f1
md"""
For this report we focused on stabilizing GPT's learning. We discovered that GPT was using an unusually low discount value for that large but unstable run that we showcased last time and found that changing this value improved stability but did not entirely fix the issueas the following graph indicates:
"""

# ╔═╡ ab745b36-164c-4534-be1f-a703a2010f3e
md"""
Note that all but one GPT seed diverges, though one does so only after nearly 1.5e8 steps. We are currently searching other hyperparameters with the updated discount value to see if these improve stability. Unfortunately all but `num_embeddings=1` and `num_embeddings=2` crashed due to memory issues on the smaller machines. Experiments with larger embedding numbers are currently running on rldl18 but have not yet converged:
"""

# ╔═╡ 9b848612-298c-4893-b8d5-15e7bf84ff24
md"""
## Status on other tasks
- Moving to Pytorch repository (almost done)
- Analyzing perception outputs (Logan)
"""

# ╔═╡ 8798b46c-1b0b-4d34-93ef-1d5062f8a632
md"""
## Architectures
This report compares the performance of a standard A2C baseline against a "GPT-in-the-middle" variant. The architectures are as follows:

#### A2C

- Convolution with
  - output-size $32$
  - kernel-shape $8\times 8$
  - stride $4\times 4$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $4\times 4$
  - stride $2\times 2$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $3\times 3$
  - stride $1\times 1$
- ReLU
- Linear layer with output-size $512$
- Parallel value head and action head

#### GPT
(differences in **bold**)

- Convolution with
  - output-size $32$
  - kernel-shape $8\times 8$
  - stride $4\times 4$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $4\times 4$
  - stride $2\times 2$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $3\times 3$
  - stride $1\times 1$
- ReLU
- **Linear layer with output-size depending on GPT size**
- **GPT torso**
- **one `action_hidden_size` layer**
- one layer projecting to action and value
"""

# ╔═╡ 623ed8e1-7f9e-4092-bf39-b9a3d90a7d23
colors(n) = cmap("I2"; N=n);

# ╔═╡ 5d3bfc57-6a5a-41e5-8f63-6902e4958936
begin
	struct HTMLDocument
		embedded
	end
	function Base.show(io::IO, mime::MIME"text/html", doc::HTMLDocument)
		println(io, "<html>")
		show(io, mime, doc.embedded)
		println(io, "</html>")
	end
end;

# ╔═╡ 41bd06a2-e7bd-46c6-9249-6e69223b0e11
begin
	HASRUA_ENDPOINT_URL = "http://rldl12.eecs.umich.edu:1200/v1/graphql"
	HASURA_ADMIN_SECRET = "hNuxFSTJMk28GMiFDpZDbviKtelFQamcb20UJiYZSIZ4LYsyLCCwQFYA6Y8HRSUo"
end;

# ╔═╡ 94185c13-b6d2-4337-b5ce-336f5e128032
function gql_query(query:: String; variables:: Dict = nothing)
	r = HTTP.request(
		"POST",
		HASRUA_ENDPOINT_URL;
		verbose=3,
		body= JSON.json(Dict("query" => query, "variables" => variables)) ,
		headers=Dict("x-hasura-admin-secret" => HASURA_ADMIN_SECRET),
	)
	@chain r begin
		_.body
		String
		JSON.parse(_, null=missing)
		_["data"]
	end
end;

# ╔═╡ 0f54629c-e045-47bf-90a5-36e44c05b8f0
function sweep_runs(sweep_ids::AbstractVector{Int}, max_step::Int)
	query = """
		query getSweepRuns(\$ids: [Int!], \$max_step: Int!) {
		  logs_less_than_step(args: {max_step: \$max_step}, where: {run: {sweep_id: {_in: \$ids}}}) {
			log
			run_id
			run {
			  metadata
			  sweep_id
			}
		  }
		}
  	"""
	rows = @chain gql_query(query; variables=Dict("ids" => sweep_ids, "max_step" => max_step)) begin
		_["logs_less_than_step"]
		map(d -> Dict(
				"run_id" => d["run_id"],
				"sweep_id" => d["run"]["sweep_id"],
				d["log"]...,
				d["run"]["metadata"]["parameters"]...
				), _)
		map(d -> Dict(
				d...,
				[name => get(d, name, nothing) for name in [
							"action_hidden_size",
							"gpt",
							"time",
							"gae",
							"gradient_clip",
							"nonlinearity",
							"normalize_observation",
							"normalize_torso_output",
							"optimizer",
							"num_embeddings",
							"save_interval",
							"save_path"
						]]...
				), _)
		collect
	end
	vcat(DataFrame.(rows)...)

end;

# ╔═╡ 67944dbb-0ebb-44d9-b6f3-79e8d5610f61
begin
	gpt_comparison_data = sweep_runs([98, 32], 300000000)
	gpt_comparison_data.gpt = replace(
		gpt_comparison_data.gpt,
		nothing => "Baseline"
	)
	plot(
        gpt_comparison_data,
		x=:step, y="episode return",
		group=:run_id,
		color=:gpt,
		Guide.xlabel("Step"),
		Guide.ylabel("Episode Return"),
		Geom.line,
		Scale.color_discrete(colors),
		Guide.colorkey(title="GPT size"),
		Guide.title("Breakout"),
		alpha=[0.5]
	) |> HTMLDocument
end

# ╔═╡ f4368e0e-dd79-4a80-bae0-84a09de65e97
begin
	plot(
        sweep_runs([99, 103], 300000000),
		x=:step, y="episode return",
		group=:run_id,
		color="num_embeddings",
		Guide.xlabel("Step"),
		Guide.ylabel("Episode Return"),
		Geom.line,
		Scale.color_discrete(colors),
		Guide.colorkey(title="#embeddings"),
		Guide.title("Breakout"),
		alpha=[0.5]
	) |> HTMLDocument
end

# ╔═╡ bbed8889-d401-4167-9370-9927aabbc83b
# @bind window_size Select(string.(collapse_runs ? (100:100:500) : (1:10)), default=collapse_runs ? "200" : "5")
window_size = "11";

# ╔═╡ Cell order:
# ╟─0be25163-f9c1-41bc-82cf-9f63ff8e609a
# ╟─654ad80a-a802-49d5-8373-b0b4056fd8f1
# ╟─67944dbb-0ebb-44d9-b6f3-79e8d5610f61
# ╟─ab745b36-164c-4534-be1f-a703a2010f3e
# ╠═f4368e0e-dd79-4a80-bae0-84a09de65e97
# ╟─9b848612-298c-4893-b8d5-15e7bf84ff24
# ╟─8798b46c-1b0b-4d34-93ef-1d5062f8a632
# ╟─623ed8e1-7f9e-4092-bf39-b9a3d90a7d23
# ╟─5d3bfc57-6a5a-41e5-8f63-6902e4958936
# ╟─6e3ef066-d115-11eb-2338-013a707dfe8a
# ╟─41bd06a2-e7bd-46c6-9249-6e69223b0e11
# ╟─94185c13-b6d2-4337-b5ce-336f5e128032
# ╟─0f54629c-e045-47bf-90a5-36e44c05b8f0
# ╟─bbed8889-d401-4167-9370-9927aabbc83b


text:
__precompile__()

module RStan

using RCall, DataFrames
export rstan

include("init.jl")
include("stan_interface.jl")

end # module


text:
@testset "basic MortalityTable" begin
    @testset "basic structure" begin
        q1 = UltimateMortality([i for i = 0:19], start_age=0)
        @test q1[0] == 0
        @test q1[1] == 1
        @test q1[0:1] == [0, 1]
        @test omega(q1) == 19
        @test_throws BoundsError q1[omega(q1) + 1]

        # non-zero start age
        q2 = UltimateMortality([i for i = 0:9], start_age=5)
        @test_throws BoundsError q2[4]
        @test q2[5] == 0

        # select strucutre
        s = [ia + d for ia = 0:5, d = 0:4]

        q3 = SelectMortality(s, q1, start_age=0)
        @test q3[0][0] == 0
        @test q3[0][0:1] == [0, 1]
        @test omega(q3[0]) == 19
        @test_throws BoundsError q3[omega(q3[0]) + 1]
        @test q3[0][19] == 19

        @test q3[5][5] == 5

        mt1 = MortalityTable(q3, q1)

        @test mt1.select[0][1] == 1
        @test mt1.ultimate[1] == 1

        mt2 = MortalityTable(q1)

        @test mt2.ultimate[0] == 0
        @test mt2[0] == 0

    end

    @testset "off-aligned select and ult" begin
        select_matrix = [(i + j - 1) / 100 for i in 0:10,j in 1:20]
        ult = UltimateMortality([i / 100 for i in 18:100], start_age=18)

        select = SelectMortality(select_matrix, ult)

        for issue_age in 0:10
            for dur in 1:30
                @test select[issue_age][issue_age + dur - 1] == (issue_age + dur - 1) / 100
            end
        end
    end

    # test time zero accumlated force
    @testset "accumulated force" begin
        q4 = UltimateMortality([0.1,0.3,0.6,1])

        @test survival(q4, 0) ≈ 1
        @test decrement(q4, 0) ≈ 0

        @test survival(q4, 1) ≈ 0.9
        @test decrement(q4, 1) ≈ 0.1

        @test survival(q4, 1, 1) ≈ 1.0
        @test survival(q4, 1, 2) ≈ 0.7
        @test decrement(q4, 1, 1) ≈ 0.0
        @test decrement(q4, 1, 2) ≈ 0.3

        @test survival(q4, 1, 4) ≈ 0.0
        @test decrement(q4, 1, 4) ≈ 1.0

        @test survival(q4, -1) ≈ 1.0
        @test survival(q4, 4, -1) ≈ 1.0
        @test decrement(q4, -1) ≈ 0.0
        @test decrement(q4, 4, -1) ≈ 0.0
    end

    @testset "Metadata" begin
        d = TableMetaData()

        @test isnothing(d.name)

        d = TableMetaData(name="test")
        @test d.name == "test"
    end

    @testset "mortality_vector" begin
        v = [i for i in 3:10]
        q = mortality_vector(v, start_age=3)
        @test q[3] == 3
        @test q[10] == 10

        q = mortality_vector(collect(0:5))
        @test q[0] == 0
        @test q[5] == 5
    end

end


text:
import NLPModels: grad, cons, jac

"""
unconstrained: return the infinite norm of the gradient of the objective function

required: state.gx (filled if nothing)
"""
function unconstrained_check(pb    :: AbstractNLPModel,
                             state :: NLPAtX;
                             pnorm :: Float64 = Inf,
                             kwargs...)

    if state.gx == nothing # should be filled if empty
        update!(state, gx = grad(pb, state.x))
    end

    res = norm(state.gx, pnorm)

    return res
end

"""
unconstrained 2nd: check the norm of the gradient and the smallest
                   eigenvalue of the hessian.

required: state.gx, state.Hx (filled if nothing)
"""
function unconstrained2nd_check(pb    :: AbstractNLPModel,
                                state :: NLPAtX;
                                pnorm :: Float64 = Inf,
                                kwargs...)

    if state.gx == nothing # should be filled if empty
        update!(state, gx = grad(pb, state.x))
    end
    if state.Hx == nothing
        update!(state, Hx = hess(pb, state.x))
    end

    res = max(norm(state.gx, pnorm),
              max(- eigmin(state.Hx + state.Hx' - diagm(0 => diag(state.Hx))), 0.0))

    return res
end

"""
optim_check_bounded: gradient of the objective function projected

required: state.gx (filled if void)
"""
function optim_check_bounded(pb    :: AbstractNLPModel,
                             state :: NLPAtX;
                             pnorm :: Float64 = Inf,
                             kwargs...)

    if state.gx == nothing # should be filled if void
     update!(state, gx = grad(pb, state.x))
    end

    proj = max.(min.(state.x - state.gx, pb.meta.uvar), pb.meta.lvar)
    gradproj = state.x - proj
    res = norm(gradproj, pnorm)

    return res
end

"""
constrained: return the violation of the KKT conditions
length(lambda) > 0
"""
function _grad_lagrangian(pb    :: AbstractNLPModel,
                          state :: NLPAtX)

 if (pb.meta.ncon == 0) & !has_bounds(pb)
  return state.gx
 elseif pb.meta.ncon == 0
  return state.gx + state.mu
 else
  return state.gx + state.mu + state.Jx' * state.lambda
 end
end

function _sign_multipliers_bounds(pb    :: AbstractNLPModel,
                                  state :: NLPAtX)
 if has_bounds(pb)
  return vcat(min.(max.( state.mu,0.0), - state.x + pb.meta.uvar),
              min.(max.(-state.mu,0.0),   state.x - pb.meta.lvar))
 else
  return zeros(0)
 end
end

function _sign_multipliers_nonlin(pb    :: AbstractNLPModel,
                                  state :: NLPAtX)
 if pb.meta.ncon == 0
  return zeros(0)
 else
  return vcat(min.(max.( state.lambda,0.0), - state.cx + pb.meta.ucon),
              min.(max.(-state.lambda,0.0),   state.cx - pb.meta.lcon))
 end
end

function _feasibility(pb    :: AbstractNLPModel,
                      state :: NLPAtX)
 if pb.meta.ncon == 0
  return vcat(max.(  state.x  - pb.meta.uvar,0.0),
              max.(- state.x  + pb.meta.lvar,0.0))
 else
  return vcat(max.(  state.cx - pb.meta.ucon,0.0),
              max.(- state.cx + pb.meta.lcon,0.0),
              max.(  state.x  - pb.meta.uvar,0.0),
              max.(- state.x  + pb.meta.lvar,0.0))
 end
end

"""
KKT: verifies the KKT conditions

required: state.gx
+ if bounds: state.mu
+ if constraints: state.cx, state.Jx, state.lambda
"""
function KKT(pb    :: AbstractNLPModel,
             state :: NLPAtX;
             pnorm :: Float64 = Inf,
             kwargs...)

    #Check the gradient of the Lagrangian
    gLagx      = _grad_lagrangian(pb, state)
    #Check the complementarity condition for the bounds
    res_bounds = _sign_multipliers_bounds(pb, state)
    #Check the complementarity condition for the constraints
    res_nonlin = _sign_multipliers_nonlin(pb, state)
    #Check the feasibility
    feas       = _feasibility(pb, state)

    res = vcat(gLagx, feas, res_bounds, res_nonlin)

    return norm(res, pnorm)
end