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

using Markdown
using InteractiveUtils

# ╔═╡ fd82b765-7a62-4f85-bf60-28c83b1731aa
begin
	import Pkg
	Pkg.develop("PoreMatMod")
end

# ╔═╡ 3bb3966e-a7dd-46c1-b649-33169ce424d2
using PoreMatMod, PlutoUI

# ╔═╡ d3eef6f4-ff15-47f3-8686-2c0cb0fb882d
using PoreMatMod.ExampleHelpers

# ╔═╡ 8e4fd703-f53e-4056-9466-66f07bacad8d
md"## Example: substructure replacement in a high-symmetry parent"

# ╔═╡ 49ce8a40-3c99-48a8-9403-f5538d1d9a67
input_file_message()

# ╔═╡ 98ee7a1e-1c27-4819-b752-01602cde2799
xtal_folder_message()

# ╔═╡ f5a08f77-9830-4e8e-b746-6a1d61f39009
rc[:paths][:crystals]

# ╔═╡ 8c4c9c65-61a3-4ac8-9413-48aab281237e
moiety_folder_message()

# ╔═╡ 4c5e0230-4db5-49de-b9d7-152ebc834d48
rc[:paths][:moieties]

# ╔═╡ e7be5fba-b120-47f5-a952-8577ee847d99
md"""
!!! example \"the task\"
	We have a crystallographic unit cell of NiPyC. We wish to append methyl groups onto the PyC linkers.

**Parent crystal structure**: first, we read in the `.cif` file describing the parent structure, which is in a high-symmetry representation, meaning this structure must be reproduced according to the symmetry rules of the unit cell's space group. Normally, structures are automatically transformed to P1 symmetry, so we must specify that we wish to preserve the original unit cell.
"""

# ╔═╡ c0e2ffbd-fb85-4187-b8b5-73edea90969a
begin
	# read in the parent xtal, keeping it in its original space group
	parent = Crystal("NiPyC_experiment.cif", convert_to_p1=false) # load cif file
	infer_bonds!(parent, true) 									  # infer bonds
	view_structure(parent)										  # view structure
end

# ╔═╡ 0bbe1bf4-89e8-4372-8f6a-93fc89334b00
fragment_construction_note()

# ╔═╡ f02fdbf1-2c43-40df-a678-e3c535751f3e
md"""
**Query fragment**: next, we construct a query fragment to match the linker in NiPyC.
"""

# ╔═╡ 16faf8c8-d5c2-4755-945b-546cdac27350
query = moiety("PyC.xyz");

# ╔═╡ 75610c4a-9218-406c-b28f-8e299e197135
view_query_or_replacement("PyC.xyz")

# ╔═╡ 5e90af51-f777-47f9-980b-19bca38de5d4
with_terminal() do
	display_query_or_replacement_file("PyC.xyz")
end

# ╔═╡ 70468997-24fa-4c1b-9f02-379823f97db8
md"""
**Replacement fragment**: then, we construct a replacement fragment as a corrected version of the query fragment (with hydrogen atoms appropriately appended).
"""

# ╔═╡ aac504f5-080e-47b8-9a33-d43c16dc87e7
replacement = moiety("PyC-CH3.xyz");

# ╔═╡ 976795e6-e696-4a37-9cd1-83b4d2695a96
view_query_or_replacement("PyC-CH3.xyz")

# ╔═╡ 7f82be68-558d-4dce-8492-bcc585ddf448
with_terminal() do
	display_query_or_replacement_file("PyC-CH3.xyz")
end

# ╔═╡ 2e3299b5-3b63-48d9-90d7-b8dde5715907
md"""
**Find and replace**: finally, we search the parent MOF for the query fragment and effect the replacement.
"""

# ╔═╡ 0fd8b373-9fec-4b18-b1a6-1b11741e5215
search = query in parent

# ╔═╡ 0b287f09-a121-430a-a749-dc93492d1680
child = substructure_replace(search, replacement, nb_loc=1, wrap=false)

# ╔═╡ d8ea8d5e-f17a-412b-8461-15ba6d9621ec
write_cif(child)

# ╔═╡ eae2225c-40f0-4d68-a9a2-43a39a82f029
view_structure(child)

# ╔═╡ 000144c6-be54-4774-b449-698e9da0741a
write_cif(child, "data/crystals/symmetry_child.cif")

# ╔═╡ 47efaf2e-6758-42c5-aab2-80c1d7725b4a
md"""
**making the super-cell** once we have the high-symmetry structure modified, we can apply the symmetry rules and get a replicated super-cell.
"""

# ╔═╡ 36a1be30-db0e-4c45-8894-859e36793482
begin
	p1_child = Crystal("symmetry_child.cif", check_overlap=false)
	p1_child = replicate(p1_child, (2, 2, 2))
	infer_bonds!(p1_child, true)
	write_cif(p1_child, "supercell.cif")
	view_structure(p1_child)
end

# ╔═╡ Cell order:
# ╟─fd82b765-7a62-4f85-bf60-28c83b1731aa
# ╟─8e4fd703-f53e-4056-9466-66f07bacad8d
# ╠═3bb3966e-a7dd-46c1-b649-33169ce424d2
# ╠═d3eef6f4-ff15-47f3-8686-2c0cb0fb882d
# ╟─49ce8a40-3c99-48a8-9403-f5538d1d9a67
# ╠═98ee7a1e-1c27-4819-b752-01602cde2799
# ╠═f5a08f77-9830-4e8e-b746-6a1d61f39009
# ╟─8c4c9c65-61a3-4ac8-9413-48aab281237e
# ╠═4c5e0230-4db5-49de-b9d7-152ebc834d48
# ╟─e7be5fba-b120-47f5-a952-8577ee847d99
# ╠═c0e2ffbd-fb85-4187-b8b5-73edea90969a
# ╟─0bbe1bf4-89e8-4372-8f6a-93fc89334b00
# ╟─f02fdbf1-2c43-40df-a678-e3c535751f3e
# ╠═16faf8c8-d5c2-4755-945b-546cdac27350
# ╟─75610c4a-9218-406c-b28f-8e299e197135
# ╟─5e90af51-f777-47f9-980b-19bca38de5d4
# ╟─70468997-24fa-4c1b-9f02-379823f97db8
# ╠═aac504f5-080e-47b8-9a33-d43c16dc87e7
# ╟─976795e6-e696-4a37-9cd1-83b4d2695a96
# ╟─7f82be68-558d-4dce-8492-bcc585ddf448
# ╟─2e3299b5-3b63-48d9-90d7-b8dde5715907
# ╠═0fd8b373-9fec-4b18-b1a6-1b11741e5215
# ╠═0b287f09-a121-430a-a749-dc93492d1680
# ╠═d8ea8d5e-f17a-412b-8461-15ba6d9621ec
# ╠═eae2225c-40f0-4d68-a9a2-43a39a82f029
# ╠═000144c6-be54-4774-b449-698e9da0741a
# ╟─47efaf2e-6758-42c5-aab2-80c1d7725b4a
# ╠═36a1be30-db0e-4c45-8894-859e36793482

ext: jl


text:
# Example GUI to look at a dataset

using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie
using GridLayoutBase
using Observables

using AlgebraOfGraphics
using AlgebraOfGraphics: linear, smooth, density, AlgebraicList, layoutplot!

function clean!(g::GridLayout)
    foreach(delete!, Any, g)
    foreach(GridLayoutBase.remove_from_gridlayout!, g.content)
end

function gui!(scene, layout, df)
    table = data(df)

    plot_options = [("none", nothing), ("scatter", Scatter), ("lines", Lines), ("barplot", BarPlot)]
    plot_menu = LMenu(scene, options = plot_options, textsize=30, width=500)

    layout[1, 1] = LText(scene, "plot type", textsize=30)
    layout[1, 2] = plot_menu

    analysis_options = [("none", nothing), ("linear", linear), ("smooth", smooth),
        ("density", density)]

    analysis_menu = LMenu(scene, options = analysis_options, textsize=30, width=500)

    layout[2, 1] = LText(scene, "analysis", textsize=30)
    layout[2, 2] = analysis_menu

    axis_options = vcat([("None", nothing)], [(s, Symbol(s)) for s in propertynames(table.data)])

    styles = ["x", "y", "color", "marker", "markersize", "linestyle", "layout_x", "layout_y"]
    style_menus = [LMenu(scene, options=axis_options, textsize=30, width=500) for style in styles]
    N = length(styles)

    for i in 1:N
        layout[2+i, 1] = LText(scene, styles[i], textsize=30)
        layout[2+i, 2] = style_menus[i]
    end

    plot_button = LButton(scene, label="Plot", textsize=30)
    add_and_plot_button = LButton(scene, label="Add and Plot", buttoncolor="powderblue", textsize=30)

    layout[1, 3] = plot_button
    layout[1, 4] = add_and_plot_button

    axs = layout[1:2+N, 5] = GridLayout(1, 2) # for plot and legend
    colsize!(axs, 2, Relative(1/5)) # make legend a fourth of the plot's width

    state = AlgebraicList()

    function add_and_plot()
        acc = table
        T = plot_menu.selection[]
        isnothing(T) || (acc *= spec(T))
        an = analysis_menu.selection[]
        isnothing(an) || (acc *= an)
        for (st, st_menu) in zip(styles, style_menus)
            sym = st === "x" ? Symbol("1") : st === "y" ? Symbol("2") : Symbol(st)
            val = st_menu.selection[]
            isnothing(val) || (acc *= style(; sym => val))
        end
        clean!(axs)
        state += acc
        layoutplot!(scene, axs, state)
    end

    on(plot_button.clicks) do _
        state = AlgebraicList()
        add_and_plot()
    end
    on(_ -> add_and_plot(), add_and_plot_button.clicks)

    return scene
end

function gui(df; kwargs...)
    scene, layout = MakieLayout.layoutscene(; kwargs...)
    return gui!(scene, layout, df)
end

# Try it on iris dataset

using RDatasets: dataset

iris = dataset("datasets", "iris")

display(gui(iris))


text:
using Persa
using Surprise
using DatasetsCF

ds = DatasetsCF.MovieLens()

for (ds_train, ds_test) in Persa.KFolds(ds, 10)
  model = Surprise.IRSVD(ds_train)

  Persa.train!(model, ds_train)

  print(Persa.aval(model, ds_test, Persa.recommendation(ds)))
end


text:
using PolyaGammaSamplers
using Documenter

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

makedocs(;
    modules=[PolyaGammaSamplers],
    authors="Iván Gutiérrez <ivangutierrez1988@gmail.com> and contributors",
    repo="https://github.com/igutierrezm/PolyaGammaSamplers.jl/blob/{commit}{path}#{line}",
    sitename="PolyaGammaSamplers.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://igutierrezm.github.io/PolyaGammaSamplers.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Getting Started" => "getting_started.md",
        "Samplers" => "samplers.md",
    ],
)

deploydocs(;
    repo="github.com/igutierrezm/PolyaGammaSamplers.jl",
)


text:
using LowRankModels
using Distributions

# This is just to create a low-rank representation of some count data, I don't care how accurate it is.
m,n = 100, 50;
k = 2;
A = rand(Poisson(),(m,n));
losses = convert(Array{Loss,1}, fill(PoissonLoss(10),n));
rx, ry = ZeroReg(), ZeroReg();
g_pre = GLRM(A, losses, rx, ry, k, scale=false, offset=false);
X_real, Y_real, ch = fit!(g_pre, params=Params(0.00001,min_stepsize=0.00000000001,max_iter=10000));

