text:
#========================================================================

	Steady state solution for an Aiyagari model with a diffusion process

	Translated julia code based on code from Ben Moll's Website:
		http://www.princeton.edu/%7Emoll/HACTproject.htm

========================================================================#
using LinearAlgebra, SparseArrays, PyPlot

# Model parameters
γ = 2 # CRRA utility
α = 0.35 # parameter for cobb-douglas production function
δ = 0.1 # captal depreciation
z_mean = 1.0 # mean of levels of diffusion process
σ_sq = (.10)^2 # also for diffusion process
corr = exp(-0.3)
ρ = 0.05

K = 3.8 #inital guess of aggregate capital, need this to be close for algorithm to converge
relax = 0.99 # a relaxation parameter
J=40

z_min = 0.5
z_max = 1.5
a_min = -1.0
a_max = 30.0
H =100

# Paramters for the simulation
maxit = 100
maxitK = 100
crit = 10^(-6)
critK = 1e-5
Δ = 1000 #step size for HJB

# Diffusion process in levels
θ = -log(corr)
Var = σ_sq/(2*θ)

# Grid spaces for variables
a = LinRange(a_min, a_max, H)
a = convert(Array,a)
da = (a_max-a_min)/(H-1)

z = LinRange(z_min, z_max, J)
z = convert(Array,z)
dz = (z_max-z_min)/(J-1)
dz_sq = dz^2

aa = a*ones(1,J)
zz = ones(H,1)*z' #check later

# Other variables
μ = θ*(z_mean.-z) #drift for HJB via Ito's lemma
s_sq = σ_sq.*ones(1,J)

# Spaces for the finite difference terms
Vaf, Vab, Vzf, Vzb, Vzz, c = [zeros(H,J) for i in 1:6]

# Construct the matrix to summarize the evolution of z over time
	# this matrix will not change over the finite differences

yy = -(s_sq/dz_sq)' - μ/dz
χ = s_sq/(2*dz_sq)
ζ = μ/dz + (s_sq/(2*dz_sq))'

# Define the Diagonals of this matrix
updiag_z = zeros(H,1)
   for j = 1:J
	   global updiag_z =[updiag_z; repeat([ζ[j]], H, 1)]
   end
updiag_z =(updiag_z[:])

centerdiag_z=repeat([χ[1]+yy[1]],H,1)
   for j = 2:J-1
	   global centerdiag_z = [centerdiag_z; repeat([yy[j]], H, 1)]
   end
centerdiag_z=[centerdiag_z; repeat([yy[J]+ζ[J]], H, 1)]
centerdiag_z = centerdiag_z[:]

lowdiag_z = repeat([χ[2]], H, 1)
   for j=3:J
	   global lowdiag_z = [lowdiag_z; repeat([χ[j]],H,1)]
   end
lowdiag_z=lowdiag_z[:]

Aswitch = sparse(Diagonal(centerdiag_z))+ [zeros(H,H*J);  sparse(Diagonal(lowdiag_z)) zeros(H*(J-1), H)]+ sparse(Diagonal(updiag_z))[(H+1):end,1:(H*J)]

# Inital Guess for the prices
r = α * K^(α-1)-δ
w = (1-α) * K^(α)
v0 = (w*zz + r.* aa).^(1-γ)/(1-γ)/ρ
global v = v0
dist = zeros(1,maxit)

# The main loop, this makes sure markets clear