# Now we do the actual model using the perfect data and try to recapture it.
A_real = impute(losses, X_real'*Y_real);

g = GLRM(A_real, losses, rx, ry, k, scale=true, offset=true);
X, Y, ch = fit!(g, params=Params(0.00001,min_stepsize=0.00000000001,max_iter=10000));
U = X'*Y;
A_imputed = impute(losses, X'*Y);

error_metric(g)
errors(Domain[l.domain for l in losses], losses, U, A_real);


text:
let
    @test sparse((jactestmeta([1,2], [3,4], 3) do x,y,t
                      t>1 ? x[1] = y[2] : x[2] = y[1]
                  end)) == sparse([2, 1], [1, 2], true)

    function h(y, x, t)
        if t > 0
            y[2, :] .= x
        else
            y[1, :] .= x
        end
    end

    @test sparse(jactestmeta(h, rand(3,3), ones(3), 3)) == sparse([1,4,7,2,5,8], [1,2,3,1,2,3], true, 9, 3)
end


text:
wrap(base) = joinpath("wrapper", base * ".jl")

# used by Clang.jl
include("CEnum.jl")
using .CEnum

# all type definitions
include(wrap("manual_commons"))
include(wrap("commons"))

# wrappers for scip headers
include(wrap("scip_bandit"))
include(wrap("scip_benders"))
include(wrap("scip_branch"))
include(wrap("scip_compr"))
include(wrap("scip_concurrent"))
include(wrap("scip_conflict"))
include(wrap("scip_cons"))
include(wrap("scip_copy"))
include(wrap("scip_cut"))
include(wrap("scip_datastructures"))
include(wrap("scip_debug"))
include(wrap("scip_dialog"))
include(wrap("scip_disp"))
include(wrap("scip_event"))
include(wrap("scip_expr"))
include(wrap("scip_general"))
include(wrap("scip_heur"))
include(wrap("scip_lp"))
include(wrap("scip_mem"))
include(wrap("scip_message"))
include(wrap("scip_nlp"))
include(wrap("scip_nodesel"))
include(wrap("scip_nonlinear"))
include(wrap("scip_numerics"))
include(wrap("scip_param"))
include(wrap("scip_presol"))
include(wrap("scip_pricer"))
include(wrap("scip_probing"))
include(wrap("scip_prob"))
include(wrap("scip_prop"))
include(wrap("scip_randnumgen"))
include(wrap("scip_reader"))
include(wrap("scip_relax"))
include(wrap("scip_reopt"))
include(wrap("scip_sepa"))
include(wrap("scip_sol"))
include(wrap("scip_solve"))
include(wrap("scip_solvingstats"))
include(wrap("scip_table"))
include(wrap("scip_timing"))
include(wrap("scip_tree"))
include(wrap("scip_validation"))
include(wrap("scip_var"))

# default SCIP plugins
include(wrap("scipdefplugins"))

# other public headers
include(wrap("pub_bandit_epsgreedy"))
include(wrap("pub_bandit_exp3"))
include(wrap("pub_bandit"))
include(wrap("pub_bandit_ucb"))
include(wrap("pub_benderscut"))
include(wrap("pub_benders"))
include(wrap("pub_branch"))
include(wrap("pub_compr"))
include(wrap("pub_conflict"))
include(wrap("pub_cons"))
include(wrap("pub_cutpool"))
include(wrap("pub_dialog"))
include(wrap("pub_disp"))
include(wrap("pub_event"))
include(wrap("pub_fileio"))
include(wrap("pub_heur"))
include(wrap("pub_history"))
include(wrap("pub_implics"))
include(wrap("pub_lp"))
include(wrap("pub_matrix"))
include(wrap("pub_message"))
include(wrap("pub_misc"))
include(wrap("pub_misc_linear"))
include(wrap("pub_misc_select"))
include(wrap("pub_misc_sort"))
include(wrap("pub_nlp"))
include(wrap("pub_nodesel"))
include(wrap("pub_paramset"))
include(wrap("pub_presol"))
include(wrap("pub_pricer"))
include(wrap("pub_prop"))
include(wrap("pub_reader"))
include(wrap("pub_relax"))
include(wrap("pub_reopt"))
include(wrap("pub_sepa"))
include(wrap("pub_sol"))
include(wrap("pub_table"))
include(wrap("pub_tree"))
include(wrap("pub_var"))

# all constraint types
include(wrap("cons_abspower"))
include(wrap("cons_and"))
include(wrap("cons_benders"))
include(wrap("cons_benderslp"))
include(wrap("cons_bivariate"))
include(wrap("cons_bounddisjunction"))
include(wrap("cons_cardinality"))
include(wrap("cons_components"))
include(wrap("cons_conjunction"))
include(wrap("cons_countsols"))
include(wrap("cons_cumulative"))
include(wrap("cons_disjunction"))
include(wrap("cons_indicator"))
include(wrap("cons_integral"))
include(wrap("cons_knapsack"))
include(wrap("cons_linear"))
include(wrap("cons_linking"))
include(wrap("cons_logicor"))
include(wrap("cons_nonlinear"))
include(wrap("cons_orbisack"))
include(wrap("cons_orbitope"))
include(wrap("cons_or"))
include(wrap("cons_pseudoboolean"))
include(wrap("cons_quadratic"))
include(wrap("cons_setppc"))
include(wrap("cons_soc"))
include(wrap("cons_sos1"))
include(wrap("cons_sos2"))
include(wrap("cons_superindicator"))
include(wrap("cons_symresack"))
include(wrap("cons_varbound"))
include(wrap("cons_xor"))

# nonlinear expressions
include(wrap("pub_expr"))
include(wrap("expr_manual"))

# SCIP_CALL: macro to check return codes, inspired by @assert
macro SC(ex)
    return :(@assert $(esc(ex)) == SCIP_OKAY)
end

const SCIP_STATISTICS_FUNCS = [
    :SCIPprintStatusStatistics,
    :SCIPprintTimingStatistics,
    :SCIPprintOrigProblemStatistics,
    :SCIPprintTransProblemStatistics,
    :SCIPprintPresolverStatistics,
    :SCIPprintConstraintStatistics,
    :SCIPprintConstraintTimingStatistics,
    :SCIPprintPropagatorStatistics,
    :SCIPprintConflictStatistics,
    :SCIPprintSeparatorStatistics,
    :SCIPprintPricerStatistics,
    :SCIPprintBranchruleStatistics,
    :SCIPprintHeuristicStatistics,
    :SCIPprintCompressionStatistics,
    :SCIPprintLPStatistics,
    :SCIPprintNLPStatistics,
    :SCIPprintRelaxatorStatistics,
    :SCIPprintTreeStatistics,
    :SCIPprintRootStatistics,
    :SCIPprintSolutionStatistics,
    :SCIPprintConcsolverStatistics,
    :SCIPprintBendersStatistics,
    :SCIPprintStatistics,
    :SCIPprintReoptStatistics,
    :SCIPprintBranchingStatistics
]


text:
#Test RefractiveIndex abstraction
using CavChaos

idx0 = UniformIdx();
idx1 = GaussianIdx();
idx2 = TaylorIdx();
idx3 = ChebFourIdx();

#RefractiveIndex test function ================================================
function testidxfunc(idx::RefractiveIndex,rval::Array{Float64,1},thetaval::Array{Float64,1})
    println("========== Testing $(summary(idx)) ==========");

    #Setup test quantities
    r = Array(Float64,length(rval),length(thetaval));
    theta = Array(Float64,length(rval),length(thetaval));
    for j=1:length(thetaval)
        for i=1:length(rval)
            r[i,j] = rval[i];
            theta[i,j] = thetaval[j];
        end
    end

    #Scalar functions
    tic();
    n0 = similar(theta);
    n1 = similar(theta); dr_n1 = similar(theta); dtheta_n1 = similar(theta);
    for i=1:length(theta)
        n0[i] = nfunc(idx,r[i],theta[i]);
        n1[i],dr_n1[i],dtheta_n1[i] = nderiv(idx,r[i],theta[i]);
    end
    println("Scalar functions compute in $(toq())s.");

    #Array functions
    tic();
    n0A = nfunc(idx,r,theta);
    n1A,dr_n1A,dtheta_n1A = nderiv(idx,r,theta);
    println("Array functions compute in $(toq())s.");

    #Test agreement
    print("Testing function agreement... ");
    if !(n0 == n1 == n0A == n1A)
        println("n agreement failed!");
        println("Returning [r theta n0 n1 n0A n1A].");
        where = find(!(n0 .== n1 .== n0A .== n1A));
        return [r[where] theta[where] n0[where] n1[where] n0A[where] n1A[where]];
    end
    if !(dr_n1 == dr_n1A)
        println("dr_n agreement failed!");
        println("Returning [r theta dr_n1 dr_n1A].");
        where = find(!(dr_n1 .== dr_n1A));
        return [r[where] theta[where] dr_n1[where] dr_n1A[where]];
    end
    if !(dtheta_n1 == dtheta_n1A)
        println("dtheta_n agreement failed!");
        println("Returning [r theta dtheta_n1 dtheta_n1A].");
        where = find(!(dtheta_n1 .== dtheta_n1A));
        return [r[where] theta[where] dtheta_n1[where] dtheta_n1A[where]];
    end
    println("passed.");

    #Estimate derivatives
    dr = rval[1]; dtheta = thetaval[2];
    dr_n0 = (n0[2:end,:]-n0[1:end-1,:])/dr;
    dtheta_n0 = (n0[:,2:end]-n0[:,1:end-1])/dtheta;
    dr_n1m = 0.5*(dr_n1[1:end-1,:]+dr_n1[2:end,:]);
    dtheta_n1m = 0.5*(dtheta_n1[:,1:end-1]+dtheta_n1[:,2:end]);
    error_r = maximum(abs(dr_n0 - dr_n1m));
    error_theta = maximum(abs(dtheta_n0 - dtheta_n1m));

    #Test accuracy
    tol = 1e-3;
    print("Testing accuracy of derivatives... ");
    if error_r > tol
        println("dr_n error is $(string(error_r))!");
        println("Returning [r theta dr_n0 dr_n1m].");
        where = find(abs(dr_n0 - dr_n1m) .> tol);
        return [r[1:end-1,:][where] theta[1:end-1,:][where] dr_n0[where] dr_n1m[where]];
    end
    if error_theta > tol
        println("dtheta_n error is $(string(error_theta))!");
        println("Returning [r theta dtheta_n0 dtheta_n1m].");
        where = find(abs(dtheta_n0 - dtheta_n1m) .> tol);
        return [r[:,1:end-1][where] theta[:,1:end-1][where] dtheta_n0[where] dtheta_n1m[where]];
    end
    println("passed.");

    return "passed";
end

#Run tests ====================================================================
function testidxs()
    r = linspace(0,1.5,101)[2:end];
    theta = linspace(0,2*pi,200);

    println("############### REFRACTIVE INDEX TEST ###############");

    testinfo0 = testidxfunc(idx0,r,theta);
    testinfo1 = testidxfunc(idx1,r,theta);
    testinfo2 = testidxfunc(idx2,r,theta);
    testinfo3 = testidxfunc(idx3,r,theta);

    print("############### ");
    if testinfo0 == testinfo1 == testinfo2 == testinfo3 == "passed"
        print("PASSED ALL");
    else
        print("FAILED");
    end
    println(" ###############\n");

    return (testinfo0,testinfo1,testinfo2,testinfo3);
end

testidxs()


text:
function pairwise{T}(v::Vector{T})

end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
using NeXLCore
using Test

@testset "Staging" begin
    stg_coords = Dict(:X=>10.0, :Y=>5.0)
    img_coords = Dict(:X=>0.027777, :Y=>-0.016666 )
    rot = deg2rad(50.0)
    i2s = image2stage(NeXLCore.DefaultStageMapping, stg_coords, img_coords, rot)
    s2i = stage2image(NeXLCore.DefaultStageMapping, stg_coords, i2s, rot)
    @test isapprox(s2i[:X], img_coords[:X], atol=1.0e-6)
    @test isapprox(s2i[:Y], img_coords[:Y], atol=1.0e-6)
end


text:
using STC.SUniward
using Base.Test
using Images

function testreflectedsetindex()
  a=zeros(20,30)
  b= padarray(a,Pad(:symmetric,[16,16],[16,16]))
  @testset begin
    for j in 1:size(a,2)
        for i in 1:size(a,1)
            a[i,j]=1
            SUniward.reflectedsetindex!(b,i,j,1)
            c = padarray(a,Pad(:symmetric,[16,16],[16,16]))
            @test sum(abs.(c.parent - b.parent)) == 0
        end
    end
  end
end

function testsuniwardcosts()
  cover=load("suniward.jld","cover")
  rho=load("suniward.jld","rho")
  @testset begin
      @test sum(abs.(rho-SUniward.suniwardcosts(cover)))<1e-6
  end
end

function testoneitemconv2_1(flt)
    # a=rand(1:30,30,30)
    a=rand(30,30)
    # i,j,v=3,3,1
    (i,j,v)=rand(1:size(a,1)),rand(1:size(a,2)),rand()
    aa = padarray(SUniward.sameconv2(a,flt),Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    a[i,j]+=v
    sum(abs.(SUniward.sameconv2(a,flt) - SUniward.oneitemconv2!(aa,flt,i,j,v).parent[size(flt,1)+1:end-size(flt,1),size(flt,2)+1:end-size(flt,2)]))
end

function testoneitemconv2_2(flt)
    a=rand(30,30)
    aa = padarray(zeros(size(a)),Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    for i in 1:size(a,1)
      for j in 1:size(a,2)
        SUniward.oneitemconv2!(aa,flt,i,j,a[i,j])
      end
    end
    sum(abs.(SUniward.sameconv2(a,flt) - aa.parent[size(flt,1)+1:end-size(flt,1),size(flt,2)+1:end-size(flt,2)]))
end

function testweightedconv(flt)
    sigma=1
    cover=rand(30,30)
    stego=rand(30,30)
    diff=cover-stego

    paddedcover= padarray(cover,Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    paddeddiff= padarray(diff,Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    invr = padarray(zeros(size(diff)),Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))

    rc = SUniward.sameconv2(paddedcover.parent, flt)
    invr.parent.=1./(abs.(rc)+sigma)

    aa = padarray(zeros(size(diff)),Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    for i in 1:size(diff,1)
      for j in 1:size(diff,2)
        SUniward.oneitemconv2!(aa,invr,flt,i,j,diff[i,j])
      end
    end

    costs=SUniward.sameconv2(diff,flt).*invr.parent[size(flt,1)+1:end-size(flt,1),size(flt,2)+1:end-size(flt,2)]
    sum(abs.(costs - aa.parent[size(flt,1)+1:end-size(flt,1),size(flt,2)+1:end-size(flt,2)]))
end

function testincrementalcosts(flt)
    sigma=1
    cover=rand(30,30)
    stego=rand(30,30)
    diff=cover-stego

    paddedcover= padarray(cover,Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    paddeddiff= padarray(diff,Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    invr = padarray(zeros(size(diff)),Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))

    rc = SUniward.sameconv2(paddedcover.parent, flt)
    invr.parent.=1./(abs(rc)+sigma)

    aa = padarray(zeros(size(diff)),Pad(:symmetric,[size(flt,1),size(flt,2)],[size(flt,1),size(flt,2)]))
    distortion=zero(eltype(cover))
    (ul,vl)=(cld(size(flt,1),2),cld(size(flt,2),2))
    (up,vp)=(fld(size(flt,1),2),fld(size(flt,2),2))
    for i in 1:size(diff,1)
      for j in 1:size(diff,2)
        idxs=(max(1,i-up):min(size(diff,1),i+ul-1),max(1,j-vp):min(size(diff,2),j+vl-1))
        # println((i,j), " ",idxs)
        distortion-=sumabs(aa[idxs...])
        SUniward.oneitemconv2!(aa,invr,flt,i,j,diff[i,j])
        distortion+=sumabs(aa[idxs...])
      end
    end

    costs=SUniward.sameconv2(diff,flt).*invr.parent[size(flt,1)+1:end-size(flt,1),size(flt,2)+1:end-size(flt,2)]
    sum(abs.(sum(abs.(costs)-distortion)))
end

function testincrementalsuniward()
    cover=UInt8.(rand(0:255,30,30))
    stego=UInt8.(rand(0:255,30,30))
    image=SUniward.IncrementalSUniward(cover)
    for i in 1:size(cover,1)
      for j in 1:size(cover,2)
        image[i,j]=stego[i,j]
      end
    end
    return(abs(SUniward.suniwardcosts(Float64.(cover),Float64.(stego))-image.distortion))
end

function testincrementalsuniward2()
    cover=UInt8.(rand(0:255,30,30))
    stego=UInt8.(rand(0:255,30,30))
    image=SUniward.IncrementalSUniward(cover)
    for i in 1:size(cover,1)
      for j in 1:size(cover,2)
        image[sub2ind(size(cover),i,j)]=stego[i,j]
      end
    end
    return(abs(SUniward.suniwardcosts(Float64.(cover),Float64.(stego))-image.distortion))
end

function testLSBcostfun()
    cover=UInt8.(rand(0:255,30,30))
    pidx=rand(1:length(cover))
    cover[pidx]=8
    stego=deepcopy(cover)
    stego[pidx]=9

    image=SUniward.IncrementalSUniward(cover);
    w0,w1=SUniward.LSBcostfun(image,pidx)

    r=sum(abs.(w1-SUniward.suniwardcosts(Float64.(cover),Float64.(stego))))

    cover=UInt8.(rand(0:255,30,30))
    pidx=rand(1:length(cover))
    cover[pidx]=9
    stego=deepcopy(cover)
    stego[pidx]=8

    image=SUniward.IncrementalSUniward(cover);
    w0,w1=SUniward.LSBcostfun(image,pidx)

    r+=sum(abs.(w0-SUniward.suniwardcosts(Float64.(cover),Float64.(stego))))
    return(r)
end

function testtrypmone()
    cover=UInt8.(rand(0:255,30,30))
    image=SUniward.IncrementalSUniward(cover);

    d1=0.0
    d2=0.0
    for pidx in 1:length(cover)
        d1 += min(SUniward.tryvalue(image,cover[pidx]+1,pidx),SUniward.tryvalue(image,cover[pidx]-1,pidx))
        d2 += min(SUniward.trypmone(image,pidx)...)

        s = rand([-1,1]);
        s = (image[pidx] == 0)? 1 : s;
        s = (image[pidx] == 255)? -1 : s;
        image[pidx]=cover[pidx] + s
    end

    return(abs(d1-d2))
end

function testbinaryembedding(s=8,h=3)
    println("binary embedding with s = $s h = $h")
    payload = 0.4;
    cover=rawview(channelview(load("1.pgm")));
    l=size(cover);
    s=div(s,2);
    cover = cover[div(l[1],2)-s+1:div(l[1],2)+s,div(l[2],2)-s+1:div(l[2],2)+s];

    numofcols=Int(round(1/payload));
    hhat=rand(1:2^h-1,numofcols);
    numofblocks=Int(floor(length(cover)/numofcols));
    message=rand(0:1,numofblocks);
    embpath = randperm(length(cover));

    stegos = SUniward.IncrementalSUniward(cover,2^h;sigma=1.0,T=Float16);
    @time (stego,distortion)=STC.variablecoding(cover,stegos,message,hhat,h,embpath);
    δ = Int.(cover) - Int.(stego.image)
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego.image),sum(δ .> 0),sum(δ .< 0)));

    stegos = SUniward.IncrementalSUniward(cover,2^h;sigma=1.0,T=Float16);
    @time (stego,distortion)=STC.variablecoding(cover,stegos,message,hhat,h,embpath);
    δ = Int.(cover) - Int.(stego.image)
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego.image),sum(δ .> 0),sum(δ .< 0)));

    stegos = [SUniward.SUniwardAdd(cover) for i in 1:2^h];
    @time (stego,distortion)=STC.variablecoding(cover,stegos,message,hhat,h,embpath);
    δ = Int.(cover) - Int.(stego.image)
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego.image),sum(δ .> 0),sum(δ .< 0)));
    println()
end

function testternaryembedding(s=8,h=3,)
    println("ternary embedding with s = $s h = $h")
    payload = 0.66*0.1;
    cover=rawview(channelview(load("1.pgm")));
    l=size(cover);
    s=div(s,2);
    cover = cover[div(l[1],2)-s+1:div(l[1],2)+s,div(l[2],2)-s+1:div(l[2],2)+s];

    numofcols=Int(round(1/payload));
    hhat=rand(1:3^h-1,numofcols);
    numofblocks=Int(floor(length(cover)/numofcols));
    message=rand(0:2,numofblocks);
    embpath = randperm(length(cover));

    stegos = SUniward.IncrementalSUniward(cover,3^h;sigma=1.0,T=Float16);
    @time (stego,distortion)=STC.variablecoding(cover,stegos,message,hhat,h,embpath,3);
    δ = Int.(cover) - Int.(stego.image)
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego.image),sum(δ .> 0),sum(δ .< 0)));

    stegos = [SUniward.SUniwardAdd(cover) for i in 1:3^h];
    @time (stego,distortion)=STC.variablecoding(cover,stegos,message,hhat,h,embpath,3);
    δ = Int.(cover) - Int.(stego.image);
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego.image),sum(δ .> 0),sum(δ .< 0)));
    println()
end

function testsaturation()
    println("testing saturation")
    for k in [1,2,254,255]
        h = 3;
        cover = UInt8.(k*ones(8,8))
        payload = 0.66*0.4;
        numofcols=Int(round(1/payload));
        hhat=rand(1:3^h-1,numofcols);
        numofblocks=Int(floor(length(cover)/numofcols));
        message=rand(0:2,numofblocks);
        embpath = randperm(length(cover));
        stegos = [SUniward.SUniwardAdd(cover) for i in 1:3^h];
        @time (stego,distortion)=STC.variablecoding(cover,stegos,message,hhat,h,embpath,3);
        δ = Int.(cover) - Int.(stego.image)
        println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego.image),sum(δ .> 0),sum(δ .< 0)));
    end
    println()
end

function testadditiveembedding(h=3,s=8,p=0.3)
    payload = p;
    cover=rawview(channelview(load("1.pgm")));
    l=size(cover);
    s=div(s,2);
    cover = cover[div(l[1],2)-s+1:div(l[1],2)+s,div(l[2],2)-s+1:div(l[2],2)+s];

    numofcols=Int(round(1/payload));
    numofblocks=Int(floor(length(cover)/numofcols));
    message=rand(0:1,numofblocks);
    embpath = randperm(length(cover));
    hhat=rand(1:2^h-1,numofcols);
    ρ = SUniward.suniwardcosts(Float64.(cover))
    costfun(img,j) = ((mod(img[j],2) == 0)? (0.0,ρ[j]) : (ρ[j],0.0))
    stego,distortion = STC.additivecoding(cover,message,hhat,h,costfun,embpath);
    δ = Int.(cover) - Int.(stego)
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego),sum(δ .> 0),sum(δ .< 0)));

    stegos = [SUniward.SUniwardAdd(cover) for i in 1:2^h]
    @time (stego,distortion)=STC.variablecoding(cover,stegos,message,hhat,h,embpath,2);
    δ = Int.(cover) - Int.(stego.image)
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",distortion,SUniward.suniwardcosts(cover,stego.image),sum(δ .> 0),sum(δ .< 0)));
end

testsuniwardcosts()
testreflectedsetindex()

@testset begin
    @test testoneitemconv2_1([1 2 3; 3 4 5; 6 7 8])<1e-6
    @test testoneitemconv2_1([1 2; 3 4; 6 7])<1e-6
    @test testoneitemconv2_1([1 2 3; 3 4 5])<1e-6
    @test testoneitemconv2_1(randn(16,16))<1e-6
end

@testset begin
    @test testoneitemconv2_2([1 2 3; 3 4 5; 6 7 8])<1e-6
    @test testoneitemconv2_2([1 2; 3 4; 6 7])<1e-6
    @test testoneitemconv2_2([1 2 3; 3 4 5])<1e-6
    @test testoneitemconv2_2(randn(16,16))<1e-6
end

@testset begin
    @test testweightedconv([1 2 3; 3 4 5; 6 7 8])<1e-6
    @test testweightedconv([1 2; 3 4; 6 7])<1e-6
    @test testweightedconv([1 2 3; 3 4 5])<1e-6
    @test testweightedconv(randn(16,16))<1e-6
end

@testset begin
    @test testincrementalsuniward()<1e-6
    @test testLSBcostfun()<1e-6
    @test testtrypmone()<1e-6
end

testbinaryembedding(32,3)
testternaryembedding(32,3)
testsaturation()


text:
abstract Evaluator

evaluation_methods = [:naive, :compiled, :flattened]

immutable Leaf
    value::Float64
end

immutable Inner
    feature::Int64
    splitValue::Float64
    left::Union(Inner, Leaf)
    right::Union(Inner, Leaf)
end

immutable Forest <: Evaluator
    trees::Vector{Inner}
end

evaluate(l::Leaf, features::Vector{Float64}) = l.value

function evaluate(n::Inner, features::Vector{Float64})
    return features[n.feature] < n.splitValue ? evaluate(n.left, features) : evaluate(n.right, features)
end

function evaluate(f::Forest, features::Vector{Float64})
    result = 0.0
    for tree in f.trees
        result += evaluate(tree, features)
    end
    return result
end


text:
"""
perfect_square(number)

Check if a number is a perfect square or not.

# Example
```jula
perfect_square(9)   # returns True
perfect_square(16)  # returns True
perfect_square(1)   # returns True
perfect_square(0)   # returns True
perfect_square(10)  # returns False
perfect_square(-9)  # returns False
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee) and [Rratic](https://github.com/Rratic)
"""
function perfect_square(number::N)where N<:Integer
    number<0&&return false
    val=isqrt(number)
    return val*val == number
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
function bar2line(opt::Dict)
    for i = 1:length(opt["series"])
        opt["series"][i]["type"] = "line"
    end

    opt["series"][1]["boundaryGap"] = false

    return opt
end

function ec_line{S<:String}(legendlist::Vector{S}, data::Array, labels=[1:size(data, 1)]; smooth::Bool = true)
    opt = ec_bar(legendlist, data, labels; smooth = smooth)
    bar2line(opt)
end

function ec_line_piled{S<:String}(legendlist::Vector{S}, data::Array, labels=[1:size(data, 1)]; smooth::Bool = true)
    opt = ec_bar_piled(legendlist, data, labels; smooth = smooth)
    bar2line(opt)
end

function ec_line_h{S<:String}(legendlist::Vector{S}, data::Array, labels=[1:size(data, 1)]; smooth::Bool = true)
    opt = ec_bar_h(legendlist, data, labels; smooth = smooth)
    bar2line(opt)
end

function ec_area{S<:String}(legendlist::Vector{S}, data::Array, labels=[1:size(data, 1)]; smooth::Bool = true)
    opt = ec_line(legendlist, data, labels; smooth = smooth)

    for i = 1:length(opt["series"])
        opt["series"][i]["itemStyle"] = Dict(
                "normal" => Dict("areaStyle" => Dict("type" => "default"))
            )
    end

    return opt
end

function ec_area_piled{S<:String}(legendlist::Vector{S}, data::Array, labels=[1:size(data, 1)]; smooth::Bool = true)
    opt = ec_area(legendlist, data, labels; smooth = smooth)

    for i = 1:length(opt["series"])
        opt["series"][i]["stack"] = "total"
    end

    return opt
end


text:
const LIBAMINDA = Base.source_dir() * "/libaminda.so"

function update_labeling(nits, ρ, σ, ϕ, d, b)
  (nx, ny, nz) = size(ϕ)
  ccall((:update_labeling, LIBAMINDA),
    Void, (Int32, Int32, Int32, Int32, Ptr{Float16}, Ptr{Float16}, Ptr{Float16}, Ptr{Float16}, Ptr{Float16}),
    nits, nx, ny, nz, ρ, σ, ϕ, d, b)
end