for iter in 1:maxitK
	println("Main Loop iteration $(iter)")

	for n in 1:maxit
		V=v

		# Forward difference
		Vaf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vaf[H,:] = (w*z .+ r.*a_max).^(-γ)

		# Backward difference
		Vab[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vab[1,:] = (w*z .+ r.*a_min).^(-γ)

		# caculate consumption and savings
		# First the forward difference case
		cf = Vaf.^(-1/γ)
		sf = w*zz + r.*aa -cf

		# consumption and savings for the backward difference
		cb = Vab.^(-1/γ)
		sb = w*zz + r.*aa -cb

		#Consumption and V' at the steady state
		c0 = w*zz + r.*aa
		Va0 = c0.^(-γ)

		# Setup to the Upwind scheme
		If = sf .> 0
		Ib = sb .< 0
		I0 = (1 .- If - Ib)

		#Now implement the upwind scheme
		Va_Upwind = Vaf.*If + Vab.*Ib + Va0.*I0

		c= Va_Upwind.^(-1/γ)
		u = c.^(1-γ)/(1-γ)

		X = -min.(sb,0)/da
		Y = -max.(sf,0)/da + min.(sb,0)/da
		Z = max.(sf,0)/da

		updiag = 0
		   for j = 1:J
			   updiag =[updiag; Z[1:H-1,j];0]
		   end
		updiag =(updiag[:])

		centerdiag=reshape(Y,H*J,1)

		lowdiag = X[2:H,1]
		   for j=2:J
			   lowdiag = [lowdiag; 0; X[2:H,j]]
		   end
		lowdiag=lowdiag[:]

		AA = sparse(Diagonal(centerdiag[:]))+ [zeros(1, H*J); sparse(Diagonal(lowdiag)) zeros(H*J-1,1)] + sparse(Diagonal(updiag))[2:end, 1:(H*J)]

		global A= AA + Aswitch
		B = (1/Δ + ρ)*sparse(I,H*J,H*J)-A

		u_stacked = reshape(u, H*J, 1)
		V_stacked = reshape(V,H*J, 1)

		b = u_stacked + V_stacked/Δ

		V_stacked = B\b

		V = reshape(V_stacked, H,J)
		V_change = V-v
		global v = V

		dist[n] = findmax(V_change)[1]

		if dist[n] < crit
			println("Value Function Converge Iteration = $(n)")
			break
		end
	end
	# Now solve the Fokker-Planck equation
	AT = A'
	b = zeros(H*J,1)

	# Dirty fix to avoid matrix singularities
	i_fix=1
	b[i_fix] =.1
	row = [zeros(1,i_fix-1) 1 zeros(1,H*J-i_fix)]
	AT[i_fix,:] = row

	# Now solve the linear system
	gg = AT\b
	g_sum = gg'*ones(H*J,1)*da*dz
	gg = gg./g_sum

	global g = reshape(gg,H,J)

	# Update aggregate capital
	S = sum(g'*a*da*dz)
	println(S)

	if abs(K-S)<critK
		break
	end

	#update prices
	global K =relax*K + (1-relax)*S
	global r = α * K^(α-1) - δ
	global w = (1-α) * K^(α)
end

# Graphs

ss = w*zz + r.*aa - c
icut = 50
acut=a[1:icut]
sscut=ss[1:icut,:]
gcut = g[1:icut,:]

surface(acut,z,sscut', α=0.7,camera=(45,25),
        colorbar=false, zrotation=90,
        xlabel="Wealth", ylabel="Density")
png("OptimalSavings")

surface(acut,z,gcut', α=0.7,camera=(45,25))
png("KFE_pdf")


text:
begin
    Δt = 5minutes

    field_names = (:b, :u, :v)

    fields_by_case = Dict(
       "free_convection" => (:b),
       "weak_wind_strong_cooling" => (:b, :u, :v),
       "strong_wind_weak_cooling" => (:b, :u, :v),
       "strong_wind" => (:b, :u, :v),
       "strong_wind_no_rotation" => (:b, :u)
    )

    transformation = (b = Transformation(normalization=ZScore()),
                      u = Transformation(normalization=ZScore()),
                      v = Transformation(normalization=ZScore()))

    parameter_set = RiBasedParameterSet

    closure = closure_with_parameters(RiBasedVerticalDiffusivity(Float64;), parameter_set.settings)

    true_parameters = parameter_set.settings

    directory = "calibrate_ri_based_to_6_day_lesbrary"
    isdir(directory) || mkpath(directory)
end

ext: jl


text:
"""
    CountModel
Abstract type which summarizes all models in the CountTimeSeries.jl package.
"""
abstract type CountModel end

"""
    INGARCH
Abstract type which summarizes all models in the INGARCH framework. Subtype of CountModel.
"""
abstract type INGARCH<:CountModel end

"""
    INARMA
Abstract type which summarizes all models in the INARMA framework. Subtype of CountModel.
"""
abstract type INARMA<:CountModel end

"""
    INGARCHModel
Struct to specify an INGARCH model. It contains information on
* `distr`: Distribution
* `link`: Link function
* `pastObs`: Past observations included in the conditional mean definition
* `pastMean`: Past mean included in the conditional mean definition
* `X`: Regressor matrix
* `external`: Indicator which regressors enter the system externally
* `zi`: Indicator. Is zero inflation modelled?
"""
mutable struct INGARCHModel<:INGARCH
    distr::String
    link::String
    pastObs::Array{T, 1} where T<:Integer
    pastMean::Array{T, 1} where T<:Integer
    X::Array{T, 2} where T<: AbstractFloat
    external::Array{Bool, 1}
    zi::Bool
end

"""
    INARCHModel
Struct to specify an INARCH model. It contains information on
* `distr`: Distribution
* `link`: Link function
* `pastObs`: Past observations included in the conditional mean definition
* `X`: Regressor matrix
* `external`: Indicator which regressors enter the system externally
* `zi`: Indicator. Is zero inflation modelled?
"""
mutable struct INARCHModel<:INGARCH
    distr::String
    link::String
    pastObs::Array{T, 1} where T<:Integer
    X::Array{T, 2} where T<: AbstractFloat
    external::Array{Bool, 1}
    zi::Bool
end

"""
    IIDModel
Struct to specify an IID model. It contains information on
* `distr`: Distribution
* `link`: Link function
* `X`: Regressor matrix
* `external`: Indicator which regressors enter the system externally
* `zi`: Indicator. Is zero inflation modelled?
"""
mutable struct IIDModel<:INGARCH
    distr::String
    link::String
    X::Array{T, 2} where T<: AbstractFloat
    external::Array{Bool, 1}
    zi::Bool
end

"""
    INARMAModel
Struct to specify an INARMA model. It contains information on
* `distr`: Distributions (vector)
* `link`: Link functions (vector)
* `pastObs`: Past observations included
* `pastMean`: Past mean included
* `X`: Regressor matrix
* `external`: Indicator which regressors enter the system externally
* `zi`: Indicator. Is zero inflation modelled?
"""
mutable struct INARMAModel<:INARMA
    distr::Array{String, 1}
    link::Array{String, 1}
    pastObs::Array{T, 1} where T<:Integer
    pastMean::Array{T, 1} where T<:Integer
    X::Array{T, 2} where T<: AbstractFloat
    external::Array{Bool, 1}
    zi::Bool
end

"""
    INARModel
Struct to specify an INAR model. It contains information on
* `distr`: Distributions (vector)
* `link`: Link functions (vector)
* `pastObs`: Past observations included
* `X`: Regressor matrix
* `external`: Indicator which regressors enter the system externally
* `zi`: Indicator. Is zero inflation modelled?
"""
mutable struct INARModel<:INARMA
    distr::Array{String, 1}
    link::Array{String, 1}
    pastObs::Array{T, 1} where T<:Integer
    X::Array{T, 2} where T<: AbstractFloat
    external::Array{Bool, 1}
    zi::Bool
end

"""
    INMAModel
Struct to specify an INMA model. It contains information on
* `distr`: Distributions (vector)
* `link`: Link functions (vector)
* `pastMean`: Past mean included
* `X`: Regressor matrix
* `external`: Indicator which regressors enter the system externally
* `zi`: Indicator. Is zero inflation modelled?
"""
mutable struct INMAModel<:INARMA
    distr::Array{String, 1}
    link::Array{String, 1}
    pastMean::Array{T, 1} where T<:Integer
    X::Array{T, 2} where T<: AbstractFloat
    external::Array{Bool, 1}
    zi::Bool
end


text:
#############################
# Example of Horner's Method and
# finding roots of polynomials
# Created for the YouTube channel
# https://www.youtube.com/OscarVeliz
# @author Oscar Veliz
#############################

# Long Evaluation
function longway(a, x)
	result = 0
	for i in 1:length(a)
		result += a[i] * x^(i-1)
	end
	return result
end

# Iterative Horner Evaluation
function horner(a, x)
	result = a[end]
	for i in length(a)-1:-1:1
		result = a[i] + x * result
	end
	return result
end

# Recursive Horner Evaluation
function rHorner(a, x, i)
	if i == length(a)
		return a[i]
	else
		return a[i] + x * rHorner(a,x,i+1)
	end
end

# Synthetic Division
# Given array of coefficients and value x return quotient
function quotient(a, x)
	q = zeros(length(a))
	q[end] = a[end]
	for i in length(a)-1:-1:1
		q[i] = x * q[i+1]+a[i]
	end
	return view(q, 2:length(q))#ignore remainder
end

# Differentiate using Power Rule
function ddx(a)
	d = zeros(length(a)-1)
	for i in 1:length(d)
		d[i] = a[i+1]*i
	end
	return d
end

# Find derivative at a point using Horner
function hornerDeriv(a, x)
	q = quotient(a, x)
	return horner(q,x)
end

# Newton's Method strictly on Polynomials using Horner
function newtonHorner(a, x)
	i = 0
	px = horner(a,x)
	# Power Rule
	#q = ddx(a)
	#qx = horner(q,x)
	# Ruffini's Rule
	qx = hornerDeriv(a,x)
	while abs(px) > 0.0000001 && i < 1000
		if qx == 0 # avoid divide by zero
			return NaN
		end
		x = x - px / qx
		px = horner(a,x)
		#qx = horner(q,x) # Power Rule
		qx = hornerDeriv(a,x) # Ruffini's Rule
		i = i + 1
	end
	if i == 1000 # hit the iteration cap (divierging or stuck in loop)
		return NaN
	end
	r = round(x) # check if rounding will give you a better answer
	if(abs(horner(a,r)) < abs(horner(a,x)))
		x = r
	end
	return x
end

# Assumes that there are no imaginary roots
function findRoots(a)
	roots = zeros(length(a)-1)
	for r in 1:length(roots)
		printPoly(a)
		root = newtonHorner(a,0) # initial guess of zero
		println("r",r," = ",root)
		roots[r] = root
		a = quotient(a,root) # deflation
	end
	return roots
end

# Print in polynomial form instead of array of coefficients
function printPoly(a)
	for i in length(a):-1:1
		# remove ".0" at end
		ai = string(a[i])
		if length(ai) >= 2 && ai[end] == '0' && ai[end-1] == '.'
			ai = ai[1:end-2]
		end
		if i == 1 # no x term
			println(ai)
		elseif i == 2 # no exponent
			print(ai,"x + ")
		else
			print(ai,"x^",(i-1)," + ")
		end
	end
end

p = [-2310, 727, 382, -72, -8, 1]
println(sort(findRoots(p)))

max_stars_repo_name: deboradeben/numerical-veliz

max_issues_repo_name: KulakovaEA/numerical-veliz

max_forks_repo_name: KulakovaEA/numerical-veliz


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays

@testset "conversion" begin
    A = 2 * rand(ComplexF64, (20, 10)) .- 1
    v = rand(ComplexF64, 10)
    w = rand(ComplexF64, 20)
    V = rand(ComplexF64, 10, 3)
    W = rand(ComplexF64, 20, 3)
    α = rand()
    β = rand()
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(M)

    @test Matrix(M) == A
    @test Array(M) == A
    @test convert(AbstractArray, M) == A
    @test convert(AbstractMatrix, M) === A
    @test convert(Matrix, M) === A
    @test convert(Array, M) === A
    @test convert(AbstractMatrix, M) == A
    @test convert(Matrix, M) == A
    @test convert(Array, M) == A
    @test Matrix(M') == A'
    @test Matrix(transpose(M)) == copy(transpose(A))
    @test convert(AbstractMatrix, M') isa Adjoint
    @test convert(Matrix, M*3I) == A*3

    # sparse matrix generation/conversion
    @test sparse(M) == sparse(Array(M))
    @test convert(SparseMatrixCSC, M) == sparse(Array(M))

    B = copy(A)
    B[rand(1:length(A), 30)] .= 0
    MS = LinearMap(B)
    @test sparse(MS) == sparse(Array(MS)) == sparse(B)

    J = LinearMap(I, 10)
    @test J isa LinearMap{Bool}
    @test sparse(J) == Matrix{Bool}(I, 10, 10)
    @test nnz(sparse(J)) == 10
end


text:
# Julia wrapper for header: /usr/include/scip/cons_sos1.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrSOS1(scip)
    ccall((:SCIPincludeConshdlrSOS1, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsSOS1(scip, cons, name, nvars, vars, weights, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsSOS1, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, weights, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicSOS1(scip, cons, name, nvars, vars, weights)
    ccall((:SCIPcreateConsBasicSOS1, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}), scip, cons, name, nvars, vars, weights)
end

function SCIPaddVarSOS1(scip, cons, var, weight)
    ccall((:SCIPaddVarSOS1, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, weight)
end

function SCIPappendVarSOS1(scip, cons, var)
    ccall((:SCIPappendVarSOS1, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}), scip, cons, var)
end

function SCIPgetNVarsSOS1(scip, cons)
    ccall((:SCIPgetNVarsSOS1, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsSOS1(scip, cons)
    ccall((:SCIPgetVarsSOS1, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetWeightsSOS1(scip, cons)
    ccall((:SCIPgetWeightsSOS1, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetConflictgraphSOS1(conshdlr)
    ccall((:SCIPgetConflictgraphSOS1, libscip), Ptr{SCIP_DIGRAPH}, (Ptr{SCIP_CONSHDLR},), conshdlr)
end

function SCIPgetNSOS1Vars(conshdlr)
    ccall((:SCIPgetNSOS1Vars, libscip), Cint, (Ptr{SCIP_CONSHDLR},), conshdlr)
end

function SCIPvarIsSOS1(conshdlr, var)
    ccall((:SCIPvarIsSOS1, libscip), UInt32, (Ptr{SCIP_CONSHDLR}, Ptr{SCIP_VAR}), conshdlr, var)
end

function SCIPvarGetNodeSOS1(conshdlr, var)
    ccall((:SCIPvarGetNodeSOS1, libscip), Cint, (Ptr{SCIP_CONSHDLR}, Ptr{SCIP_VAR}), conshdlr, var)
end

function SCIPnodeGetVarSOS1(conflictgraph, node)
    ccall((:SCIPnodeGetVarSOS1, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_DIGRAPH}, Cint), conflictgraph, node)
end

function SCIPmakeSOS1sFeasible(scip, conshdlr, sol, changed, success)
    ccall((:SCIPmakeSOS1sFeasible, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{SCIP_SOL}, Ptr{UInt32}, Ptr{UInt32}), scip, conshdlr, sol, changed, success)
end


text:
using ASTInterpreter
using Base.Test

# write your own tests here
include("asttest.jl")
include("stepping.jl")
include("evaling.jl")


text:
using TestThings
using Test

@testset "TestThings.jl" begin
    @test foo([1,-1,3,4]) == 4
    @test bar([1,-1,3,4]) == -1
    @test TestThings.greet() == nothing
    @test TestThings.foobar([1]) == (1,1)
end


text:
@nospecialize
function extract_row(rowt::RowT)
    function extract_row_(fields::TypeScope, rowt::RowT)
        @match rowt begin
            RowCons(k, _, _) && if haskey(fields, k) end =>
                throw(RowFieldDuplicate(k))
            RowCons(k, v, rowt) =>
                let fields = ImDict(fields, k => v)
                    extract_row_(fields, rowt)
                end
            RowMono => (fields, nothing)
            RowPoly(Record(rowt)) =>
                extract_row_(fields, rowt)
            RowPoly(t) => (fields, Some(t))
        end
    end
    extract_row_(TypeScope(), rowt)
end

function subst(freshmap::Dict{A, B}, root::HMT) where {A<:HMT, B<:HMT}
    previsit(nothing, root) do ::Nothing, root::HMT
        if root isa A
            nothing, get(freshmap, root, root)
        else
            nothing, root
        end
    end
end

function mk_tcstate(tctx::Vector{HMT}, new_tvar_hook::Union{Nothing, Function}=nothing)
    genvars = Genvar[]
    genvar_links = Set{UInt}[]

    function new_genvar(s::Symbol)::Var
        genlevel = length(genvars) + 1
        genvar = Genvar(genlevel, s)
        push!(genvars, genvar)
        push!(genvar_links, Set{UInt}())
        Var(genvar)
    end

    function unlink(maxlevel :: Integer)
        while true
            level = length(genvars)
            if level <= maxlevel
                break
            end
            pop!(genvars)
            vars = pop!(genvar_links) :: Set{UInt}
            for typevar_id in vars
                tctx[typevar_id] = Var(Refvar(typevar_id))
            end

        end
    end

    function unlink(f::Function, maxlevel :: Integer)
        while true
            level = length(genvars)
            if level <= maxlevel
                break
            end
            gen = pop!(genvars)
            vars = pop!(genvar_links) :: Set{UInt}
            for typevar_id in vars
                v = Var(Refvar(typevar_id))
                tctx[typevar_id] = v
                f(Var(gen), v)
            end
        end
    end

    function fresh_visit(freshmap::TypeScope, a::HMT)
        @match a begin
            Fresh(s) => (freshmap, get(freshmap, s, a))
            Forall(ns, _) => (mk_type_scope(Pair{Symbol, HMT}[n => t for (n, t) in freshmap if !(n in ns)]), a)
            _ => (freshmap, a)
        end
    end
    function fresh(freshmap::TypeScope, a::HMT)
        previsit(fresh_visit, freshmap, a)
    end

    function tvar_of_int(i::Integer)
        Var(Refvar(i))
    end
    function new_tvar()::HMT
        vid = UInt(length(tctx) + 1)
        if new_tvar_hook !== nothing
            new_tvar_hook(vid)
        end
        tvar = tvar_of_int(vid)
        push!(tctx, tvar)
        tvar
    end

    function int_of_tvar(x::Var)::Optional{UInt}
        @match x begin
            Var(Refvar(i)) => Some(i)
            _ => Nothing
        end
    end

    function occur_in(i::Refvar, ty::HMT)
        @switch ty begin
            @case Var(Genvar(_, _))
            return false
            @case Var(i′) && if i′ === i end
            return false
            @case _
            visit_func(x) =
            @match x begin
                Var(i′) && if i′ === i end => false
                _ => true
            end
            return !visit_check(visit_func, ty)
        end
    end

    function prune(x)
        function vfunc(::Nothing, a::HMT)
            nothing, @match a begin
                Var(Refvar(i)) =>
                @match tctx[i] begin
                    Var(Refvar(i′)) &&
                    if i′ === i end => a

                    a =>
                    let t = prune(a)
                        tctx[i] = t
                        t
                    end
                end
                _ => a
            end
        end
        previsit(vfunc, nothing, x)
    end

    function type_less(lhs::HMT, rhs::HMT)
        lhs = prune(lhs)
        rhs = prune(rhs)
        lhs === rhs && return true

        (@match lhs, rhs begin
            (_, Var(_) || Nom(_) || Fresh(_)) => unify(rhs, lhs)

            (Forall(ns1, p1), Forall(ns2, p2)) =>
                (begin
                    pt = Pair{Symbol, HMT}
                    subst1 = mk_type_scope(pt[a => new_genvar(a) for a in ns1])
                    subst2 = mk_type_scope(pt[a => new_tvar() for a in ns2])

                    type_less(fresh(subst1, p1), fresh(subst2, p2))
                end)
            (Var(_), _) => unify(lhs, rhs)
            (_, Forall(ns2, p2)) =>
                (begin
                    pt = Pair{Symbol, HMT}
                    subst_v = mk_type_scope(pt[a => new_tvar() for a in ns2])
                    type_less(lhs, fresh(subst_v, p2))
                end)
            (Nom(_) || Fresh(_), _) => unify(lhs, rhs)

            # A: (forall a. a -> a) -> [int]
            # B: (int -> int) -> [int]
            # A : B
            (Arrow(a1, r1), Arrow(a2, r2)) =>
                type_less(a2, a1) && type_less(r1, r2)

            # A: (forall a. a) int
            # B: list a
            # A : B
            (App(f1, a1), App(f2, a2)) =>
                type_less(f1, f2) && type_less(a1, a2)

            (Tup(xs1), Tup(xs2)) =>
                all(zip(xs1, xs2)) do (lhs, rhs)
                    type_less(lhs, rhs)
                end

            (Record(a), Record(b)) =>
                (begin
                    (m1, ex1) = extract_row(a)
                    (m2, ex2) = extract_row(b)
                    common_keys =
                        intersect(keys(m1), keys(m2))
                    only_by_1 = [k => v for (k, v) in m1 if !(k in common_keys)]
                    only_by_2 = [k => v for (k, v) in m2 if !(k in common_keys)]

                    all(common_keys) do k
                        type_less(m1[k], m2[k])
                    end || return false

                    function row_check(row1, row2, only_by_1, only_by_2)
                        @match row1, row2 begin
                            (nothing, nothing) => isempty(only_by_1) && isempty(only_by_2)
                            (Some(_), nothing) => false
                            (nothing, Some(row2)) => true
                            (Some(row1), Some(row2)) => isempty(only_by_1) && isempty(only_by_2)
                        end
                    end
                    row_check(ex1, ex2, only_by_1, only_by_2)

                end)
            _ => false
        end)
    end

    function unify(lhs::HMT, rhs::HMT)
        lhs = prune(lhs)
        rhs = prune(rhs)
        lhs === rhs && return true

        (@match lhs, rhs begin
            (Forall{N1}(ns1, p1) where N1, Forall{N2}(ns2, p2) where N2) =>
                N1 === N2 &&
                (begin
                    pt = Pair{Symbol, HMT}
                    subst1 = mk_type_scope(pt[a => new_tvar() for a in ns1])
                    genlevel = length(genvar_links)
                    subst2 = mk_type_scope(pt[a => new_genvar(a) for a in ns2])
                    generic = fresh(subst2, p2)
                    unify(fresh(subst1, p1), generic) || return false

                    generic = prune(generic)
                    remap2 = Dict{Var, Fresh}(v => Fresh(k) for (k, v) in subst2)
                    remap1 = Dict{Var, Fresh}()
                    for (k, v) in subst1
                        v = prune(v)
                        haskey(remap2, v) || return false
                        haskey(remap1, v) && return false
                        remap1[v] = Fresh(k)
                    end
                    backmap1 = Dict{Var, Fresh}()
                    backmap2 = Dict{Var, Fresh}()

                    unlink(genlevel) do v::Var, i::Var
                        backmap1[i] = remap1[v]
                        backmap2[i] = remap2[v]
                    end

                    unify(subst(backmap2, p2), p2) &&
                    unify(subst(backmap1, p1), p1)
                end)
            (Var(Refvar(i) && ai), b) =>
            if occur_in(ai, b)
                throw(IllFormedType("a = a -> b"))
            else
                @match b begin
                    Genvar(genlevel, _) =>
                        push!(genvar_links[genlevel], i)
                    _ => nothing
                end
                tctx[i] = b
                true
            end
            (Var(Genvar(_, _)), _) => false
            (a, (Var(_) && b)) => unify(b, a)

            (Arrow(a1, r1), Arrow(a2, r2)) =>
                unify(a1, a2) && unify(r1, r2)

            (App(f1, a1), App(f2, a2)) =>
                unify(f1, f2) && unify(a1, a2)

            (Tup(xs1), Tup(xs2)) =>
                all(zip(xs1, xs2)) do (lhs, rhs)
                    unify(lhs, rhs)
                end

            (Record(a), Record(b)) =>
                (begin
                    (m1, ex1) = extract_row(a)
                    (m2, ex2) = extract_row(b)
                    common_keys =
                        intersect(keys(m1), keys(m2))
                    only_by_1 = [k => v for (k, v) in m1 if !(k in common_keys)]
                    only_by_2 = [k => v for (k, v) in m2 if !(k in common_keys)]

                    all(common_keys) do k
                        unify(m1[k], m2[k])
                    end &&
                    (begin
                        function row_check(row1, row2, only_by_1, only_by_2)
                            @match row1, row2 begin
                                (nothing, nothing) => isempty(only_by_1) && isempty(only_by_2)
                                (Some(_), nothing) => row_check(row2, row1, only_by_2, only_by_1)
                                (nothing, Some(row2)) => isempty(only_by_2) &&
                                    unify(row2, Record(row_of_pairs(only_by_1, RowMono)))
                                (Some(row1), Some(row2)) =>
                                    begin
                                        polyrow = RowPoly(new_tvar())
                                        ex2 = Record(row_of_pairs(only_by_1, polyrow))
                                        ex1 = Record(row_of_pairs(only_by_2, polyrow))
                                        unify(row1, ex1) &&
                                            unify(row2, ex2)
                                    end
                            end
                        end
                        row_check(ex1, ex2, only_by_1, only_by_2)
                    end)
                end)
            _ => false
        end)
    end

    function instantiate(hmt::HMT)
        @match hmt begin
            Forall(ns, t) => begin
                pt = Pair{Symbol, HMT}
                subst = mk_type_scope(pt[a => new_tvar() for a in ns])
                fresh(subst, t)
            end
            _ => hmt
        end
    end

    function generalise(genmap::Dict{UInt, Symbol}, hmt::HMT)
        function visitor(@nospecialize(_), hmt::HMT)
            @match hmt begin
                Var(Refvar(i)) =>
                    let x = get(genmap, i, nothing)
                        x === nothing && return hmt
                        nothing, Fresh(x)
                    end
                _ => (nothing, hmt)
            end
        end

        Forall(Tuple(values(genmap)), previsit(visitor, genmap, hmt))
    end

    (unify = unify,
        type_less = type_less,
        tctx = tctx,
        instantiate = instantiate,
        generalise = generalise,
        new_tvar = new_tvar,
        new_genvar = new_genvar,
        genvar_links = genvar_links,
        genvars = genvars,
        unlink = unlink,
        tvar_of_int = tvar_of_int,
        int_of_tvar = int_of_tvar,
        fresh = fresh,
        occur_in = occur_in,
        prune = prune,
        extract_row = extract_row)
end

function subst_tvar_visit(ctx::Function, root::HMT)
    @match root begin
        Var(Refvar(i)) => (ctx, ctx(i))
        Var(Genvar(_, _)) => error("unexpected genvar outside its scope")
        _ => (ctx, root)
    end
end

(lhs::HMT) ⪯ (rhs::HMT)  = begin
   small_tc = mk_tcstate(HMT[])
   subst_table = Dict{UInt, HMT}()
   function subst(i::UInt)
        get!(subst_table, i) do
            small_tc.new_tvar()
        end
   end
   lhs = previsit(subst_tvar_visit, subst, lhs)
   rhs = previsit(subst_tvar_visit, subst, rhs)
   small_tc.type_less(lhs, rhs)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: resource_groups
using AWS.Compat
using AWS.UUIDs
"""
    CreateGroup()

Creates a resource group with the specified name and description. You can optionally include a resource query, or a service configuration.

# Required Parameters
- `Name`: The name of the group, which is the identifier of the group in other operations. You can't change the name of a resource group after you create it. A resource group name can consist of letters, numbers, hyphens, periods, and underscores. The name cannot start with AWS or aws; these are reserved. A resource group name must be unique within each AWS Region in your AWS account.

# Optional Parameters
- `Configuration`: A configuration associates the resource group with an AWS service and specifies how the service can interact with the resources in the group. A configuration is an array of GroupConfigurationItem elements.  You can specify either a Configuration or a ResourceQuery in a group, but not both.
- `Description`: The description of the resource group. Descriptions can consist of letters, numbers, hyphens, underscores, periods, and spaces.
- `ResourceQuery`: The resource query that determines which AWS resources are members of this group.  You can specify either a ResourceQuery or a Configuration, but not both.
- `Tags`: The tags to add to the group. A tag is key-value pair string.
"""

create_group(Name; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/groups", Dict{String, Any}("Name"=>Name); aws_config=aws_config)
create_group(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/groups", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name), args)); aws_config=aws_config)

"""
    DeleteGroup()

Deletes the specified resource group. Deleting a resource group does not delete any resources that are members of the group; it only deletes the group structure.

# Optional Parameters
- `Group`: The name or the ARN of the resource group to delete.
- `GroupName`: Don't use this parameter. Use Group instead.
"""

delete_group(; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/delete-group"; aws_config=aws_config)
delete_group(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/delete-group", args; aws_config=aws_config)

"""
    GetGroup()

Returns information about a specified resource group.

# Optional Parameters
- `Group`: The name or the ARN of the resource group to retrieve.
- `GroupName`: Don't use this parameter. Use Group instead.
"""

get_group(; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/get-group"; aws_config=aws_config)
get_group(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/get-group", args; aws_config=aws_config)

"""
    GetGroupConfiguration()

Returns the service configuration associated with the specified resource group. AWS Resource Groups supports configurations for the following resource group types:    AWS::EC2::CapacityReservationPool - Amazon EC2 capacity reservation pools. For more information, see Working with capacity reservation groups in the EC2 Users Guide.

# Optional Parameters
- `Group`: The name or the ARN of the resource group.
"""

get_group_configuration(; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/get-group-configuration"; aws_config=aws_config)
get_group_configuration(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/get-group-configuration", args; aws_config=aws_config)

"""
    GetGroupQuery()

Retrieves the resource query associated with the specified resource group.

# Optional Parameters
- `Group`: The name or the ARN of the resource group to query.
- `GroupName`: Don't use this parameter. Use Group instead.
"""

get_group_query(; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/get-group-query"; aws_config=aws_config)
get_group_query(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/get-group-query", args; aws_config=aws_config)

"""
    GetTags()

Returns a list of tags that are associated with a resource group, specified by an ARN.

# Required Parameters
- `Arn`: The ARN of the resource group whose tags you want to retrieve.

"""

get_tags(Arn; aws_config::AWSConfig=global_aws_config()) = resource_groups("GET", "/resources/$(Arn)/tags"; aws_config=aws_config)
get_tags(Arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("GET", "/resources/$(Arn)/tags", args; aws_config=aws_config)

"""
    GroupResources()

Adds the specified resources to the specified group.

# Required Parameters
- `Group`: The name or the ARN of the resource group to add resources to.
- `ResourceArns`: The list of ARNs for resources to be added to the group.

"""

group_resources(Group, ResourceArns; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/group-resources", Dict{String, Any}("Group"=>Group, "ResourceArns"=>ResourceArns); aws_config=aws_config)
group_resources(Group, ResourceArns, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/group-resources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Group"=>Group, "ResourceArns"=>ResourceArns), args)); aws_config=aws_config)

"""
    ListGroupResources()

Returns a list of ARNs of the resources that are members of a specified resource group.

# Optional Parameters
- `Filters`: Filters, formatted as ResourceFilter objects, that you want to apply to a ListGroupResources operation. Filters the results to include only those of the specified resource types.    resource-type - Filter resources by their type. Specify up to five resource types in the format AWS::ServiceCode::ResourceType. For example, AWS::EC2::Instance, or AWS::S3::Bucket.    When you specify a resource-type filter for ListGroupResources, AWS Resource Groups validates your filter resource types against the types that are defined in the query associated with the group. For example, if a group contains only S3 buckets because its query specifies only that resource type, but your resource-type filter includes EC2 instances, AWS Resource Groups does not filter for EC2 instances. In this case, a ListGroupResources request returns a BadRequestException error with a message similar to the following:  The resource types specified as filters in the request are not valid.  The error includes a list of resource types that failed the validation because they are not part of the query associated with the group. This validation doesn't occur when the group query specifies AWS::AllSupported, because a group based on such a query can contain any of the allowed resource types for the query type (tag-based or AWS CloudFormation stack-based queries).
- `Group`: The name or the ARN of the resource group
- `GroupName`: Don't use this parameter. Use Group instead.
- `MaxResults`: The total number of results that you want included on each page of the response. If you do not include this parameter, it defaults to a value that is specific to the operation. If additional items exist beyond the maximum you specify, the NextToken response element is present and has a value (is not null). Include that value as the NextToken request parameter in the next call to the operation to get the next part of the results. Note that the service might return fewer results than the maximum even when there are more results available. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: The parameter for receiving additional results if you receive a NextToken response in a previous request. A NextToken response indicates that more output is available. Set this parameter to the value provided by a previous call's NextToken response to indicate where the output should continue from.
"""

list_group_resources(; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/list-group-resources"; aws_config=aws_config)
list_group_resources(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/list-group-resources", args; aws_config=aws_config)

"""
    ListGroups()

Returns a list of existing resource groups in your account.

# Optional Parameters
- `Filters`: Filters, formatted as GroupFilter objects, that you want to apply to a ListGroups operation.    resource-type - Filter the results to include only those of the specified resource types. Specify up to five resource types in the format AWS::ServiceCode::ResourceType . For example, AWS::EC2::Instance, or AWS::S3::Bucket.    configuration-type - Filter the results to include only those groups that have the specified configuration types attached. The current supported values are:   AWS:EC2::CapacityReservationPool
- `maxResults`: The total number of results that you want included on each page of the response. If you do not include this parameter, it defaults to a value that is specific to the operation. If additional items exist beyond the maximum you specify, the NextToken response element is present and has a value (is not null). Include that value as the NextToken request parameter in the next call to the operation to get the next part of the results. Note that the service might return fewer results than the maximum even when there are more results available. You should check NextToken after every operation to ensure that you receive all of the results.
- `nextToken`: The parameter for receiving additional results if you receive a NextToken response in a previous request. A NextToken response indicates that more output is available. Set this parameter to the value provided by a previous call's NextToken response to indicate where the output should continue from.
"""

list_groups(; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/groups-list"; aws_config=aws_config)
list_groups(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/groups-list", args; aws_config=aws_config)

"""
    SearchResources()

Returns a list of AWS resource identifiers that matches the specified query. The query uses the same format as a resource query in a CreateGroup or UpdateGroupQuery operation.

# Required Parameters
- `ResourceQuery`: The search query, using the same formats that are supported for resource group definition. For more information, see CreateGroup.

# Optional Parameters
- `MaxResults`: The total number of results that you want included on each page of the response. If you do not include this parameter, it defaults to a value that is specific to the operation. If additional items exist beyond the maximum you specify, the NextToken response element is present and has a value (is not null). Include that value as the NextToken request parameter in the next call to the operation to get the next part of the results. Note that the service might return fewer results than the maximum even when there are more results available. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: The parameter for receiving additional results if you receive a NextToken response in a previous request. A NextToken response indicates that more output is available. Set this parameter to the value provided by a previous call's NextToken response to indicate where the output should continue from.
"""

search_resources(ResourceQuery; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/resources/search", Dict{String, Any}("ResourceQuery"=>ResourceQuery); aws_config=aws_config)
search_resources(ResourceQuery, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/resources/search", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceQuery"=>ResourceQuery), args)); aws_config=aws_config)

"""
    Tag()

Adds tags to a resource group with the specified ARN. Existing tags on a resource group are not changed if they are not specified in the request parameters.  Do not store personally identifiable information (PII) or other confidential or sensitive information in tags. We use tags to provide you with billing and administration services. Tags are not intended to be used for private or sensitive data.

# Required Parameters
- `Arn`: The ARN of the resource group to which to add tags.
- `Tags`: The tags to add to the specified resource group. A tag is a string-to-string map of key-value pairs.

"""

tag(Arn, Tags; aws_config::AWSConfig=global_aws_config()) = resource_groups("PUT", "/resources/$(Arn)/tags", Dict{String, Any}("Tags"=>Tags); aws_config=aws_config)
tag(Arn, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("PUT", "/resources/$(Arn)/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Tags"=>Tags), args)); aws_config=aws_config)

"""
    UngroupResources()

Removes the specified resources from the specified group.

# Required Parameters
- `Group`: The name or the ARN of the resource group from which to remove the resources.
- `ResourceArns`: The ARNs of the resources to be removed from the group.

"""

ungroup_resources(Group, ResourceArns; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/ungroup-resources", Dict{String, Any}("Group"=>Group, "ResourceArns"=>ResourceArns); aws_config=aws_config)
ungroup_resources(Group, ResourceArns, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/ungroup-resources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Group"=>Group, "ResourceArns"=>ResourceArns), args)); aws_config=aws_config)

"""
    Untag()

Deletes tags from a specified resource group.

# Required Parameters
- `Arn`: The ARN of the resource group from which to remove tags. The command removed both the specified keys and any values associated with those keys.
- `Keys`: The keys of the tags to be removed.

"""

untag(Arn, Keys; aws_config::AWSConfig=global_aws_config()) = resource_groups("PATCH", "/resources/$(Arn)/tags", Dict{String, Any}("Keys"=>Keys); aws_config=aws_config)
untag(Arn, Keys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("PATCH", "/resources/$(Arn)/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Keys"=>Keys), args)); aws_config=aws_config)

"""
    UpdateGroup()

Updates the description for an existing group. You cannot update the name of a resource group.

# Optional Parameters
- `Description`: The new description that you want to update the resource group with. Descriptions can contain letters, numbers, hyphens, underscores, periods, and spaces.
- `Group`: The name or the ARN of the resource group to modify.
- `GroupName`: Don't use this parameter. Use Group instead.
"""

update_group(; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/update-group"; aws_config=aws_config)
update_group(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/update-group", args; aws_config=aws_config)

"""
    UpdateGroupQuery()

Updates the resource query of a group.

# Required Parameters
- `ResourceQuery`: The resource query to determine which AWS resources are members of this resource group.

# Optional Parameters
- `Group`: The name or the ARN of the resource group to query.
- `GroupName`: Don't use this parameter. Use Group instead.
"""

update_group_query(ResourceQuery; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/update-group-query", Dict{String, Any}("ResourceQuery"=>ResourceQuery); aws_config=aws_config)
update_group_query(ResourceQuery, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups("POST", "/update-group-query", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceQuery"=>ResourceQuery), args)); aws_config=aws_config)


text:
#####
##### constructors
#####

ChainRules.@non_differentiable (::Type{T} where {T<:Array})(::UndefInitializer, args...)

function rrule(::Type{T}, x::AbstractArray) where {T<:Array}
    project_x = ProjectTo(x)
    Array_pullback(ȳ) = (NoTangent(), project_x(ȳ))
    return T(x), Array_pullback
end

#####
##### `vect`
#####

@non_differentiable Base.vect()

# Case of uniform type `T`: the data passes straight through,
# so no projection should be required.
function rrule(::typeof(Base.vect), X::Vararg{T, N}) where {T, N}
    vect_pullback(ȳ) = (NoTangent(), NTuple{N}(ȳ)...)
    return Base.vect(X...), vect_pullback
end

# Numbers and arrays are often promoted, to make a uniform vector.
# ProjectTo here reverses this
function rrule(
    ::typeof(Base.vect),
    X::Vararg{Union{Number,AbstractArray{<:Number}}, N},
) where {N}
    projects = map(ProjectTo, X)
    function vect_pullback(ȳ)
        X̄ = ntuple(n -> projects[n](ȳ[n]), N)
        return (NoTangent(), X̄...)
    end
    return Base.vect(X...), vect_pullback
end

# Data is unmodified, so no need to project.
function rrule(::typeof(Base.vect), X::Vararg{Any,N}) where {N}
    vect_pullback(ȳ) = (NoTangent(), ntuple(n -> ȳ[n], N)...)
    return Base.vect(X...), vect_pullback
end

#####
##### `reshape`
#####

function rrule(::typeof(reshape), A::AbstractArray, dims::Tuple{Vararg{Union{Colon,Int}}})
    A_dims = size(A)
    function reshape_pullback(Ȳ)
        return (NoTangent(), reshape(Ȳ, A_dims), NoTangent())
    end
    return reshape(A, dims), reshape_pullback
end

function rrule(::typeof(reshape), A::AbstractArray, dims::Union{Colon,Int}...)
    A_dims = size(A)
    function reshape_pullback(Ȳ)
        ∂A = reshape(Ȳ, A_dims)
        ∂dims = broadcast(_ -> NoTangent(), dims)
        return (NoTangent(), ∂A, ∂dims...)
    end
    return reshape(A, dims...), reshape_pullback
end

#####
##### `repeat`
#####
function rrule(::typeof(repeat), xs::AbstractArray; inner=ntuple(_->1, ndims(xs)), outer=ntuple(_->1, ndims(xs)))

    project_Xs = ProjectTo(xs)
    S = size(xs)
    function repeat_pullback(ȳ)
        dY = unthunk(ȳ)
        Δ′ = zero(xs)
        # Loop through each element of Δ, calculate source dimensions, accumulate into Δ′
        for (dest_idx, val) in pairs(IndexCartesian(), dY)
            # First, round dest_idx[dim] to nearest gridpoint defined by inner[dim], then
            # wrap around based on original size S.
            src_idx = [mod1(div(dest_idx[dim] - 1, inner[dim]) + 1, S[dim]) for dim in 1:length(S)]
            Δ′[src_idx...] += val
        end
        x̄ = project_Xs(Δ′)
        return (NoTangent(), x̄)
    end

    return repeat(xs; inner = inner, outer = outer), repeat_pullback
end

function rrule(::typeof(repeat), xs::AbstractArray, counts::Integer...)

    project_Xs = ProjectTo(xs)
    S = size(xs)
    function repeat_pullback(ȳ)
        dY = unthunk(ȳ)
        size2ndims = ntuple(d -> isodd(d) ? get(S, 1+d÷2, 1) : get(counts, d÷2, 1), 2*ndims(dY))
        reduced = sum(reshape(dY, size2ndims); dims = ntuple(d -> 2d, ndims(dY)))
        x̄ = project_Xs(reshape(reduced, S))
        return (NoTangent(), x̄, map(_->NoTangent(), counts)...)
    end
    return repeat(xs, counts...), repeat_pullback
end

#####
##### `hcat`
#####

function rrule(::typeof(hcat), Xs::Union{AbstractArray, Number}...)
    Y = hcat(Xs...)  # note that Y always has 1-based indexing, even if X isa OffsetArray
    ndimsY = Val(ndims(Y))  # this avoids closing over Y, Val() is essential for type-stability
    sizes = map(size, Xs)   # this avoids closing over Xs
    project_Xs = map(ProjectTo, Xs)
    function hcat_pullback(ȳ)
        dY = unthunk(ȳ)
        hi = Ref(0)  # Ref avoids hi::Core.Box
        dXs = map(project_Xs, sizes) do project, sizeX
            ndimsX = length(sizeX)
            lo = hi[] + 1
            hi[] += get(sizeX, 2, 1)
            ind = ntuple(ndimsY) do d
                if d==2
                    d > ndimsX ? lo : lo:hi[]
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            dX = if ndimsX > 0
                # Here InplaceableThunk breaks @inferred, removed for now
                # InplaceableThunk(dX -> dX .+= view(dY, ind...), @thunk(dY[ind...]))
                dY[ind...]
            else
                # This is a hack to perhaps avoid GPU scalar indexing
                sum(view(dY, ind...))
            end
            return project(dX)
        end
        return (NoTangent(), dXs...)
    end
    return Y, hcat_pullback
end

function rrule(::typeof(reduce), ::typeof(hcat), As::AbstractVector{<:AbstractVecOrMat})
    widths = map(A -> size(A,2), As)
    function reduce_hcat_pullback_2(dY)
        hi = Ref(0)
        dAs = map(widths) do w
            lo = hi[]+1
            hi[] += w
            dY[:, lo:hi[]]
        end
        return (NoTangent(), NoTangent(), dAs)
    end
    return reduce(hcat, As), reduce_hcat_pullback_2
end

function rrule(::typeof(reduce), ::typeof(hcat), As::AbstractVector{<:AbstractVector})
    axe = axes(As,1)
    function reduce_hcat_pullback_1(dY)
        hi = Ref(0)
        dAs = map(_ -> dY[:, hi[]+=1], axe)
        return (NoTangent(), NoTangent(), dAs)
    end
    return reduce(hcat, As), reduce_hcat_pullback_1
end

#####
##### `vcat`
#####

function rrule(::typeof(vcat), Xs::Union{AbstractArray, Number}...)
    Y = vcat(Xs...)
    ndimsY = Val(ndims(Y))
    sizes = map(size, Xs)
    project_Xs = map(ProjectTo, Xs)
    function vcat_pullback(ȳ)
        dY = unthunk(ȳ)
        hi = Ref(0)
        dXs = map(project_Xs, sizes) do project, sizeX
            ndimsX = length(sizeX)
            lo = hi[] + 1
            hi[] += get(sizeX, 1, 1)
            ind = ntuple(ndimsY) do d
                if d==1
                    d > ndimsX ? lo : lo:hi[]
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            dX = if ndimsX > 0
                # InplaceableThunk(@thunk(dY[ind...]), dX -> dX .+= view(dY, ind...))
                dY[ind...]
            else
                sum(view(dY, ind...))
            end
            return project(dX)
        end
        return (NoTangent(), dXs...)
    end
    return Y, vcat_pullback
end

function rrule(::typeof(reduce), ::typeof(vcat), As::AbstractVector{<:AbstractVecOrMat})
    Y = reduce(vcat, As)
    ndimsY = Val(ndims(Y))
    heights = map(A -> size(A,1), As)
    function reduce_vcat_pullback(dY)
        hi = Ref(0)
        dAs = map(heights) do z
            lo = hi[]+1
            hi[] += z
            ind = ntuple(d -> d==1 ? (lo:hi[]) : (:), ndimsY)
            dY[ind...]
        end
        return (NoTangent(), NoTangent(), dAs)
    end
    return Y, reduce_vcat_pullback
end

#####
##### `cat`
#####

_val(::Val{x}) where {x} = x

function rrule(::typeof(cat), Xs::Union{AbstractArray, Number}...; dims)
    Y = cat(Xs...; dims=dims)
    cdims = dims isa Val ? Int(_val(dims)) : dims isa Integer ? Int(dims) : Tuple(dims)
    ndimsY = Val(ndims(Y))
    sizes = map(size, Xs)
    project_Xs = map(ProjectTo, Xs)
    function cat_pullback(ȳ)
        dY = unthunk(ȳ)
        prev = fill(0, _val(ndimsY))  # note that Y always has 1-based indexing, even if X isa OffsetArray
        dXs = map(project_Xs, sizes) do project, sizeX
            ndimsX = length(sizeX)
            index = ntuple(ndimsY) do d
                if d in cdims
                    d > ndimsX ? (prev[d]+1) : (prev[d]+1:prev[d]+sizeX[d])
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            for d in cdims
                prev[d] += get(sizeX, d, 1)
            end
            dX = if ndimsX > 0
                # InplaceableThunk(@thunk(dY[index...]), dX -> dX .+= view(dY, index...))
                dY[index...]
            else
                sum(view(dY, index...))
            end
            return project(dX)
        end
        return (NoTangent(), dXs...)
    end
    return Y, cat_pullback
end

#####
##### `hvcat`
#####

function rrule(::typeof(hvcat), rows, values::Union{AbstractArray, Number}...)
    Y = hvcat(rows, values...)
    cols = size(Y,2)
    ndimsY = Val(ndims(Y))
    sizes = map(size, values)
    project_Vs = map(ProjectTo, values)
    function hvcat_pullback(dY)
        prev = fill(0, 2)
        dXs = map(project_Vs, sizes) do project, sizeX
            ndimsX = length(sizeX)
            index = ntuple(ndimsY) do d
                if d in (1, 2)
                    d > ndimsX ? (prev[d]+1) : (prev[d]+1:prev[d]+sizeX[d])
                else
                    d > ndimsX ? 1 : (:)
                end
            end
            prev[2] += get(sizeX, 2, 1)
            if prev[2] == cols
                prev[2] = 0
                prev[1] += get(sizeX, 1, 1)
            end
            project(dY[index...])
        end
        return (NoTangent(), NoTangent(), dXs...)
    end
    return Y, hvcat_pullback
end

#####
##### `reverse`
#####

# 1-dim case allows start/stop, N-dim case takes dims keyword
# whose defaults changed in Julia 1.6... just pass them all through:

function frule((_, xdot), ::typeof(reverse), x::AbstractArray, args...; kw...)
    return reverse(x, args...; kw...), reverse(xdot, args...; kw...)
end

function rrule(::typeof(reverse), x::AbstractArray, args...; kw...)
    project = ProjectTo(x)
    nots = map(_ -> NoTangent(), args)
    function reverse_pullback(dy)
        dx = @thunk project(reverse(unthunk(dy), args...; kw...))
        return (NoTangent(), dx, nots...)
    end
    return reverse(x, args...; kw...), reverse_pullback
end

#####
##### `circshift`
#####

function frule((_, xdot), ::typeof(circshift), x::AbstractArray, shifts)
    return circshift(x, shifts), circshift(xdot, shifts)
end

function rrule(::typeof(circshift), x::AbstractArray, shifts)
    project = ProjectTo(x)
    function circshift_pullback(dy)
        dx = @thunk project(circshift(unthunk(dy), map(-, shifts)))
        # Note that circshift! is useless for InplaceableThunk, as it overwrites completely
        return (NoTangent(), dx, NoTangent())
    end
    return circshift(x, shifts), circshift_pullback
end

#####
##### `fill`
#####

function frule((_, xdot), ::typeof(fill), x::Any, dims...)
    return fill(x, dims...), fill(xdot, dims...)
end

function rrule(::typeof(fill), x::Any, dims...)
    project = x isa Union{Number, AbstractArray{<:Number}} ? ProjectTo(x) : identity
    nots = map(_ -> NoTangent(), dims)
    fill_pullback(Ȳ) = (NoTangent(), project(sum(Ȳ)), nots...)
    return fill(x, dims...), fill_pullback
end


text:
@testset "Quadratic" begin
    for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy))
        isinplace = constructor == interpolate!
        f(x) = sin((x-3)*2pi/9 - 1)
        xmax = 10
        A = Float64[f(x) for x in 1:xmax]
        for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
            itp1 = @inferred(constructor(copier(A), BSpline(Quadratic(BC(GT())))))
            ax1 = axes(itp1)[1]
            @test Interpolations.lbounds(itp1) == (GT == OnGrid ? (first(ax1),) : (first(ax1) - 0.5,))
            @test Interpolations.ubounds(itp1) == (GT == OnGrid ? (last(ax1),) : (last(ax1) + 0.5,))
            @test_throws ArgumentError parent(itp1)
            check_axes(itp1, A, isinplace)
            check_inbounds_values(itp1, A)
            check_oob(itp1)
            can_eval_near_boundaries(itp1)

            # test that inner region is close to data
            for x in 3.1:.2:8.1
                @test f(x) ≈ itp1(x) atol=abs(0.1 * f(x))
            end
        end

        f(x,y) = sin(x/10)*cos(y/6)
        xmax, ymax = 30,10
        A = Float64[f(x,y) for x in 1:xmax, y in 1:ymax]

        # test that inner region is close to data
        for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
            itp2 = @inferred(constructor(copier(A), BSpline(Quadratic(BC(GT())))))
            check_axes(itp2, A, isinplace)
            check_inbounds_values(itp2, A)
            check_oob(itp2)
            can_eval_near_boundaries(itp2)

            for x in 3.1:.2:xmax-3, y in 3.1:2:ymax-3
                @test f(x,y) ≈ itp2(x,y) atol=abs(0.1 * f(x,y))
            end
        end
    end

    # InPlace
    let
        f(x) = sin((x-3)*2pi/9 - 1)
        xmax = 10
        A = Float64[f(x) for x in 1:xmax]
        itp1 = interpolate!(copy(A), BSpline(Quadratic(InPlace(OnCell()))))
        @test axes(itp1) == axes(A)
        check_inbounds_values(itp1, A)

        f(x,y) = sin(x/10)*cos(y/6)
        xmax, ymax = 30,10
        A = Float64[f(x,y) for x in 1:xmax, y in 1:ymax]
        itp2 = interpolate!(copy(A), BSpline(Quadratic(InPlace(OnCell()))))
        @test axes(itp2) == axes(A)
        check_inbounds_values(itp2, A)
    end
end


text:
using Plots
using DifferentialEquations

labs = ["Prey", "Predator"];
function f(du,u,p,t)
    du[1] = (p[1]*u[1]) - (p[2]*u[1]*u[2])
    du[2] = (p[2]*u[1]*u[2]) - (p[3]*u[2])
end
u0 = [10.0;4.0]
p = [0.1;0.01;0.01]
tspan = (0.0, 1000.0)
prob = ODEProblem(f,u0,tspan,p)

# You need to solve the problem.  i.e. you need to run the simulation:
res = solve(prob, reltol=1e-14,abstol=1e-14)

# Then you need to get the simulated values in a nice format for plotting
t = res.t
prey = [values[1] for values in res.u]
pred = [values[2] for values in res.u]

# Plot the dynamics as usual
p1 = plot(t,prey, label="Prey")
plot!(t,pred, title="Predator Prey Dynamics", yaxis="Population Size", xaxis="Time (Years)", label="Predator",legend=:topright)

# Plot the phase-space diagram
p2 = plot(pred,prey, title="Phase Space Diagram", yaxis="Prey Population", xaxis="Predator Population", legend=false)

# Combine the two plots side by side
plot(p1,p2,layout=(1,2),linewidth=3 )


text:
using EngEconomics

costA = 15000
omA = 55
benefitA = 450
costB = 10000
omB = 35
benefitB = 200
n = 30
i = 0.05

# Consider do nothing option and pick best from b/c analysis
# 1.) Rank options by smallest initial cost to biggest initial cost
#     (1) B
#     (2) A
# 2.) Compare 1 and 2
ΔEUAB_21 = (benefitB - benefitA) - (omB - omA)
ΔEUAC_21 = (costB - costA) * capitalRecoveryFactor(i, n)

ΔBCRatio_21 = ΔEUAB_21 / ΔEUAC_21

if ΔBCRatio_21 >= 1
	println("Pick the one with the bigger initial cost, i.e. B")
else
	println("Pick the one with the smaller initial cost, i.e. Do Nothing")
end

disbenefitA = 400
disbenefitB = 500

ΔEUAB_21 = (benefitB - benefitA) - (omB - omA) - (disbenefitB - disbenefitA)
ΔEUAC_21 = (costB - costA) * capitalRecoveryFactor(i, n)

ΔBCRatio_21 = ΔEUAB_21 / ΔEUAC_21

if ΔBCRatio_21 >= 1
	println("Pick the one with the bigger initial cost, i.e. B")
else
	println("Pick the one with the smaller initial cost, i.e. Do Nothing")
end


text:
module OpenEEW

using Dates
using AWSCore, AWSS3, AWSSDK, DataFrames, Interpolations, JSON2, SeisIO, StructArrays
include("io.jl")
include("aws.jl")
include("processing.jl")
end


text:
using Anderson
using Test

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


text:
mutable struct DistMatrix{T} <: ElementalMatrix{T}
	obj::Ptr{Cvoid}
    g::Grid # keep the grid around to avoid that it's freed before the matrix
end

for (elty, ext) in ((:ElInt, :i),
                    (:Float32, :s),
                    (:Float64, :d),
                    (:ComplexF32, :c),
                    (:ComplexF64, :z))
    @eval begin

        # destructor to be used in finalizer. Don't call explicitly
        function destroy(A::DistMatrix{$elty})
            ElError(ccall(($(string("ElDistMatrixDestroy_", ext)), libEl), Cuint,
                (Ptr{Cvoid},), A.obj))
            return nothing
        end

        function DistMatrix(::Type{$elty}, colDist::Dist = MC, rowDist::Dist = MR, grid::Grid = DefaultGrid[])
            obj = Ref{Ptr{Cvoid}}(C_NULL)
            ElError(ccall(($(string("ElDistMatrixCreateSpecific_", ext)), libEl), Cuint,
                (Cint, Cint, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                colDist, rowDist, grid.obj, obj))
            A = DistMatrix{$elty}(obj[], grid)
            finalizer(destroy, A)
            return A
        end

        # Probably not necesary to have this function as we carry around a reference to the Grid.
        # function Grid(A::DistMatrix{$elty})
        #     g = Grid()
        #     ElError(ccall(($(string("ElDistMatrixGrid_", ext)), libEl), Cuint,
        #         (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
        #         A.obj, Ref{Ptr{Cvoid}}(g.obj)))
        #     return g
        # end

        function comm(A::DistMatrix{$elty})
            cm = Ref{ElComm}()
            ElError(ccall(($(string("ElDistMatrixDistComm_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElComm}),
                A.obj, cm))
            return cm[]
        end

        function get(A::DistMatrix{$elty}, i::Integer, j::Integer)
            rv = Ref{$elty}(0)
            ElError(ccall(($(string("ElDistMatrixGet_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt, Ref{$elty}),
                A.obj, i - 1, j - 1, rv))
            return rv[]
        end

        function getLocal(A::DistMatrix{$elty}, i::Integer, j::Integer)
            rv = Ref{$elty}(0)
            ElError(ccall(($(string("ElDistMatrixGetLocal_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt, Ref{$elty}),
                A.obj, i - 1, j - 1, rv))
            return rv[]
        end

        function globalCol(A::DistMatrix{$elty}, i::Integer)
            rv = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistMatrixGlobalCol_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, Ref{ElInt}),
                A.obj, i - 1, rv))
            return rv[] + 1
        end

        function globalRow(A::DistMatrix{$elty}, i::Integer)
            rv = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistMatrixGlobalRow_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, Ref{ElInt}),
                A.obj, i - 1, rv))
            return rv[] + 1
        end

        function height(A::DistMatrix{$elty})
            rv = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistMatrixHeight_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, rv))
            return rv[]
        end

        function localHeight(A::DistMatrix{$elty})
            rv = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistMatrixLocalHeight_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, rv))
            return rv[]
        end

        function localWidth(A::DistMatrix{$elty})
            rv = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistMatrixLocalWidth_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, rv))
            return rv[]
        end

        function processPullQueue(A::DistMatrix{$elty}, buf::Array{$elty,2})
            ElError(ccall(($(string("ElDistMatrixProcessPullQueue_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ptr{$elty}),
                A.obj, buf))
            return buf
        end

        function processQueues(A::DistMatrix{$elty})
            ElError(ccall(($(string("ElDistMatrixProcessQueues_", ext)), libEl), Cuint,
                (Ptr{Cvoid},), A.obj))
            return A
        end

        function queuePull(A::DistMatrix{$elty}, i::Integer, j::Integer)
            ElError(ccall(($(string("ElDistMatrixQueuePull_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt),
                A.obj, i - 1, j - 1))
            return nothing
        end

        function queueUpdate(A::DistMatrix{$elty}, i::Integer, j::Integer, value::$elty)
            ElError(ccall(($(string("ElDistMatrixQueueUpdate_", ext)), libEl), Cuint,
              (Ptr{Cvoid}, ElInt, ElInt, $elty),
              A.obj, i - 1, j - 1, value))
            return nothing
        end

        function reserve(A::DistMatrix{$elty}, numEntries::Integer)
            ElError(ccall(($(string("ElDistMatrixReserve_", ext)), libEl), Cuint,
              (Ptr{Cvoid}, ElInt),
              A.obj, numEntries))
            return nothing
        end

        function width(A::DistMatrix{$elty})
            rv = Ref{ElInt}(0)
            ElError(ccall(($(string("ElDistMatrixWidth_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, Ref{ElInt}),
                A.obj, rv))
            return rv[]
        end

        function resize!(A::DistMatrix{$elty}, i::Integer, j::Integer = 1) # to mimic vector behavior
            ElError(ccall(($(string("ElDistMatrixResize_", ext)), libEl), Cuint,
                (Ptr{Cvoid}, ElInt, ElInt),
                A.obj, i, j))
            return A
        end
    end
end

DistMatrix() = DistMatrix(Float64)

#########################
### Julia convenience ###
#########################

# Do I want to provide this function? It's an invitation to be slow
getindex(A::DistMatrix, i::Integer, j::Integer) = get(A, i, j)

# This might be wrong. Should consider how to extract distributions properties of A
function similar(::DistMatrix, ::Type{T}, sz::Dims) where {T}
    A = DistMatrix(T)
    resize!(A, sz...)
    return A
end

function getindex(A::DistMatrix, iInd::Colon, jInd::UnitRange)
    B = DistMatrix(eltype(A))
    zeros!(B, size(A, 1), length(jInd))
    for j = jInd
        for i = 1:size(A, 1)
            queueUpdate(B, i, j, A[i,j])
        end
    end
    processQueues(B)
    return B
end

# FixMe! Should this one handle vectors of matrices?
function hcat(x::Vector{DistMatrix{T}}) where {T}
    l    = length(x)
    if l == 0
        throw(ArgumentError("cannot flatten empty vector"))
    else
        x1   = x[1]
        m, n = size(x1, 1), size(x1, 2)
        if n != 1
            throw(ArgumentError("elements have to be vectors, i.e. the second dimension has to have size one"))
        end
        A    = DistMatrix(T)
        zeros!(A, m, l*n)
        for j = 1:l
            xj = x[j]
            for k = 1:localWidth(xj)
                for i = 1:localHeight(xj)
                    xji = getLocal(xj, i, 1)
                    queueUpdate(A, globalRow(xj, i), j, xji)
                end
            end
        end
        processQueues(A)
        return A
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE

"""
    project_from_slave_to_master(Val{:Seg2}, xs, ns, xm1, xm2)

Find the projection of a slave node `xs`, having normal vector `ns`, onto master
elements with nodes (`xm1`, `xm2`). Returns master element dimensionless parameter
xi, that is,

    xm = 1/2*(1-xi)*xm1 + 1/2*(1+xi)*xm2

# Example

```jldoctest
xm1 = [7.0, 2.0]
xm2 = [4.0, -2.0]
xs = [0.0, 0.0]
ns = [3.0/5.0, -4.0/5.0]
xi2 = project_from_slave_to_master(Val{:Seg2}, xs, ns, xm1, xm2)
round(xi2, digits=6)

# output

1.833333
```

"""
function project_from_slave_to_master(::Type{Val{:Seg2}}, xs, ns, xm1, xm2)
    nom = ns[1]*(xm1[2] + xm2[2] - 2*xs[2]) - ns[2]*(xm1[1] + xm2[1] - 2*xs[1])
    denom = ns[1]*(xm1[2] - xm2[2]) + ns[2]*(xm2[1] - xm1[1])
    return nom/denom
end

"""
    project_from_master_to_slave(Val{:Seg2}, xm, xs1, xs2, ns1, ns2)

Find the projection of a master node `xm`, to the slave surface with nodes
(`xs1`, `xs2`), in direction of slave surface normal defined by (`ns1`, `ns2`).
Returns slave element dimensionless parameter, that is, to find coordinates
in slave side:

    xs = 1/2*(1-xi)*xs1 + 1/2*(1+xi)*xs2

# Example

```jldoctest
xm = [4.0, -2.0]
xs1 = [0.0, 0.0]
xs2 = [4.0, 3.0]
ns1 = [3.0/5.0, -4.0/5.0]
ns2 = sqrt(2)/2*[1, -1]
xi1 = project_from_master_to_slave(Val{:Seg2}, xm, xs1, xs2, ns1, ns2)
round(xi1, digits=6)

# output

-0.281575
```

"""
function project_from_master_to_slave(::Type{Val{:Seg2}}, xm, xs1, xs2, ns1, ns2)

    # Special case when normal is constant. Then we have
    # unique solution and linear equation to solve.
    if isapprox(ns1, ns2)
        nom = -2*ns1[1]*xm[2] + ns1[1]*xs1[2] + ns1[1]*xs2[2] + 2*ns1[2]*xm[1] - ns1[2]*xs1[1] - ns1[2]*xs2[1]
        denom = ns1[1]*xs1[2] - ns1[1]*xs2[2] - ns1[2]*xs1[1] + ns1[2]*xs2[1]
        return nom/denom
    end

    nom_b = 2*ns1[1]*xm[2] - 2*ns1[1]*xs1[2] - 2*ns1[2]*xm[1] + 2*ns1[2]*xs1[1] - 2*ns2[1]*xm[2] + 2*ns2[1]*xs2[2] + 2*ns2[2]*xm[1] - 2*ns2[2]*xs2[1]
    denom_b = ns1[1]*xs1[2] - ns1[1]*xs2[2] - ns1[2]*xs1[1] + ns1[2]*xs2[1] - ns2[1]*xs1[2] + ns2[1]*xs2[2] + ns2[2]*xs1[1] - ns2[2]*xs2[1]
    nom_c = -2*ns1[1]*xm[2] + ns1[1]*xs1[2] + ns1[1]*xs2[2] + 2*ns1[2]*xm[1] - ns1[2]*xs1[1] - ns1[2]*xs2[1] - 2*ns2[1]*xm[2] + ns2[1]*xs1[2] + ns2[1]*xs2[2] + 2*ns2[2]*xm[1] - ns2[2]*xs1[1] - ns2[2]*xs2[1]
    denom_c = ns1[1]*xs1[2] - ns1[1]*xs2[2] - ns1[2]*xs1[1] + ns1[2]*xs2[1] - ns2[1]*xs1[2] + ns2[1]*xs2[2] + ns2[2]*xs1[1] - ns2[2]*xs2[1]

    a = 1.0
    b = nom_b / denom_b
    c = nom_c / denom_c
    d = b^2 - 4*a*c
    if d < 0
        @warn("Mortar2D.calculate_projections(): negative discriminant $d")
        @warn("xm = $xm, xs1 = $xs1, xs2 = $xs2, ns1 = $ns1, ns2 = $ns2")
    end
    sols = [-b + sqrt(d), -b - sqrt(d)]/(2.0*a)
    return sols[argmin(abs.(sols))]
end