function map_to_plane(img, hom, nx, ny)
  @assert size(hom) == (3,3)
  (nc, inx, iny) = size(img)
  tgt = zeros(Float32, (nc, nx, ny))
  ccall((:map_to_plane, LIBAMINDA),
    Void, (Ptr{Float32}, Int32, Int32, Int32, Ptr{Float32}, Int32, Int32, Ptr{Float32}),
    hom', nc, inx, iny, img, nx, ny, tgt)
  return tgt
end

function mean_and_inverse_deviation(img, w)
  (nc, nx, ny) = size(img)
  mn = zeros(Float32, (nc, nx, ny))
  ccall((:mean_and_inverse_deviation, LIBAMINDA),
    Void, (Int32, Int32, Int32, Int32, Ptr{Float32}, Ptr{Float32}),
    w, nc, nx, ny, img, mn)
  return mn
end

function normalized_cross_correlation(im1, mn1, im2, mn2, w)
  (nc, nx, ny) = size(im1)
  nxcorr = fill(nx, ny, -1.0)
  ccall((:normalized_cross_correlation, LIBAMINDA),
    Void, (Int32, Int32, Int32, Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
    w, nc, nx, ny, im1, mn1, im2, mn2, nxcorr)
  return nxcorr
end

function maximum(a,b)
  (nx, ny) = size(a)
  mx = Array(Float32, (nx, ny))
  ccall((:maximum, LIBAMINDA),
    Void, (Int32, Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
    nx, ny, a, b, mx)
  return mx
end

function fill(nx, ny, x)
  a = Array(Float32, (nx, ny))
  ccall((:fill, LIBAMINDA),
    Void, (Int32, Int32, Ptr{Float32}, Float32),
    nx, ny, a, x)
  return a
end

function update_depth(nxc, mxnxc, zs, z)
  (nx, ny) = size(nxc)
  ccall((:update_depth, LIBAMINDA),
    Void, (Int32, Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Float32),
    nx, ny, nxc, mxnxc, zs, z)
end

max_stars_repo_name: twadleigh/aminda

max_issues_repo_name: twadleigh/aminda

max_forks_repo_name: twadleigh/aminda


text:
using Nextion
using Test

@testset "text" begin

    println("Create NexSerial")
    nexSerial = NexSerial("/dev/ttyUSB0")

    println("Init")
    init(nexSerial)

    println("Create objects")
    #nexText = NexText(nexSerial, Name("t1"), pid=PageID(2), cid=ComponentID(1))
    #nexText = NexText(nexSerial, Name("t1"))
    nexText = NexText(nexSerial, Name("t1"), cid=ComponentID(1))

    println("Reset")
    reset(nexSerial)

    sleep(1)

    send(nexSerial, "page pg_text")

    sleep(1)

    nexText.text = "Hello"

    sleep(1)
    @test nexText.text == "Hello"  # ToDo get text

    #=
    sleep(0.5)

    nexText.backcolor = Colour.BLUE
    nexText.text = "1"

    sleep(0.5)

    nexText.backcolor = Colour.WHITE
    nexText.text = "2"

    sleep(0.5)

    nexText.backcolor = Colour.RED
    nexText.text = "3"

    sleep(0.5)

    nexText.text = "Bye nexText"

    sleep(0.5)

    nexText.visible = false
    nexText.backcolor = Colour.WHITE

    sleep(2)

    nexText.visible = true

    sleep(0.5)

    nexText.alignment.horizontal = Alignment.Horizontal.RIGHT
    nexText.alignment.vertical = Alignment.Vertical.DOWN
    nexText.forecolor = Colour.BLUE

    # Change fontid from 0 to 1
    nexText.font = Font(1)

    #reset(nexSerial)
    #sleep(0.5)
    =#

    close(nexSerial)

    @test true

end


text:
using Api
Api.main()


text:
using Documenter
using Discord

makedocs(;
    modules=[Discord],
    format=:html,
    pages=[
        "Home"     => "index.md",
        "Client"   => "client.md",
        "Events"   => "events.md",
        "REST API" => "rest.md",
        "Helpers"  => "helpers.md",
        "Types"    => "types.md",
        "Tutorial" => "tutorial.md",
    ],
    repo="https://github.com/PurgePJ/Discord.jl/blob/{commit}{path}#L{line}",
    sitename="Discord.jl",
    authors="PurgePJ <sindur.esl@gmail.com>, christopher-dG <chrisadegraaf@gmail.com>",
    assets=[
        "assets/logo.png",
    ],
)

deploydocs(;
    repo="github.com/PurgePJ/Discord.jl",
    target="build",
    julia="1.0",
    deps=nothing,
    make=nothing,
)


text:
using GeoPhyInv
using LinearAlgebra
using Random

mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

tgrid=range(0, stop=1.0, step=0.1);

srcwav=SrcWav(tgrid, ageom, [:P]);

Random.randn!(srcwav[3][:P]);

x=randn(length(tgrid));
update!(srcwav, [:P,], x);

x1=randn(length(tgrid));
x2=randn(length(tgrid));
update!(srcwav[1], [:P,], x1);
update!(srcwav[2], [:P,], x2);

rmul!(srcwav, 2.0);


text:
module ReservoirManagement

  using TimeSeries
  using Base.Dates
  using DataStructures
  using SIUnits
  using SIUnits.ShortUnits

  using JuMP
  using Gurobi

  using StatsBase
  using HypothesisTests
  using Distributions
  using ExtremeValueDistributions # Strange dependency (not published in METADATA), try to reimplement somewhere
  # (and make it fully compatible with Julia 0.4, without warnings).

  using DataFrames
  using Colors
  using Plots
  using StatPlots

  # The generalised distributions are defined in both Distributions and ExtremeValueDistributions, but those are identical.
  GeneralizedExtremeValue = Distributions.GeneralizedExtremeValue
  GeneralizedPareto = Distributions.GeneralizedPareto

  # Conflict between standard Dates and SIUnits.
  Second = Dates.Second

  # To allow overriding those functions.
  import Base: copy, isempty, start, next, done, eltype, length

  # Order of inclusion is important!
  include("utils.jl")
  include("reservoir/_ds.jl")
  include("read.jl")
  include("stats/_stats.jl")
  include("scenariosgeneration.jl")
  include("solvers/_solvers.jl")
  include("plots_dataframe.jl")
  include("plots_plots.jl")
  include("report.jl")

  # Export things that should be exported.
  export
    ## utils.jl
    Height, Volume, Discharge, Power,
    SeasonsSeparator, HydrologicalWesternEurope, seasons, isSeasonKnown, ensureSeasonKnown, begins, ends, length, seasonBoundariesGuarantee, separate,
    NonLeapYearIterator,
    ArrayAsIterableReservoirScenarios, eachscenario,
    ArrayAsIterableReservoirScenario, totalInflow,
    toLength,

    ## reservoir/
    # bathymetry.jl
    ReservoirBathymetry, NullReservoirBathymetry, computeVolume, volume, computeLevel, level,
    WhiteBoxReservoirBathymetry, iswhitebox, isblackbox,
    PolynomialReservoirBathymetry, ispolynomial, isPolynomial, coefficient,
    LinearReservoirBathymetry, getIntercept, getSlope, islinear, isLinear, intercept, slope,
    QuadraticReservoirBathymetry, isquadratic, isQuadratic,
    PiecewiseLinearReservoirBathymetry, ispiecewiselinear, isPiecewiseLinear, isPWL,
    SquareRootReservoirBathymetry, issquareroot, isSquareRoot,

    # hydropower.jl
    HydropowerUnit,
    getId, getPurposeId, getDamOutputId, getMaximumDischarge, getMaximumPower, getEfficiency,
    id, purposeId, damOutputId, maxDischarge, maxPower, efficiency,

    # purpose.jl
    Purpose, getId, isDeterministic, id,
    WaterWithdrawalPurpose, isWaterWithdrawal,
    DeterministicWaterWithdrawalPurpose, getNeed, need, DeterministicDrinkingWater, DeterministicEnvironmentalFlow,
    DeterministicHydropower, getHydropowerUnit, isHydropower, hydropower,
    DeterministicPurposes,

    # output.jl
    DamOutput, getId, getCondition, id, condition, hasCondition, isConditional,
    ConveyanceDamOutput, hasConveyanceFactor, getConveyanceFactor, conveyance,
    NoCondition, isempty,
    ConstantDamOutput, getDischarge, isConstant, discharge,
    HydropowerDamOutput, getId, getDischarge, isHydropower, getHydropowerUnit, hydropower,
    ConveyanceSpillway, ConveyancePipe,
    MinimumReservoirLevelCondition, getMinimumLevel, minLevel,
    ConstantDamOutputs,

    # river.jl
    River, getName, getScenarios, getScenarioTS, getScenario,
    countScenarios, countTimeSteps, getPeriod, getTimePoints,
    name, scenarios, scenario, period,

    # river_natural.jl
    NaturalRiver,

    # river_diverted.jl
    DivertedRiver,
    getEnvironmentalFlow, getMaximumFlow, getMaximumAllowableFlow,
    environmentalFlow, maxFlow, maxAllowableFlow,

    # reservoir.jl
    Reservoir, getName, getVariant, getMinimumCapacity, getMaximumCapacity, getHydropowerUnits,
    hasHydropower, getHydropowerIds, getPenstockHydropower,
    getPurposes, getOutputs, getBathymetry, hasPurpose, getPurpose, getPurposeIds, getOutputIds, hasOutput, getOutput,
    name, variant, minCapacity, maxCapacity, penstockHydropower, hydropowerUnits, outputs, outputsId, output, bathymetry, purpose,
    countScenarios, countTimeSteps, getPeriod, getTimePoints,
    rivers, naturalRivers, divertedRivers, getNaturalRivers, getDivertedRivers, getRivers, getNaturalRiver, getDivertedRiver, hasRiver, getRiver,
    getDivertedRiverNames, getNaturalRiverNames, getRiverNames,
    totalInflow, removeScenarios, keepTimeSteps,
    NaturalRiverIterator, eachnaturalriver,
    DivertedRiverIterator, eachdivertedriver,
    RiverIterator, eachriver,
    ScenarioReservoirIterator, eachscenario,

    # basin.jl
    # Not for now: nothing uses it. For future multireservoir implementation.

    # bathymetry_fwd.jl
    level,

    # output_fwd.jl
    getMinimumVolume, minVolume,

    # misc.jl
    countScenarios, countTimeSteps,

    ## read.jl
    makeScenarios,

    ## stats/
    # base.jl
    computeCI,
    RiverScenarioIterator, eachscenario,
    generateScenarios,
    ReservoirScenarioIterator, eachscenario,

    # utils.jl
    splitIntoSets,

    # adam2015.jl
    printReport,

    # tstudent.jl
    # Nothing, no user-visible interface (only through the functions of base.jl).

    ## scenariosgeneration.jl
    scenarioDuplication, scenarioConcatenation, scenarioMerging, scenarioMixing,
    addScenarioGenerationAlgorithm!,
    scenarioGeneration, scenarioDuplication, scenarioConcatenation, scenarioConcatenation, scenarioMerging, scenarioMixing,
    scenarioShift,

    ## solvers/
    # ds.jl
    Options, Solution, SolverStatistics, ReservoirSolution,
    RuleCurveOptions, RuleCurveSolution, EmptyRuleCurveSolution, RuleCurveSolverStatistics,
    MinimumRuleCurveOptions, MinimumRuleCurveSolution, EmptyMinimumRuleCurveSolution, MinimumRuleCurveSolverStatistics,
    MaximumRuleCurveOptions, MaximumRuleCurveSolution, EmptyMaximumRuleCurveSolution, MaximumRuleCurveSolverStatistics,
    EventOptions, EventSolution, EmptyEventSolution, EventSolverStatistics,
    InverseOptions, InverseSolution, EmptyInverseSolution, InverseSolverStatistics,
    EvaluationOptions, EvaluationSolution, EmptyEvaluationSolution, EmptyEvaluationSolution, EvaluationSolverStatistics,

    # ds_accessors.jl
    getReservoir, isFeasible, getSolutionObject, hasName, getName, getSolverIdentifier, getSolverDescription, getSolverOptions, getSolverStatistics,
    isInfeasible, getNameOrElse, getNameOrElse,
    getPeriod, getVariant,
    getSolution, countTimeSteps,

    # none.jl
    NoSolverOptions, NoSolverSolution, NoSolverStatistics, wrapExistingSolution,

    # _helpers_inflow.jl
    # _helpers_output.jl
    # _helpers_bathymetry.jl
    # _helpers_purpose.jl
    # For now, helpers API not old enough to be exposed to the public.

    # evaluate_crossvalidate.jl
    FeasibilityCrossValidationOptions, FeasibilityCrossValidationSolution, FeasibilityCrossValidationSolverStatistics,
    crossValidatedFeasibilityProbability,

    # evaluate_deviation.jl
    DeviationOptions, DeviationSolution, DeviationSolverStatistics,
    maximumDeviation,

    # evaluate_mc_feasibility.jl
    FeasibilityMonteCarloOptions, FeasibilityMonteCarloSolution, FeasibilityMonteCarloSolverStatistics,
    feasibilityProbability,

    # evaluate_purposeshortage.jl
    ShortageOptions, ShortageSolution, ShortageSolverStatistics,
    computePurposeShortageIndex,

    # evaluate_hydropower.jl
    # Nothing for now, far from working state.

    # inverse_functions.jl
    FunctionsInverseOptions, FunctionsInverseSolution, FunctionsInverseSolverStatistics,
    worstCaseFunctions,

    # rulecurve_min_default.jl
    DefaultMinimumRuleCurveOptions, DefaultMinimumRuleCurveSolution, DefaultMinimumRuleCurveSolverStatistics,
    minimumRuleCurve,

    # rulecurve_min_safe.jl
    SafetyMinimumRuleCurveOptions, SafetyMinimumRuleCurveSolution, SafetyMinimumRuleCurveSolverStatistics,
    safeMinimumRuleCurve,

    ## plots_
    # dataframes.jl
    getDataFrameForRiver, getDataFrameForSolutions,

    # plots.jl
    plotRiverAlongReservoirs, plotSolutionVsScenarios, plotSolutions, plotRule,

    ## report.jl
    report
end


text:
using Documenter

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

using Turkie

makedocs(;
    sitename="Turkie",
    format=Documenter.HTML(),
    modules=[Turkie],
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/theogf/Turkie.jl.git", push_preview=true
)


text:
module Test_rmapps

using Random
using LinearAlgebra
using MatrixPencils
using Polynomials
using Test

# λ = Polynomial([0,1],:λ)
# @time info, iz, nfp, ip = rmkstruct(λ, 1, fast = true)
Random.seed!(21235)

@testset "Rational Matrix Applications" begin

fast = true
for fast in (false,true)

# P = 0

@time info, iz, nfp, ip = rmkstruct(0, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [0], [], 0, 0) && iz == [] && ip == [] && nfp == 0

@time info, iz, nfp, ip = rmkstruct(0, 1, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [0], [], 0, 0) && iz == [] && ip == [] && nfp == 0

@time val, iz, info = rmzeros(0, fast = fast)
@test val ≈ Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [0], [], 0, 0)

@time val, iz, info = rmzeros(0, 1, fast = fast)
@test val ≈ Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [0], [], 0, 0)

@time val, ip, id = rmpoles(0, fast = fast)
@test val ≈ Float64[] && ip == [] && id == []

@time val, ip, id = rmpoles(0, 1, fast = fast)
@test val ≈ Float64[] && ip == [] && id == []

@test rmrank(0) == 0 && rmrank(0,fastrank=false) == 0

@test rmrank(0,1) == 0 && rmrank(0,1,fastrank=false) == 0

# P = 1
@time info, iz, nfp, ip = rmkstruct(1, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [], 0, 1) && iz == [] && ip == [] && nfp == 0

@time info, iz, nfp, ip = rmkstruct(1, 1, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [1], 0, 1) && iz == [] && ip == [] && nfp == 0

@time val, iz, info = rmzeros(1, fast = fast)
@test val ≈ Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [1], 0, 1)

@time val, iz, info = rmzeros(1, 1, fast = fast)
@test val ≈ Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [1], 0, 1)

@time val, ip, id = rmpoles(1, fast = fast)
@test val ≈ Float64[] && ip == [] && id == []

@time val, ip, id = rmpoles(1, 1, fast = fast)
@test val ≈ Float64[] && ip == [] && id == []

@test rmrank(1) == 1 && rmrank(1,fastrank=false) == 1

@test rmrank(1,1) == 1 && rmrank(1,1,fastrank=false) == 1

# P = λ
λ = Polynomial([0,1],:λ)

@time info, iz, nfp, ip = rmkstruct(λ, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [], 1, 1) && iz == [] && ip == [1] && nfp == 0

@time info, iz, nfp, ip = rmkstruct(λ, 1, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [1, 1], 1, 1) && iz == [] && ip == [1] && nfp == 0

@time info, iz, nfp, ip = rmkstruct(1, λ, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [2], 0, 1) && iz == [1] && ip == [] && nfp == 1

@time val, iz, info = rmzeros(λ, fast = fast)
@test val ≈ [0.] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [1,1], 1, 1)

@time val, iz, info = rmzeros(λ, 1, fast = fast)
@test val ≈ [0.] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [1,1], 1, 1)

@time val, iz, info = rmzeros(1, λ, fast = fast)
@test val == [Inf] && iz == [1] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [2], 0, 1)

@time val, ip, id = rmpoles(λ, fast = fast)
@test val == [Inf] && ip == [1] && id == [2]

@time val, ip, id = rmpoles(λ, 1, fast = fast)
@test val == [Inf] && ip == [1] && id == [2]

@test rmrank(λ) == 1 && rmrank(λ,fastrank=false) == 1

@test rmrank(λ,1) == 1 && rmrank(λ,1,fastrank=false) == 1

# P = [λ 1]
λ = Polynomial([0,1],:λ)
N = [λ one(λ)]
D = [one(λ) one(λ)]

@time  info, iz, nfp, ip = rmkstruct(N, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [], [], 0, 1) && iz == [] && ip == [1] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(D, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [], [], 0, 1) && iz == [] && ip == [] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(N, D, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [], [1, 1], 0, 1) && iz == [] && ip == [1] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(D, N, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [], [1], 0, 1) && iz == [] && ip == [] && nfp == 1

@time val, iz, info = rmzeros(N, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [], [1, 1], 0, 1)

@time val, iz, info = rmzeros(N, D, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [], [1, 1], 0, 1)

@time val, iz, info = rmzeros(D, N, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [], [1], 0, 1)

@time val, ip, id = rmpoles(N, fast = fast)
@test val == [Inf] && ip == [1] && id == [2]

@time val, ip, id = rmpoles(N, D, fast = fast)
@test val == [Inf] && ip == [1] && id == [2]

@time val, ip, id = rmpoles(D, N, fast = fast)
@test val ≈ [0] && ip == [] && id == []

@test rmrank(N) == 1 && rmrank(N,fastrank=false) == 1

@test rmrank(N,D) == 1 && rmrank(N,D,fastrank=false) == 1

@test rmrank(D,N) == 1 && rmrank(D,N,fastrank=false) == 1

# Example 3: P = [λ^2 λ; λ 1] DeTeran, Dopico, Mackey, ELA 2009
λ = Polynomial([0,1],:λ)
N = [λ^2 λ; λ 1]
D = [1 1; 1 1]

@time  info, iz, nfp, ip = rmkstruct(N, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [], 0, 1) && iz == [] && ip == [2] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(N, D, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [1, 1, 1], 0, 1) && iz == [] && ip == [2] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(poly2pm(N), D, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [1, 1, 1], 0, 1) && iz == [] && ip == [2] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(N, poly2pm(D), fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [1, 1, 1], 0, 1) && iz == [] && ip == [2] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(D, N, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [1], 0, 1) && iz == [] && ip == [] && nfp == 2

@time val, iz, info = rmzeros(N, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [1, 1, 1], 0, 1)

@time val, iz, info = rmzeros(N, D, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [1, 1, 1], 0, 1)

@time val, iz, info = rmzeros1(N, D, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [], 0, 1)

@time val, iz, info = rmzeros(N, poly2pm(D), fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [1, 1, 1], 0, 1)

@time val, iz, info = rmzeros(D, N, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [1], 0, 1)

@time val, ip, id = rmpoles(N, fast = fast)
@test val == [Inf, Inf] && ip == [2] && id ==  [1, 3]

@time val, ip, id = rmpoles(N, D, fast = fast)
@test val == [Inf, Inf] && ip == [2] && id ==  [1, 3]

@time val, ip, id = rmpoles(N, poly2pm(D), fast = fast)
@test val == [Inf, Inf] && ip == [2] && id ==  [1, 3]

@time val, ip, id = rmpoles(D, N, fast = fast)
@test val ≈ [0, 0] && ip == [] && id ==  []

@test rmrank(N) == 1 && rmrank(N, D, fastrank=false) == 1

# Example 4: Van Dooren, Dewilde, LAA, 1983
N = zeros(Int,3,3,3);
N[:,:,1] = [1 2 -2; 0 -1 -2; 0 0 0];
N[:,:,2] = [1 3 0; 1 4 2; 0 -1 -2];
N[:,:,3] = [1 4 2; 0 0 0; 1 4 2];
D = ones(3,3);

@time val, iz, info = rmzeros(N, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 1, 1, 1], 1, 2)

@time val, iz, info = rmzeros(N, D, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 1, 1, 1], 1, 2)

@time val, ip, id = rmpoles(N, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 3]

@time val, ip, id = rmpoles(N, D, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 3]

@time val, iz, info = rmzeros1(N, D, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1], 1, 2)

@time val, iz, info = rmzeros1(N, poly2pm(D), fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1], 1, 2)

@time val, iz, info = rmzeros1(pm2poly(N), reshape(D,3,3,1), fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1], 1, 2)

@time val, iz, info  = rmzeros1(N, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1], 1, 2)

@time val, ip, id = rmpoles1(N, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 1, 1, 1, 3]

@time val, ip, id = rmpoles1(N, D, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 1, 1, 1, 3]

@test rmrank(N) == 2 && rmrank(N, D, fastrank=false) == 2

λ = Polynomial([0,1],:λ);
N = [λ^2 + λ + 1 4λ^2 + 3λ + 2 2λ^2 - 2;
λ 4λ - 1 2λ - 2;
λ^2 4λ^2 - λ 2λ^2 - 2λ];
D = ones(3,3);

@time  info, iz, nfp, ip = rmkstruct(N, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [1], [2], 1, 2)  && iz == [] && ip == [2] && nfp == 0

@time  info, iz, nfp, ip = rmkstruct(N, D, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [1], [1, 1, 1, 1], 1, 2)  && iz == [] && ip == [2] && nfp == 0

@time info, iz, nfp, ip = rmkstruct(D, N, fast = fast, atol = 1.e-7)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([3], [1], [2, 5], 0, 2)  && iz == [1, 4] && ip == [] && nfp == 9

@time val, iz, info = rmzeros(N, fast = fast, atol = 1.e-7)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 1, 1, 1], 1, 2)

@time val, iz, info = rmzeros(N, D, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 1, 1, 1], 1, 2)

@time val, iz, id = rmpoles(N, fast = fast)
@test val ≈ [Inf, Inf] && iz == [2] && id == [1, 3]

@time val, ip, id = rmpoles(N, D, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 3]

@time val, iz, info = rmzeros1(N, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1], 1, 2)

@time val, ip, id = rmpoles1(N, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 1, 1, 1, 3]

@test rmrank(N) == 2 && rmrank(N,D,fastrank=false) == 2

@test rmrank(rand(2,3,5)) == 2 && rmrank(rand(2,3,5),rand(2,3,5),fastrank=false) == 2

@test rmrank(rand(3,2,5)) == 2 && rmrank(rand(3,2,5),rand(3,2,5),fastrank=false) == 2

# Example Vlad Ionescu
s = Polynomial([0, 1],:s);
N = [s^2+3*s+3 1; -1 2*s^2+7*s+4];
D = [(s+1)^2 s+2; (s+1)^3 (s+1)*(s+2)]
@time info, iz, nfp, ip = rmkstruct(N, D, fast = fast, atol=1.e-7)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == (Int64[], Int64[], [1, 1], 4, 2)  && iz == [] && ip == [] && nfp == 4

@time val, ip, id = rmpoles(N, D, fast = fast, atol=1.e-7)
@test fromroots(val) ≈ Polynomial([2, 7, 9, 5, 1]) && ip == [] && id == []

@time val, ip, id = rmpoles1(N, D, fast = fast, atol=1.e-7)
@test fromroots(val) ≈ Polynomial([2, 7, 9, 5, 1]) && ip == [] && id == [1, 1, 1, 1]

@time val, ip, id = rmpoles1(poly2pm(N), D, fast = fast, atol=1.e-7)
@test fromroots(val) ≈ Polynomial([2, 7, 9, 5, 1]) && ip == [] && id == [1, 1, 1, 1]

@time val, ip, id = rmpoles1(N, poly2pm(D), fast = fast, atol=1.e-7)
@test fromroots(val) ≈ Polynomial([2, 7, 9, 5, 1]) && ip == [] && id == [1, 1, 1, 1]

# Example 1 (Polynomial): Dopico et al. 2020
c1 = [0.06155268911372547, 0.006545378521362721, 0.40445990039119284, 0.8829892254580274, 0.7573496766341161, 0.4094804382958148];
c2 = [0.4059208742974014, 0.48094667571705574];
e1 = Polynomial(c1)
e2 = Polynomial(c2)
P = [e1 0; 0 e2]

@time zer, iz, info = rmzeros(P,fast = fast, atol=1.e-9)
@test fromroots(zer) ≈ e1*e2/(e1[end]*e2[end]) && iz == [] && (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [1, 1, 1, 1], 6, 2)

# Example 2 (Rational): Dopico et al. 2020
c1 = [0.06155268911372547, 0.006545378521362721, 0.40445990039119284, 0.8829892254580274, 0.7573496766341161, 0.4094804382958148];
c2 = [0.4059208742974014, 0.48094667571705574];
e1 = Polynomial(c1)
e2 = Polynomial(c2)
N = [e1 0; 1 e2]
D = [1 1; Polynomial([0, 1]) 1]

@time zer, iz, info = rmzeros(N, D, fast = fast, atol=1.e-9)
@test fromroots(zer) ≈ e1*e2/(e1[end]*e2[end])*Polynomial([0,1]) && iz == [] && (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [1, 1, 1, 1], 7, 2)

@time pol, ip, id = rmpoles(N, D, fast = fast, atol=1.e-7)
@test pol ≈ [0, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [1, 5] && id == [2, 6]

# Example 4.26 (Rational): Dopico et al. 2020 (submitted)
s = Polynomial([0,1]);
N = [s^2+s-1 -1; -1 -s^2+s-2];
D = [s s; 1 1];
@time zer, iz, info = rmzeros(N, D, fast = fast, atol=1.e-7)
@test fromroots(zer) ≈ Polynomial([-1, 3, 0, 0, 1]) && iz == [] && (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [1, 1, 1, 1], 4, 2)

@time pol, ip, id = rmpoles(N, D, fast = fast, atol=1.e-7)
@test pol ≈ [0, Inf, Inf, Inf] && ip == [1, 2] && id == [2, 3]

Ty = Float64
for Ty in (Float64, Complex{Float64})

N = rand(Ty,2,3,5);
N1 = rand(Ty,3,2,5);
D = ones(2,3);
D1 = ones(3,2);

abstol = sqrt(eps(one(real(Ty))))

@time info, iz, nfp, ip = rmkstruct(N, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([8], [], [], 0, 2) && iz == [] && ip ==[4, 4] && nfp == 0

@time info, iz, nfp, ip = rmkstruct(N1, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [8], [], 0, 2) && iz == [] && ip ==[4, 4] && nfp == 0

@time info, iz, nfp, ip = rmkstruct(N, N, fast = fast, atol = abstol)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0, 0], [0], [1], 0, 1) && iz == [] && ip ==[] && nfp == 0

@time val, iz, info = rmzeros(N, fast = fast)
@test val == Ty[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [1, 1, 1, 1], 0, 2)

@time val, ip, id = rmpoles(N, fast = fast)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [5, 5]

@time val, iz, info = rmzeros1(N, fast = fast)
@test val == Ty[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [], 0, 2)

@time val, iz, info = rmzeros1(N1, fast = fast)
@test val == Ty[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [8], [], 0, 2)

@time val, ip, id = rmpoles1(N, fast = fast, atol = abstol)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [1, 5, 5]

@time val, ip, id = rmpoles1(N1, fast = fast, atol = abstol)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [1, 5, 5]

@test rmrank(N) == 2 && rmrank(N,fastrank=false) == 2
@test rmrank(N,D) == 2 && rmrank(N,D,fastrank=false) == 2

@test rmrank(N1) == 2 && rmrank(N1,fastrank=false) == 2
@test rmrank(N1,D1) == 2 && rmrank(N1,D1,fastrank=false) == 2

N = rand(Ty,2,2,4);
D = rand(Ty,2,2,2);

@time info, iz, nfp, ip = rmkstruct(N, D, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [1, 1, 1, 1], 8, 2) && iz == [] && ip ==[2, 2] && nfp == 4

@time val, iz, info = rmzeros(N, D, fast = fast)
@test iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [1, 1, 1, 1], 8, 2)

@time val, ip, id = rmpoles(N, D, fast = fast)
@test real(val[isinf.(val)]) == Ty[Inf, Inf, Inf, Inf] && ip == [2, 2] && id == [3, 3]

@time info, iz, nfp, ip = rmkstruct(D, N, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [3, 3], 8, 2) && iz == [2, 2] && ip ==[] && nfp == 12

@time val, iz, info = rmzeros(D, N, fast = fast)
@test real(val[isinf.(val)]) == Ty[Inf, Inf, Inf, Inf] && iz == [2, 2] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [3, 3], 8, 2)

@time val, ip, id = rmpoles(D, N, fast = fast)
@test length(val) == 12 && ip == [] && id == []

end
end
#end

end  # applications testset

end  # module


text:
mutable struct Experiment
	id::Int64
	title::String
  date::String # e.g. "20211121"
  category::String # e.g "Running"
  rating::Int64 # e.g "0"
  userid::Int64 # e.g "46"
  locked::Int8 # e.g "0"
  canread::String # e.g "team"
  canwrite::String # e.g "user"
  lastchange::DateTime
  up_item_id::Union{Nothing,Any}
  has_attachment::Union{Nothing,Any}
  next_step::Union{Nothing,Any}
  category_id::Int64 # e.g "10"
  color::String # e.g "29AEB9"
  fullname::String # e.g "John Doe"
  recent_comment::Union{Nothing,Any}
  has_comment::Int8 # e.g "0"
  timestamped::Int8 # e.g "0"
end

# we need this to parse json strings as Julia structs
StructTypes.StructType(::Type{Experiment}) = StructTypes.Struct()

Base.show(io::IO, x::Experiment) = print(io,
"""
$(typeof(x)):
id: $(x.id) title: $(x.title)
""")


text:
# TernaryLayer and TernaryMERA types, and methods thereof.
# To be `included` in MERAKit.jl.

"""
    TernaryLayer{ST, ET, Tan} <: SimpleLayer

The type for layers of a ternary MERA.

Each layer consists of two tensors, a 2-to-2 disentangler, often called `u`, and a 3-to-1
isometry, often called `w`.

The type parameters are `ST` for space type, e.g. `ComplexSpace` or `SU2Space`; `ET` for
element type, e.g. `Complex{Float64}`; and `Tan` for whether this layer is a tangent layer
or not.  If `Tan = false`, the layer in question is an actual MERA layer. If `Tan = true` it
consists, instead of the actual tensors, of Stiefel/Grassmann tangent vectors of these
tensors.

`TernaryLayer` implements the iteration interface, returning first the disentangler, then
the isometry.

Index numbering convention is as follows, where the physical indices are at the bottom:
Disentangler:
```
 3│    4│
 ┌┴─────┴┐
 │   u   │
 └┬─────┬┘
 1│    2│
```

Isometry:
```
    4│
 ┌───┴───┐
 │   w   │
 └┬──┬──┬┘
 1│ 2│ 3│
```
"""
struct TernaryLayer{ST, ET, Tan} <: SimpleLayer
    # Even though the types are here marked as any, they are restricted to be specific,
    # concrete types dependent on ST, ET and Tan, both in the constructor and in
    # getproperty.
    disentangler::Any
    isometry::Any

    function TernaryLayer{ST, ET, Tan}(disentangler, isometry) where {ST, ET, Tan}
        DisType = disentangler_type(ST, ET, Tan)
        IsoType = ternaryisometry_type(ST, ET, Tan)
        disconv = convert(DisType, disentangler)::DisType
        isoconv = convert(IsoType, isometry)::IsoType
        return new{ST, ET, Tan}(disconv, isoconv)
    end
end

function TernaryLayer(disentangler::DisType, isometry::IsoType) where {
    ST,
    DisType <: AbstractTensorMap{ST, 2, 2},
    IsoType <: AbstractTensorMap{ST, 3, 1}
}
    ET = eltype(DisType)
    @assert eltype(IsoType) === ET
    return TernaryLayer{ST, ET, false}(disentangler, isometry)
end

function TernaryLayer(disentangler::DisTanType, isometry::IsoTanType) where {
    ST,
    DisType <: AbstractTensorMap{ST, 2, 2},
    IsoType <: AbstractTensorMap{ST, 3, 1},
    DisTanType <: Stiefel.StiefelTangent{DisType},
    IsoTanType <: Grassmann.GrassmannTangent{IsoType},
}
    ET = eltype(DisType)
    @assert eltype(IsoType) === ET
    return TernaryLayer{ST, ET, true}(disentangler, isometry)
end

function Base.getproperty(l::TernaryLayer{ST, ET, Tan}, sym::Symbol) where {ST, ET, Tan}
    if sym === :disentangler
        T = disentangler_type(ST, ET, Tan)
    elseif sym === :isometry
        T = ternaryisometry_type(ST, ET, Tan)
    else
        T = Any
    end
    return getfield(l, sym)::T
end

"""
    TernaryMERA{N}

A ternary MERA is a MERA consisting of `TernaryLayer`s.
"""
TernaryMERA{N} = GenericMERA{N, T, O} where {T <: TernaryLayer, O}
layertype(::Type{TernaryMERA}) = TernaryLayer
#Base.show(io::IO, ::Type{TernaryMERA}) = print(io, "TernaryMERA")
#Base.show(io::IO, ::Type{TernaryMERA{N}}) where {N} = print(io, "TernaryMERA{($N)}")

# Implement the iteration and indexing interfaces.
# See simplelayer.jl for details.
_tuple(layer::TernaryLayer) = (layer.disentangler, layer.isometry)

function operatortype(::Type{TernaryLayer{ST, ET, false}}) where {ST, ET}
    return tensormaptype(ST, 2, 2, ET)
end
operatortype(::Type{TernaryLayer{ST, ET, true}}) where {ST, ET} = Nothing

TensorKit.spacetype(::Type{TernaryLayer{ST, ET, Tan}}) where {ST, ET, Tan} = ST

scalefactor(::Type{<:TernaryLayer}) = 3

causal_cone_width(::Type{<:TernaryLayer}) = 2

Base.eltype(::Type{TernaryLayer{ST, ET, Tan}}) where {ST, ET, Tan} = ET

outputspace(layer::TernaryLayer) = space(layer.disentangler, 1)
inputspace(layer::TernaryLayer) = space(layer.isometry, 4)'
internalspace(layer::TernaryLayer) = space(layer.isometry, 1)
internalspace(m::TernaryMERA, depth) = internalspace(getlayer(m, depth))

function expand_inputspace(layer::TernaryLayer, V_new)
    u, w = layer
    w = pad_with_zeros_to(w, 4 => V_new')
    return TernaryLayer(u, w)
end

function expand_outputspace(layer::TernaryLayer, V_new)
    u, w = layer
    u = pad_with_zeros_to(u, 1 => V_new, 2 => V_new)
    w = pad_with_zeros_to(w, 2 => V_new)
    return TernaryLayer(u, w)
end

function expand_internalspace(layer::TernaryLayer, V_new)
    u, w = layer
    u = pad_with_zeros_to(u, 3 => V_new', 4 => V_new')
    w = pad_with_zeros_to(w, 1 => V_new, 3 => V_new)
    return TernaryLayer(u, w)
end

function randomlayer(
    ::Type{TernaryLayer}, ::Type{T}, Vin, Vout, Vint = Vout; random_disentangler = false
) where {T}
    w = randomisometry(T, Vint ⊗ Vout ⊗ Vint, Vin)
    u = initialize_disentangler(T, Vout, Vint, random_disentangler)
    return TernaryLayer(u, w)
end

function ascending_fixedpoint(layer::TernaryLayer)
    width = causal_cone_width(layer)
    Vtotal = ⊗(ntuple(n->inputspace(layer), Val(width))...)
    return id(storagetype(operatortype(layer)), Vtotal)
end

function scalingoperator_initialguess(l::TernaryLayer, irreps...)
    width = causal_cone_width(l)
    V = inputspace(l)
    interlayer_space = ⊗(ntuple(n->V, Val(width))...)
    outspace = interlayer_space
    inspace = interlayer_space
    for irrep in irreps
        inspace = inspace ⊗ spacetype(V)(irrep => 1)
    end
    typ = eltype(l)
    t = TensorMap(randn, typ, outspace ← inspace)
    return t
end

function gradient(layer::TernaryLayer, env::TernaryLayer; metric = :euclidean)
    u, w = layer
    uenv, wenv = env
    # The environment is the partial derivative. We need to turn that into a tangent vector
    # of the Stiefel manifold point u or w.
    # The factor of two is from the partial_x + i partial_y derivative of the cost function,
    # and how it depends on both v and v^dagger.
    ugrad = Stiefel.project!(2*uenv, u; metric = metric)
    wgrad = Grassmann.project!(2*wenv, w)
    return TernaryLayer(ugrad, wgrad)
end

function precondition_tangent(layer::TernaryLayer, tan::TernaryLayer, rho)
    u, w = layer
    utan, wtan = tan
    @planar rho_wl[-1; -11] := rho[-1 1; -11 1]
    @planar rho_wr[-1; -11] := rho[1 -1; 1 -11]
    rho_w = (rho_wl + rho_wr) / 2.0
    @planar(
        rho_u[-1 -2; -11 -12] :=
        w'[12; 1 2 -11] * w'[22; -12 3 4] *
        rho[11 21; 12 22] *
        w[1 2 -1; 11] * w[-2 3 4; 21]
    )
    utan_prec = precondition_tangent(utan, rho_u)
    wtan_prec = precondition_tangent(wtan, rho_w)
    return TernaryLayer(utan_prec, wtan_prec)
end

# # # Invariants

function space_invar_intralayer(layer::TernaryLayer)
    u, w = layer
    matching_bonds = (
        (space(u, 3)', space(w, 3)),
        (space(u, 4)', space(w, 1)),
    )
    allmatch = all(pair -> ==(pair...), matching_bonds)
    # Check that the dimensions are such that isometricity can hold.
    allmatch &= all((u, w)) do v
        codom, dom = fuse(codomain(v)), fuse(domain(v))
        infimum(dom, codom) == dom
    end
    return allmatch
end

function space_invar_interlayer(layer::TernaryLayer, next_layer::TernaryLayer)
    u, w = layer.disentangler, layer.isometry
    unext, wnext = next_layer.disentangler, next_layer.isometry
    matching_bonds = (
        (space(w, 4)', space(unext, 1)),
        (space(w, 4)', space(unext, 2)),
        (space(w, 4)', space(wnext, 2)),
    )
    allmatch = all(pair -> ==(pair...), matching_bonds)
    return allmatch
end

# # # Ascending and descending superoperators

# TODO Fix this to work with anyons.
#"""
#Return the ascending superoperator of the one site in the middle of the isometries in a
#TernaryMERA, as a TensorMap. Unlike most ascending superoperators, this one is actually
#affordable to construct as a full tensor.
#"""
#ascending_superop_onesite(m::TernaryMERA) = ascending_superop_onesite(getlayer(m, Inf))
#
#function ascending_superop_onesite(layer::TernaryLayer)
#    w = layer.isometry
#    @planar superop[-1 -2; -11 -12] := w'[-11; 1 -1 2] * w[1 -2 2; -12]
#    return superop
#end

const TernaryOperator{S} = AbstractTensorMap{S, 2, 2}
const ChargedTernaryOperator{S} = AbstractTensorMap{S, 2, 3}
const DoubleChargedTernaryOperator{S} = AbstractTensorMap{S, 2, 4}

function ascend(
    op::Union{TernaryOperator, ChargedTernaryOperator, DoubleChargedTernaryOperator},
    layer::TernaryLayer
) where {S1}
    u, w = layer
    l = ascend_left(op, layer)
    r = ascend_right(op, layer)
    m = ascend_mid(op, layer)
    scaled_op = (l+r+m)/3
    return scaled_op
end

function ascend(op::SquareTensorMap{1}, layer::TernaryLayer)
    return ascend(expand_support(op, causal_cone_width(TernaryLayer)), layer)
end

# TODO Think about how to best remove the code duplication of having the separate methods
# for ordinary, charged, and double charged operators. Note also that there's a mix of using
# @tensor and @planar.
function ascend_left(op::ChargedTernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @tensor(
        scaled_op[-100 -200; -300 -400 -1000] :=
        w[51 52 53; -300] * w[54 11 12; -400] *
        u[41 42; 53 54] *
        op[31 32; 52 41 -1000] *
        u'[21 55; 32 42] *
        w'[-100; 51 31 21] * w'[-200; 55 11 12]
    )
    return scaled_op
end

function ascend_right(op::ChargedTernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @tensor(
        scaled_op[-100 -200; -300 -400 -1000] :=
        w[11 12 65; -300] * w[63 61 62; -400] *
        u[51 52; 65 63] *
        op[31 41; 52 61 -1000] *
        u'[64 21; 51 31] *
        w'[-100; 11 12 64] * w'[-200; 21 41 62]
    )
    return scaled_op
end

function ascend_mid(op::ChargedTernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 6X^6
    @tensor(
        scaled_op[-100 -200; -300 -400 -1000] :=
        w[31 32 41; -300] * w[51 21 22; -400] *
        u[1 2; 41 51] *
        op[11 12; 1 2 -1000] *
        u'[42 52; 11 12] *
        w'[-100; 31 32 42] * w'[-200; 52 21 22]
    )
    return scaled_op
end

function ascend_left(op::TernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        scaled_op[-100 -200; -300 -400] :=
        w[51 52 53; -300] * w[54 11 12; -400] *
        u[41 42; 53 54] *
        op[31 32; 52 41] *
        u'[21 55; 32 42] *
        w'[-100; 51 31 21] * w'[-200; 55 11 12]
    )
    return scaled_op
end

function ascend_right(op::TernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        scaled_op[-100 -200; -300 -400] :=
        w[11 12 65; -300] * w[63 61 62; -400] *
        u[51 52; 65 63] *
        op[31 41; 52 61] *
        u'[64 21; 51 31] *
        w'[-100; 11 12 64] * w'[-200; 21 41 62]
    )
    return scaled_op
end

function ascend_mid(op::TernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 6X^6
    @planar(
        scaled_op[-100 -200; -300 -400] :=
        w[31 32 41; -300] * w[51 21 22; -400] *
        u[1 2; 41 51] *
        op[11 12; 1 2] *
        u'[42 52; 11 12] *
        w'[-100; 31 32 42] * w'[-200; 52 21 22]
    )
    return scaled_op
end

function ascend_left(op::DoubleChargedTernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        temp1[-1 -2 -3 -4; -11 -12 -13 -14] :=
        op[31 32; -1 -11 -12 -13] *
        u'[21 -4; 32 -14] *
        w'[-3; -2 31 21]
    )
    temp2 = braid(temp1, (11, 12, 13, 14, 15, 1, 100, 16), (1, 2, 3, 7, 6, 4), (5, 8))
    @planar(
        temp3[-100 -1000 -2000 -200; -300 -400] :=
        w[51 52 53; -300] * w[54 11 12; -400] *
        u[41 42; 53 54] *
        temp2[52 51 -100 -1000 -2000 55; 41 42] *
        w'[-200; 55 11 12]
    )
    scaled_op = braid(temp3, (11, 100, 1, 12, 13, 14), (1, 4), (5, 6, 3, 2))
    return scaled_op
end

function ascend_right(op::DoubleChargedTernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        temp1[-1 -2 -3 -4 -5 -6; -11 -12] :=
        op[31 41; -12 -6 -5 -4] *
        u'[-1 21; -11 31] *
        w'[-2; 21 41 -3]
    )
    temp2 = braid(temp1, (11, 12, 13, 100, 1, 14, 15, 16), (1, 2, 4, 5, 3, 6), (7, 8))
    @planar(
        temp3[-100 -200 -1000 -2000; -300 -400] :=
        w[11 12 65; -300] * w[63 61 62; -400] *
        u[51 52; 65 63] *
        temp2[64 -200 -1000 -2000 62 61; 51 52] *
        w'[-100; 11 12 64]
    )
    scaled_op = braid(temp3, (11, 12, 100, 1, 13, 14), (1, 2), (5, 6, 4, 3))
    return scaled_op
end

function ascend_mid(op::DoubleChargedTernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 6X^6
    @planar(
        temp1[-1 -2; -3 -4 -5 -6] :=
        u[1 2; -3 -4] *
        op[11 12; 1 2 -5 -6] *
        u'[-1 -2; 11 12]
    )
    temp2 = braid(temp1, (11, 12, 13, 14, 1, 100), (1, 6, 5, 2), (3, 4))
    @planar(
        temp3[-100 -1000 -2000 -200; -300 -400] :=
        w[31 32 41; -300] * w[51 21 22; -400] *
        temp2[42 -1000 -2000 52; 41 51] *
        w'[-100; 31 32 42] * w'[-200; 52 21 22]
    )
    scaled_op = braid(temp3, (11, 100, 1, 12, 13, 14), (1, 4), (5, 6, 3, 2))
    return scaled_op
end

function descend_left(rho::TernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        scaled_rho[-100 -200; -300 -400] :=
        u'[61 62; -400 63] *
        w'[51; 52 -300 61] * w'[21; 62 11 12] *
        rho[42 22; 51 21] *
        w[52 -100 41; 42] * w[31 11 12; 22] *
        u[-200 63; 41 31]
    )
    return scaled_rho
end

function descend_right(rho::TernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        scaled_rho[-100 -200; -300 -400] :=
        u'[62 61; 63 -300] *
        w'[21; 11 12 62] * w'[51; 61 -400 52] *
        rho[22 42; 21 51] *
        w[11 12 41; 22] * w[31 -200 52; 42] *
        u[63 -100; 41 31]
    )
    return scaled_rho
end

function descend_mid(rho::TernaryOperator, layer::TernaryLayer)
    u, w = layer
    # Cost: 6X^6
    @planar(
        scaled_rho[-100 -200; -300 -400] :=
        u'[61 62; -300 -400] *
        w'[21; 11 12 61] * w'[41; 62 31 32] *
        rho[22 42; 21 41] *
        w[11 12 51; 22] * w[52 31 32; 42] *
        u[-100 -200; 51 52]
    )
    return scaled_rho
end

function descend(rho::TernaryOperator, layer::TernaryLayer)
    u, w = layer
    l = descend_left(rho, layer)
    r = descend_right(rho, layer)
    m = descend_mid(rho, layer)
    scaled_rho = (l+r+m)/3
    return scaled_rho
end

# # # Optimization

function environment(op, layer::TernaryLayer, rho; vary_disentanglers = true)
    if vary_disentanglers
        env_u = environment_disentangler(op, layer, rho)
    else
        # The adjoint is just for type stability.
        env_u = zero(layer.disentangler')'
    end
    env_w = environment_isometry(op, layer, rho)
    return TernaryLayer(env_u, env_w)
end

function minimize_expectation_ev(
    layer::TernaryLayer, env::TernaryLayer; vary_disentanglers = true
)
    u = if vary_disentanglers
        projectisometric(env.disentangler; alg = Polar())
    else
        layer.disentangler
    end
    w = projectisometric(env.isometry; alg = Polar())
    return TernaryLayer(u, w)
end

function environment_disentangler(h::TernaryOperator, layer, rho)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        env1[-1 -2; -3 -4] :=
        rho[31 21; 63 22] *
        w[61 51 41; 31] * w[42 11 12; 21] *
        u[52 -2; 41 42] *
        h[62 -1; 51 52] *
        w'[63; 61 62 -3] * w'[22; -4 11 12]
    )

    # Cost: 6X^6
    @planar(
        env2[-1 -2; -3 -4] :=
        rho[41 51; 42 52] *
        w[21 22 61; 41] * w[62 31 32; 51] *
        u[11 12; 61 62] *
        h[-1 -2; 11 12] *
        w'[42; 21 22 -3] * w'[52; -4 31 32]
    )

    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        env3[-1 -2; -3 -4] :=
        rho[21 31; 22 63] *
        w[12 11 42; 21] * w[41 51 61; 31] *
        u[-1 52; 42 41] *
        h[-2 62; 52 51] *
        w'[22; 12 11 -3] * w'[63; -4 62 61]
    )

    env = (env1 + env2 + env3)/3
    return env
end

function environment_isometry(h::TernaryOperator, layer, rho)
    u, w = layer
    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        env1[-1 -2 -3; -4] :=
        rho[81 84; 82 -4] *
        w[62 41 31; 81] * w[83 -2 -3; 84] *
        u[42 52; 31 83] *
        h[61 51; 41 42] *
        u'[63 -1; 51 52] *
        w'[82; 62 61 63]
    )

    # Cost: 6X^6
    @planar(
        env2[-1 -2 -3; -4] :=
        rho[41 62; 42 -4] *
        w[11 12 52; 41] * w[61 -2 -3; 62] *
        u[31 32; 52 61] *
        h[21 22; 31 32] *
        u'[51 -1; 21 22] *
        w'[42; 11 12 51]
    )

    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        env3[-1 -2 -3; -4] :=
        rho[31 33; 32 -4] *
        w[21 11 51; 31] * w[41 61 -3; 33] *
        u[72 62; 51 41] *
        h[71 -2; 62 61] *
        u'[73 -1; 72 71] *
        w'[32; 21 11 73]
    )

    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        env4[-1 -2 -3; -4] :=
        rho[33 31; -4 32] *
        w[-1 61 41; 33] * w[51 11 21; 31] *
        u[62 72; 41 51] *
        h[-2 71; 61 62] *
        u'[-3 73; 71 72] *
        w'[32; 73 11 21]
    )

    # Cost: 6X^6
    @planar(
        env5[-1 -2 -3; -4] :=
        rho[62 41; -4 42] *
        w[-1 -2 61; 62] * w[52 12 11; 41] *
        u[32 31; 61 52] *
        h[22 21; 32 31] *
        u'[-3 51; 22 21] *
        w'[42; 51 12 11]
    )

    # Cost: 2X^8 + 2X^7 + 2X^6
    @planar(
        env6[-1 -2 -3; -4] :=
        rho[84 81; -4 82] *
        w[-1 -2 83; 84] * w[31 41 62; 81] *
        u[52 42; 83 31] *
        h[51 61; 42 41] *
        u'[-3 63; 52 51] *
        w'[82; 63 61 62]
    )

    env = (env1 + env2 + env3 + env4 + env5 + env6)/3
    return env
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 14c4144e-abc0-4802-8b26-f8f6f6a0ddc0
using Omega

# ╔═╡ bc0da6fe-5b09-4d49-9bd8-e0530fac0186
using Distributions

# ╔═╡ ad003e86-23b7-4c8e-a6eb-009b1ec73d25
md"# Design

There are a number of key concepts in Omega
- The Ω object
- Random variable classes
- Random variables
- The conditioning operator
- The intervention operator
- Deep interventions
- The random conditional operator
- The random interventional operator
"

# ╔═╡ 50bfb49f-7081-49e7-bda2-503804599d5e
md""" # Omega
The Ω object in Omega.jl is of special importance.

There are a few perspectives on what Ω is:

1. The sample space in probability theory (often denoted $\Omega$)
2. Space of exogenous variables in as causal graphical modelling nomenclature (often denoted $U$)
3. The random number generator `rng::Random.AbstractRNG`.  Implementation wise, Omega.jl tries to be compatible with any model written as a function `f(rng::AbstractRNG)`
"""

# ╔═╡ 99fc08f8-d2b2-11eb-0599-fba4c4ec5a1a
md"## Conditional Independence"

# ╔═╡ eaa2ba5f-e319-4f74-9adc-87be3bc319d2
md"""
A basic object in Omega is a random variable class.  Each class represents a sequence of random variables.  A class is analogous to a *plate* in Bayesian networks.
"""

# ╔═╡ 69ab85c4-2d16-4e67-88a3-6d50d64e5995
md"""### Primitive Classes
Omega comes with a set of built-in primitive random variable classes, such as `StdNormal` and `StdUniform`.  There are parameterless infinite sets of random variables.
"""

# ╔═╡ e1d86ea5-ace0-456d-8ea6-115ae19a8531
As = StdNormal{Float64}()

# ╔═╡ 9cd46715-b64b-4767-adec-2748560189f7
md"To get a member of a random variable class -- a random variable -- we can use the function `nth`, which gives the `nth` memeber of a class."

# ╔═╡ b625bd45-2289-41e3-89d2-613e0313ede7
A1 = nth(As, 1)

# ╔═╡ 135bde41-b4e8-493f-8e84-f367962aaa6c
A2 = nth(As, 2);

# ╔═╡ c92aaaf2-b1c7-4950-b5cb-836bf7f8fc3a
A3 = nth(As, 3);

# ╔═╡ e8fa68c4-8771-4f4f-92e1-35121e0e726a
md"Or equivalently, use `~`:"

# ╔═╡ bced6a62-c745-4c84-8834-9e2f4d1c62e2
A1_ = 1 ~ As;

# ╔═╡ 0515706c-f4c1-44ac-aab6-58da8c006b1d
A2_ = 2 ~ As;

# ╔═╡ a12ce183-2623-4582-86c1-e1cb206bf78c
A3_ = 3 ~ As;

# ╔═╡ 359a7e70-4662-4942-af96-441e24e15839
md"__Possibly talk about:__ IDs, automatic generation"

# ╔═╡ 0d853129-a457-4650-843c-6443d4797b53
md"#### Sampling
You can sample from a random variable using `randsample`"

# ╔═╡ 4a62d4ff-582e-4dc2-ac44-84c48433e078
randsample(A1)

# ╔═╡ 0d5b28d1-5290-4266-b0ee-4c8924c0ec0e
md"""Random variables are simply functions of ω.  Conceptually, `randsample` then does two distinct steps:
1. Construct a random $ω : Ω$
2. Apply $ω$ to $A$
"""

# ╔═╡ 1408729d-8870-4116-834e-7ebb8fcf0732
ω = defω(); # Step 1 default ω (there are many types)

# ╔═╡ cbcd86a4-cece-4847-8cce-55e768ac63c6
A1(ω) # Step 2

# ╔═╡ 6bd0f9bc-95d5-46c8-be44-edefeb561b0e
md"Note that random variables are pure: reapplication produces the same result:"

# ╔═╡ d92b8bc8-28ca-4944-ba4a-a357f4139014
A1(ω)

# ╔═╡ 0587686f-9087-491c-8301-2f79bbaf2834
md"""__Omega__: In terms of representaiton terms, it is useful to think of Omega is a mapping from a primitive set of random variables to values in their domains."""

# ╔═╡ a5e82bc2-fcce-448c-99ed-a01cda32ac7f
md"### Composite Classes
A class in Omega is actually just function of the form `f(id, ω)`.
Of course, you can specify your own classes simply by constructing a function."

# ╔═╡ 4ccae103-1a6d-42cf-8546-9ccdae847ba5
μ = 0 ~ StdNormal{Float64}()

# ╔═╡ 521de9f7-e941-481e-aa0d-21a263ae19b7
Xs(id, ω) = (id ~ Normal(μ(ω), 1))(ω);

# ╔═╡ fcbe23a1-6180-4d93-bcb7-17f92668cff6
x1 = 1 ~ Xs

# ╔═╡ d6ea7aa6-568f-4f72-a152-83a60ddac0e3
x2 = 2 ~ Xs

# ╔═╡ c8d985a9-9834-427a-be94-a4777c0cc6bc
ω_ = defω()

# ╔═╡ f64a3156-b9bb-4c51-b76b-1effe872b16b
joint_ = @joint x1 x2 μ

# ╔═╡ ee2e69da-58b0-442b-ba37-da42531082fd
joint_(ω_)

# ╔═╡ d8e425ec-c385-4632-9a6f-a392bb84ede0
ω_

# ╔═╡ 3cea826b-e938-4754-b931-cf1645b77c81
randsample(x1)

# ╔═╡ a877b127-d2d8-48d1-8cde-d819adeaa811
md"To construct a random variable that produces the output of many variables, we simply create a new function"

# ╔═╡ 667dd097-cd63-4642-a838-b706e430a1c6
randsample(ω -> (x1 = x1(ω), x2 = x2(ω), μ = μ(ω)))

# ╔═╡ 759ab4f5-5e7a-42be-9de9-f7f4c0f93045
md"Or for convenience we can just use the `@joint` macro, which does the same thing"

# ╔═╡ e2d93606-afcc-4f8e-a636-e7d6f4167250
joint = @joint x1 x2 μ

# ╔═╡ c301d41f-5b20-4017-a3f4-85dce22e0a90
randsample(joint)

# ╔═╡ 6e43e42f-1633-42d1-a695-71851d123705
md"It is important to note that `x1` and `x2` are __conditionally independent__ given μ"

# ╔═╡ 9f7c47a3-bb65-4c27-a120-6cac16e727fe
md"#### Multivariante Distributions

If we want to construct a array-valued random variable from a class we use `Mv`
"

# ╔═╡ 64cdcb20-3480-4c9c-9865-310729df9b9a
X⃗ = Mv(1:10, Xs)

# ╔═╡ e2bf387e-7c8f-4ab7-af70-8d8c887c3574
randsample(X⃗)

# ╔═╡ 44604569-5983-4f77-a884-acef6fad75c6
md"""## (Full) Independence

Sometimes we need to construct random variables that are independent.  The function `iid` constructs a class of random variables that are independent."""

# ╔═╡ 7ef6607a-9048-4b41-ab33-9d942b7c5e19
iidclass = iid(x1)

# ╔═╡ 3987183c-1dd2-4007-966d-979265f8bda5
x1_iid = 1 ~ iidclass

# ╔═╡ a12cc8e9-1209-433c-939d-5dd7027b3f2c
x2_iid = 2 ~ iidclass

# ╔═╡ ef711fcd-4f01-4aeb-94ce-84b374b5aee4
randsample(@joint x1_iid x2_iid)

# ╔═╡ b849e448-54c5-4a93-a77f-714320e42e11
w_2 = defω()

# ╔═╡ e3437550-fad4-4bc5-b1e4-2a30d92f9991
joint2 = @joint  x1_iid x2_iid

# ╔═╡ f1de6343-81c7-4f2f-a6db-c77de43a74f3
joint2(w_2)

# ╔═╡ 51d2ed1f-ee25-429c-a334-2021d1367083
w_2

# ╔═╡ e77a0d45-7773-4df0-881b-3d89c6cda44e
md"It is important to note that `x1_iid` and `x2_iid` are __independent__ (not conditionally independent) and identically distributed"

# ╔═╡ 4de86493-b2bd-4af5-955b-a6f076ba31eb
md"""__Possibly talk about__: How this works.  Sample space projections"""

# ╔═╡ 585d5f13-bd99-45b1-8dcd-1dae0807b0c2
md"### Random Variable Transformations

There is no explicit concept of a probabilistic model in Omega, but we use the term conceptually to refer to the collection of random variables constructed in a program.

Cosntructing a model usually involves defining new primitive random variables of transforming them.
"

# ╔═╡ 1a423af6-d27d-4420-8bf2-f331607c1d20
md"""Transformations of random variables are themselves random variables.  We can do this explicitly by just creating a new function"""

# ╔═╡ 5f89a9a6-2f01-4444-a372-592dfc6eb963
X⃗sum(ω) = sum(X⃗(ω))

# ╔═╡ 22d84114-43f0-4edd-8588-c13cf558bdd2
randsample(@joint X⃗sum X⃗)

# ╔═╡ 8149622b-106a-46ff-8cc3-cea63d40b672
md"""A particularly important kind of random variable are those which output a Boolean.  They're important, as we shall we in a moment, because they are used for conditioning.  They are constructed in the same way as any other"""

# ╔═╡ 959e14b8-9c3f-4e89-b628-47bb27c97275
evidence(ω) = X⃗sum(ω) < -10.0

# ╔═╡ 5e675bd5-9a83-4944-9700-40ae99ef71ad
evidence_2 = X⃗sum <ₚ -10.0

# ╔═╡ 1732577b-e2c6-45a4-8e2c-ecacdcc548cd
randsample(@joint(evidence, evidence_2), 2)

# ╔═╡ f2a5f9ce-ebd1-49a7-a3bf-17ac0b37bfa7
md"#### Conditioning
Conditioning is then an operation which transforms random variables into conditional random variables"

# ╔═╡ 051a770f-4288-4f6a-b25b-b1d4834501c9
μ_posterior = cnd(μ, evidence)

# ╔═╡ 92067906-96bd-4e43-bf37-4f8da067e3f5
randsample(@joint(μ_posterior,evidence), 3; alg = RejectionSample)

# ╔═╡ c35c4b0c-86ec-4c8f-9444-b641e5c604f1
md"__Possibly talk about__: Likelhood based inference and/or predicate relaxation."

# ╔═╡ acebf7b1-d0de-430e-9cf9-fc264897feb0
md"## Interventions

An interventions in Omega.jl is a kind of program.  Specifically, an intervention transforms one random variable into another.
It is done with `intwrvene`
"

# ╔═╡ 3d12697e-7413-4794-92ea-dd4315d773bd
μ_orig = 1 ~ Normal(0, 1)

# ╔═╡ d99811f7-cfa2-42c8-b5c8-1185dd605aff
X = 2 ~ Normal(μ_orig, 1)

# ╔═╡ ca8b38b4-45b5-4177-a54a-0bcc85e50e91
X_intervene = intervene(X, μ_orig => 200)

# ╔═╡ 443a2396-c95f-4f2a-9fb2-18c0776c1715
randsample(@joint μ_orig X X_intervene)

# ╔═╡ 3775ad8a-13e4-4aba-b57d-0b84082ef855
md"We can also intervene a random variable to be another random variable"

# ╔═╡ bd99057b-dc1b-4cd8-a9d3-c2c9d3b4d952
X_intervene_2 = intervene(X, μ_orig => A1);

# ╔═╡ 99728528-28ed-4ca1-bd91-ea56544067f0
randsample(@joint μ_orig X X_intervene X_intervene_2)

# ╔═╡ 84aa55b1-6b62-4eeb-8be2-9d3d141499c6
md"We can also intervene a random variable to be another a function of itself"

# ╔═╡ 10810035-14b9-4fe9-8180-e08e61d9d199
X_intervene_3 = intervene(X, μ_orig => μ_orig +ₚ 10.0);

# ╔═╡ 54b2923e-0bb7-44b4-839d-d03cdb13c16d
md"For syntactic convenience we can use `|ᵈ`"

# ╔═╡ 5df91e58-e78d-4dd5-b183-5596cb9fa674
X |ᵈ (μ_orig => μ_orig +ₚ 10.0);

# ╔═╡ 3814bf03-c6ea-466d-b057-b8ecc9c2e4bc
treatment(ω) = X(ω) - X_intervene_3(ω)

# ╔═╡ 9eefe547-a8b1-446d-9659-5163647cfe0f
randsample(treatment)

# ╔═╡ 32ed4695-ba78-475b-817b-b6a7d749919a
randsample(@joint μ_orig X X_intervene X_intervene_2 X_intervene_3)

# ╔═╡ e0104f49-2f4c-474a-be93-521c55841833
md"""__Possibly talk about__: How interventions work in Omega.jl -- nonstandard execution"""

# ╔═╡ b8f306ef-3fc1-415d-9830-7b71048f41d9
md"#### Counterfactuals

There's no primitive counterfactual operator in Omega.  Rather, they are constructed by composes `intervene` with `cnd`
"

# ╔═╡ f44bdb57-d04c-488e-bd19-db5986ba4ae4
evidence_X = X >ₚ 0.0

# ╔═╡ f074c347-fe70-43ad-82cf-0befc91a0e33
counterfactual = X_intervene_3 |ᶜ evidence_X

# ╔═╡ 35556a02-d304-46a5-aeff-1f28dab98512
(X |ᵈ (μ_orig => μ_orig +ₚ 10.0)) |ᶜ evidence_X

# ╔═╡ 3cb01bb7-adc8-406b-8711-54cfc830f414
(X |ᶜ evidence_X) |ᵈ (μ_orig => μ_orig +ₚ 10.0)

# ╔═╡ fb8cbcdc-eb16-42cd-81bd-fcae78732f2e
X_posterior = (X |ᶜ evidence_X)

# ╔═╡ a395eb05-09ab-484b-83ce-8b62332cc9c0
@joint X X_posterior

# ╔═╡ ce89139c-ae17-4342-9931-85770fe03491
randsample(counterfactual, 100)

# ╔═╡ a2a0bb64-be53-4056-a1af-37e3fbde9392
md"""__Possibly talk about__: The nesting of the intervention and conditoning may seem counter-intuitive or even wrong, but it is not wrong"""

# ╔═╡ 42bf8f3b-c779-4727-b050-151cadee99ff
md"#### Deep Interventions

Sometimes normal interventions don't quite capture what we want.  For example, in causal discovery we might construct a random variable `X` that is a higher-order random variable in the sense that the values it takes are themselves random variables.  Then, rather than say intervene `X` to `x`, we instead want to say intervene the random variable that is outputted by `X` to `x`.  Omega allows this."

# ╔═╡ 7cb712ba-c478-4010-8604-ae30ffd47e4c
μ_hi = 1 ~ Normal(0, 1);

# ╔═╡ 736d14aa-2104-49e4-a966-c1fc288fbd47
σ_hi = 2 ~ Uniform(1, 3);

# ╔═╡ 5388b8a3-58e4-41b4-b6c9-9c36e3eb757e
y_hi = 3 ~ Normal(μ_hi, σ_hi);

# ╔═╡ f2bd636b-6ba8-4807-a08d-12cc3b4cc206
md"`choice` below flips a coin and returns either `μ_hi` or `σ_hi` -- the random variable itself, not a value from the domain"

# ╔═╡ b955bd92-6cb0-41e7-bc92-0de8dbd497d4
choice(ω) = ifelse((4 ~ Bernoulli(0.5))(ω), μ_hi, σ_hi)

# ╔═╡ 77eda7a7-443f-4c0e-a3d9-49e5baf83f79
randsample(choice)

# ╔═╡ 243cfa89-b756-418f-ab2a-8b015927d9c2
int_dist(ω) = ValueIntervention(choice(ω), 5.0)

# ╔═╡ 86170150-13c7-47a6-a258-51e6d4f59f33
randsample(int_dist)

# ╔═╡ cd0fe624-c809-4b5b-8d85-3f236f034a67
joint_hi = @joint(y_hi, μ_hi, σ_hi, choice)

# ╔═╡ 21971f4c-f973-4ef4-baeb-bd85bef8b167
md"Here we construct the higher order intervention."

# ╔═╡ 4e9b0394-5619-4a65-bbc7-427e2208a2fb
joint_hi_ = joint_hi |ᵈ hi(int_dist)

# ╔═╡ ef243319-4446-49e5-9770-aa3a77bc7ca0
randsample(joint_hi_)

# ╔═╡ 516c1e9e-534f-4acb-a6d5-d5a111c50e37
md"### Higher-Order operators

Omega has two so called higher order extensions of the conditioning operator and the intervention operator.  They are called, respectively, the random conditional distribution operator and the random interventional distribution operator.

In short these allow you to easily:
1. Construct conditional expectations, conditional variances, etc as first class random variables.  Collectively we call these distributional properties as they are not functions of individual samples but of the entire distribution.
2. Condition these distributional properties.  This allows one to add statistical knowledge into a model in a principled way.  We call this distributional inference
"

# ╔═╡ 7f868e1b-78f0-477a-b94e-900ecefd4380
θ = 1 ~ Beta(2.0, 2.0)

# ╔═╡ a04a3fd3-c174-4550-930d-57756e797e7e
x = 2 ~ Bernoulli(θ)

# ╔═╡ 3dbe06f6-6e13-4500-9525-12420b1fd236
ridxθ = rid(x, θ)

# ╔═╡ 89e217e7-b367-4f56-9f07-73fdec358fdd
md"Below are samples from the distribution *over expectations*, i.e. the conditional expectation of `x` given `θ`"

# ╔═╡ 5c20809d-7547-471e-8f4d-6337f1ddb71e
samplemean(x; n = 1000) = mean(randsample(x, n))

# ╔═╡ 071f08a1-f4ab-4d11-be9e-0502475f90cd
meandist(ω) = samplemean(ridxθ(ω))

# ╔═╡ 82974f06-19de-4249-a398-53197b04aeb6
randsample(meandist, 100)

# ╔═╡ Cell order:
# ╠═14c4144e-abc0-4802-8b26-f8f6f6a0ddc0
# ╠═bc0da6fe-5b09-4d49-9bd8-e0530fac0186
# ╟─ad003e86-23b7-4c8e-a6eb-009b1ec73d25
# ╟─50bfb49f-7081-49e7-bda2-503804599d5e
# ╟─99fc08f8-d2b2-11eb-0599-fba4c4ec5a1a
# ╟─eaa2ba5f-e319-4f74-9adc-87be3bc319d2
# ╟─69ab85c4-2d16-4e67-88a3-6d50d64e5995
# ╠═e1d86ea5-ace0-456d-8ea6-115ae19a8531
# ╟─9cd46715-b64b-4767-adec-2748560189f7
# ╠═b625bd45-2289-41e3-89d2-613e0313ede7
# ╠═135bde41-b4e8-493f-8e84-f367962aaa6c
# ╠═c92aaaf2-b1c7-4950-b5cb-836bf7f8fc3a
# ╟─e8fa68c4-8771-4f4f-92e1-35121e0e726a
# ╠═bced6a62-c745-4c84-8834-9e2f4d1c62e2
# ╠═0515706c-f4c1-44ac-aab6-58da8c006b1d
# ╠═a12ce183-2623-4582-86c1-e1cb206bf78c
# ╟─359a7e70-4662-4942-af96-441e24e15839
# ╟─0d853129-a457-4650-843c-6443d4797b53
# ╠═4a62d4ff-582e-4dc2-ac44-84c48433e078
# ╟─0d5b28d1-5290-4266-b0ee-4c8924c0ec0e
# ╠═1408729d-8870-4116-834e-7ebb8fcf0732
# ╠═cbcd86a4-cece-4847-8cce-55e768ac63c6
# ╟─6bd0f9bc-95d5-46c8-be44-edefeb561b0e
# ╠═d92b8bc8-28ca-4944-ba4a-a357f4139014
# ╟─0587686f-9087-491c-8301-2f79bbaf2834
# ╟─a5e82bc2-fcce-448c-99ed-a01cda32ac7f
# ╠═4ccae103-1a6d-42cf-8546-9ccdae847ba5
# ╠═521de9f7-e941-481e-aa0d-21a263ae19b7
# ╠═fcbe23a1-6180-4d93-bcb7-17f92668cff6
# ╠═d6ea7aa6-568f-4f72-a152-83a60ddac0e3
# ╠═c8d985a9-9834-427a-be94-a4777c0cc6bc
# ╠═f64a3156-b9bb-4c51-b76b-1effe872b16b
# ╠═ee2e69da-58b0-442b-ba37-da42531082fd
# ╠═d8e425ec-c385-4632-9a6f-a392bb84ede0
# ╠═3cea826b-e938-4754-b931-cf1645b77c81
# ╟─a877b127-d2d8-48d1-8cde-d819adeaa811
# ╠═667dd097-cd63-4642-a838-b706e430a1c6
# ╟─759ab4f5-5e7a-42be-9de9-f7f4c0f93045
# ╠═e2d93606-afcc-4f8e-a636-e7d6f4167250
# ╠═c301d41f-5b20-4017-a3f4-85dce22e0a90
# ╟─6e43e42f-1633-42d1-a695-71851d123705
# ╟─9f7c47a3-bb65-4c27-a120-6cac16e727fe
# ╠═64cdcb20-3480-4c9c-9865-310729df9b9a
# ╠═e2bf387e-7c8f-4ab7-af70-8d8c887c3574
# ╟─44604569-5983-4f77-a884-acef6fad75c6
# ╠═7ef6607a-9048-4b41-ab33-9d942b7c5e19
# ╠═3987183c-1dd2-4007-966d-979265f8bda5
# ╠═a12cc8e9-1209-433c-939d-5dd7027b3f2c
# ╠═ef711fcd-4f01-4aeb-94ce-84b374b5aee4
# ╠═b849e448-54c5-4a93-a77f-714320e42e11
# ╠═e3437550-fad4-4bc5-b1e4-2a30d92f9991
# ╠═f1de6343-81c7-4f2f-a6db-c77de43a74f3
# ╠═51d2ed1f-ee25-429c-a334-2021d1367083
# ╟─e77a0d45-7773-4df0-881b-3d89c6cda44e
# ╟─4de86493-b2bd-4af5-955b-a6f076ba31eb
# ╟─585d5f13-bd99-45b1-8dcd-1dae0807b0c2
# ╟─1a423af6-d27d-4420-8bf2-f331607c1d20
# ╠═5f89a9a6-2f01-4444-a372-592dfc6eb963
# ╠═22d84114-43f0-4edd-8588-c13cf558bdd2
# ╟─8149622b-106a-46ff-8cc3-cea63d40b672
# ╠═959e14b8-9c3f-4e89-b628-47bb27c97275
# ╠═5e675bd5-9a83-4944-9700-40ae99ef71ad
# ╠═1732577b-e2c6-45a4-8e2c-ecacdcc548cd
# ╟─f2a5f9ce-ebd1-49a7-a3bf-17ac0b37bfa7
# ╠═051a770f-4288-4f6a-b25b-b1d4834501c9
# ╠═92067906-96bd-4e43-bf37-4f8da067e3f5
# ╟─c35c4b0c-86ec-4c8f-9444-b641e5c604f1
# ╟─acebf7b1-d0de-430e-9cf9-fc264897feb0
# ╠═3d12697e-7413-4794-92ea-dd4315d773bd
# ╠═d99811f7-cfa2-42c8-b5c8-1185dd605aff
# ╠═ca8b38b4-45b5-4177-a54a-0bcc85e50e91
# ╠═443a2396-c95f-4f2a-9fb2-18c0776c1715
# ╟─3775ad8a-13e4-4aba-b57d-0b84082ef855
# ╠═bd99057b-dc1b-4cd8-a9d3-c2c9d3b4d952
# ╠═99728528-28ed-4ca1-bd91-ea56544067f0
# ╟─84aa55b1-6b62-4eeb-8be2-9d3d141499c6
# ╠═10810035-14b9-4fe9-8180-e08e61d9d199
# ╟─54b2923e-0bb7-44b4-839d-d03cdb13c16d
# ╠═5df91e58-e78d-4dd5-b183-5596cb9fa674
# ╠═3814bf03-c6ea-466d-b057-b8ecc9c2e4bc
# ╠═9eefe547-a8b1-446d-9659-5163647cfe0f
# ╠═32ed4695-ba78-475b-817b-b6a7d749919a
# ╟─e0104f49-2f4c-474a-be93-521c55841833
# ╟─b8f306ef-3fc1-415d-9830-7b71048f41d9
# ╠═f44bdb57-d04c-488e-bd19-db5986ba4ae4
# ╠═f074c347-fe70-43ad-82cf-0befc91a0e33
# ╠═35556a02-d304-46a5-aeff-1f28dab98512
# ╠═3cb01bb7-adc8-406b-8711-54cfc830f414
# ╠═fb8cbcdc-eb16-42cd-81bd-fcae78732f2e
# ╠═a395eb05-09ab-484b-83ce-8b62332cc9c0
# ╠═ce89139c-ae17-4342-9931-85770fe03491
# ╟─a2a0bb64-be53-4056-a1af-37e3fbde9392
# ╟─42bf8f3b-c779-4727-b050-151cadee99ff
# ╠═7cb712ba-c478-4010-8604-ae30ffd47e4c
# ╠═736d14aa-2104-49e4-a966-c1fc288fbd47
# ╠═5388b8a3-58e4-41b4-b6c9-9c36e3eb757e
# ╟─f2bd636b-6ba8-4807-a08d-12cc3b4cc206
# ╠═b955bd92-6cb0-41e7-bc92-0de8dbd497d4
# ╠═77eda7a7-443f-4c0e-a3d9-49e5baf83f79
# ╠═243cfa89-b756-418f-ab2a-8b015927d9c2
# ╠═86170150-13c7-47a6-a258-51e6d4f59f33
# ╠═cd0fe624-c809-4b5b-8d85-3f236f034a67
# ╟─21971f4c-f973-4ef4-baeb-bd85bef8b167
# ╠═4e9b0394-5619-4a65-bbc7-427e2208a2fb
# ╠═ef243319-4446-49e5-9770-aa3a77bc7ca0
# ╟─516c1e9e-534f-4acb-a6d5-d5a111c50e37
# ╠═7f868e1b-78f0-477a-b94e-900ecefd4380
# ╠═a04a3fd3-c174-4550-930d-57756e797e7e
# ╠═3dbe06f6-6e13-4500-9525-12420b1fd236
# ╠═071f08a1-f4ab-4d11-be9e-0502475f90cd
# ╟─89e217e7-b367-4f56-9f07-73fdec358fdd
# ╠═82974f06-19de-4249-a398-53197b04aeb6
# ╠═5c20809d-7547-471e-8f4d-6337f1ddb71e

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
import Base: +, *

# =========================
# Addition operations
# =========================

+(M1::IntervalMatrix, M2::IntervalMatrix) = IntervalMatrix(M1.mat + M2.mat)

+(x::Interval, M::IntervalMatrix) = IntervalMatrix(x .+ M.mat)
+(M::IntervalMatrix, x::Interval) = IntervalMatrix(x .+ M.mat)

+(x::Number, M::IntervalMatrix) = Interval(x) + M
+(M::IntervalMatrix, x::Number) = Interval(x) + M

# =========================
# Multiplication operations
# =========================

*(M1::IntervalMatrix, M2::IntervalMatrix) = IntervalMatrix(M1.mat * M2.mat)

*(x::Interval, M::IntervalMatrix) = IntervalMatrix(x .* M.mat)
*(M::IntervalMatrix, x::Interval) = IntervalMatrix(x .* M.mat)

*(x::Number, M::IntervalMatrix) = Interval(x) * M
*(M::IntervalMatrix, x::Number) = Interval(x) * M