text:
# Types for module InterestRates

"""
The type `DayCountConvention` sets the convention on how to count the number of days between dates, and also how to convert that number of days into a year fraction.

Given an initial date `D1` and a final date `D2`, here's how the distance between `D1` and `D2` are mapped into a year fraction for each supported day count convention:

* *Actual360* : `(D2 - D1) / 360`
* *Actual365* : `(D2 - D1) / 365`
* *BDays252* : `bdays(D1, D2) / 252`, where `bdays` is the business days between `D1` and `D2` from [BusinessDays.jl package](https://github.com/felipenoris/BusinessDays.jl).
"""
abstract DayCountConvention
type Actual360 <: DayCountConvention end
type Actual365 <: DayCountConvention end

type BDays252 <: DayCountConvention
	hc::HolidayCalendar
end

"""
The type `CompoundingType` sets the convention on how to convert a yield into an Effective Rate Factor.

Given a yield `r` and a maturity year fraction `t`, here's how each supported compounding type maps the yield to Effective Rate Factors:

* *ContinuousCompounding* : `exp(r*t)`
* *SimpleCompounding* : `(1+r*t)`
* *ExponentialCompounding* : `(1+r)^t`
"""
abstract CompoundingType
type ContinuousCompounding <: CompoundingType end   # exp(r*t)
type SimpleCompounding <: CompoundingType end       # (1+r*t)
type ExponentialCompounding <: CompoundingType end  # (1+r)^t

"""
This package provides the following curve methods.

**Interpolation Methods**

* **Linear**: provides Linear Interpolation on rates.
* **FlatForward**: provides Flat Forward interpolation, which is implemented as a Linear Interpolation on the *log* of discount factors.
* **StepFunction**: creates a step function around given data points.
* **CubicSplineOnRates**: provides *natural cubic spline* interpolation on rates.
* **CubicSplineOnDiscountFactors**: provides *natural cubic spline* interpolation on discount factors.
* **CompositeInterpolation**: provides support for different interpolation methods for: (1) extrapolation before first data point (`before_first`), (2) interpolation between the first and last point (`inner`), (3) extrapolation after last data point (`after_last`).

For *Interpolation Methods*, the field `dtm` holds the number of days between `date` and the maturity of the observed yield, following the curve's day count convention, which must be given in advance, when creating an instance of the curve. The field `zero_rates` holds the yield values for each maturity provided in `dtm`. All yields must be anual based, and must also be given in advance, when creating the instance of the curve.

**Term Structure Models**

* **NelsonSiegel**: term structure model based on *Nelson, C.R., and A.F. Siegel (1987), Parsimonious Modeling of Yield Curve, The Journal of Business, 60, 473-489*.
* **Svensson**: term structure model based on *Svensson, L.E. (1994), Estimating and Interpreting Forward Interest Rates: Sweden 1992-1994, IMF Working Paper, WP/94/114*.

For *Term Structure Models*, the field `parameters` holds the constants defined by each model, as described below. They must be given in advance, when creating the instance of the curve.

For **NelsonSiegel** method, the array `parameters` holds the following parameters from the model:
* **beta1** = parameters[1]
* **beta2** = parameters[2]
* **beta3** = parameters[3]
* **lambda** = parameters[4]

For **Svensson** method, the array `parameters` hold the following parameters from the model:
* **beta1** = parameters[1]
* **beta2** = parameters[2]
* **beta3** = parameters[3]
* **beta4** = parameters[4]
* **lambda1** = parameters[5]
* **lambda2** = parameters[6]

### Methods hierarchy

As a summary, curve methods are organized by the following hierarchy.

* `<<CurveMethod>>`
	* `<<Interpolation>>`
		* `<<DiscountFactorInterpolation>>`
			* `CubicSplineOnDiscountFactors`
			* `FlatForward`
		* `<<RateInterpolation>>`
			* `CubicSplineOnRates`
			* `Linear`
			* `StepFunction`
		* `CompositeInterpolation`
	* `<<Parametric>>`
		* `NelsonSiegel`
		* `Svensson`
"""
abstract CurveMethod
abstract Parametric <: CurveMethod
abstract Interpolation <: CurveMethod

abstract DiscountFactorInterpolation <: Interpolation
abstract RateInterpolation <: Interpolation

type CubicSplineOnRates <: RateInterpolation end
type CubicSplineOnDiscountFactors <: DiscountFactorInterpolation end
type FlatForward <: DiscountFactorInterpolation end
type Linear <: RateInterpolation end
type NelsonSiegel <: Parametric end
type Svensson <: Parametric end
type StepFunction <: RateInterpolation end

type CompositeInterpolation <: Interpolation
	before_first::Interpolation # Interpolation method to be applied before the first point
	inner::Interpolation
	after_last::Interpolation # Interpolation method to be applied after the last point
end

abstract AbstractIRCurve

type IRCurve <: AbstractIRCurve
	name::String
	daycount::DayCountConvention
	compounding::CompoundingType
	method::CurveMethod
	date::Date
	dtm::Vector{Int} # for interpolation methods, stores days_to_maturity on curve's daycount convention.
	zero_rates::Vector{Float64} # for interpolation methods, parameters[i] stores yield for maturity dtm[i].
	parameters::Vector{Float64} # for parametric methods, parameters stores model's constant parameters.
	dict::Dict{Symbol, Any}		# holds pre-calculated values for optimization, or additional parameters.

	# Constructor for Interpolation methods
	IRCurve{M<:Interpolation}(name::AbstractString, _daycount::DayCountConvention,
		compounding::CompoundingType, method::M,
		date::Date, dtm::Vector{Int},
		zero_rates::Vector{Float64}, parameters = Array(Float64,0), dict = Dict{Symbol, Any}()) = begin

		isempty(dtm) && error("Empty days-to-maturity vector")
		isempty(zero_rates) && error("Empty zero_rates vector")
		(length(dtm) != length(zero_rates)) && error("dtm and zero_rates must have the same length")
		(!issorted(dtm)) && error("dtm should be sorted before creating IRCurve instance")

		new(String(name), _daycount, compounding, method, date, dtm, zero_rates, parameters, dict)
	end

	# Constructor for Parametric methods
	IRCurve{M<:Parametric}(name::AbstractString, _daycount::DayCountConvention,
		compounding::CompoundingType, method::M,
		date::Date,
		parameters::Vector{Float64},
		dict = Dict{Symbol, Any}()) = begin
		isempty(parameters) && error("Empty yields vector")
		new(String(name), _daycount, compounding, method, date, Array(Int,0), Array(Float64,0), parameters, dict)
	end
end

# Interface for curve types
"""
	curve_get_name(curve::AbstractIRCurve) → String

Returns the name of the curve.
"""
curve_get_name(curve::AbstractIRCurve) = error("method not defined")

"""
	curve_get_daycount(curve::AbstractIRCurve) → DayCountConvention

Returns the DayCountConvention used by the curve. See DayCountConvention documentation.
"""
curve_get_daycount(curve::AbstractIRCurve) = error("method not defined")

"""
	curve_get_compounding(curve::AbstractIRCurve) → CompoundingType

Returns the CompoundingType used by the curve. See CompoundingType documentation.
"""
curve_get_compounding(curve::AbstractIRCurve) = error("method not defined")

"""
	curve_get_method(curve::AbstractIRCurve) → CurveMethod

Returns the CurveMethod used by the curve. See CurveMethod documentation.
"""
curve_get_method(curve::AbstractIRCurve) = error("method not defined")

"""
	curve_get_date(curve::AbstractIRCurve) → Date

Returns the date when the curve is observed. All zero rate calculation will be performed based on this date.
"""
curve_get_date(curve::AbstractIRCurve) = error("method not defined")

"""
	curve_get_dtm(curve::AbstractIRCurve) → Vector{Int}

Used for interpolation methods, returns `days_to_maturity` on curve's daycount convention.
"""
curve_get_dtm(curve::AbstractIRCurve) = error("method not defined")

"""
	curve_get_zero_rates(curve::AbstractIRCurve) → Vector{Float64}

Used for interpolation methods, parameters[i] returns yield for maturity dtm[i].
"""
curve_get_zero_rates(curve::AbstractIRCurve) = error("method not defined")

"""
	curve_get_model_parameters(curve::AbstractIRCurve) → Vector{Float64}

Used for parametric methods, returns model's constant parameters.
"""
curve_get_model_parameters(curve::AbstractIRCurve) = error("method not defined")

# AbstractIRCurve interface implementation for IRCurve type
curve_get_name(curve::IRCurve) = curve.name
curve_get_daycount(curve::IRCurve) = curve.daycount
curve_get_compounding(curve::IRCurve) = curve.compounding
curve_get_method(curve::IRCurve) = curve.method
curve_get_date(curve::IRCurve) = curve.date
curve_get_dtm(curve::IRCurve) = curve.dtm
curve_get_zero_rates(curve::IRCurve) = curve.zero_rates
curve_get_model_parameters(curve::IRCurve) = curve.parameters
curve_get_dict_parameter(curve::IRCurve, sym::Symbol) = curve.dict[sym]

function curve_set_dict_parameter!(curve::IRCurve, sym::Symbol, value)
	curve.dict[sym] = value
end


text:
using Knet: Knet, conv4, pool, mat, nll, accuracy, progress, sgd, param, param0, dropout, relu, bnmoments, bnparams
using Knet: bnmoments, bnparams, deconv4,tanh

#=
Conv3d
=#
struct Conv3d;
    w;
    b;
    padding;
    stride;
end

(c::Conv3d)(x) = conv4(c.w, x; padding=c.padding, stride=c.stride, mode=1) .+ c.b # size(x) = (a, b, c, out_channels, batch)
Conv3d(in_channels::Int, out_channels::Int, kernel_size::Int, padding::Int, stride::Int) =
    Conv3d(param(kernel_size, kernel_size, kernel_size, in_channels, out_channels),
    param0(1,1,1,out_channels,1), padding, stride)

#=
BatchNorm3d
=#
struct BatchNorm3d;
    moments;
    params;
end

#=
Conv3d + batchnorm3d + relu
=#
(c::BatchNorm3d)(x) = batchnorm(x, c.moments, c.params) # size(x) = (a, b, c, out_channels, batch)
BatchNorm3d(out_channels::Int) = BatchNorm3d(bnmoments(), bnparams(out_channels))

struct BatchConv3d;
    conv3d;
    batch3d;
end

(c::BatchConv3d)(x) = relu.(c.batch3d(c.conv3d(x)))
BatchConv3d(in_channels::Int, out_channels::Int, kernel_size::Int, padding::Int, stride::Int) =
    BatchConv3d(Conv3d(in_channels, out_channels, kernel_size, padding, stride), BatchNorm3d(out_channels))

#=
Conv3d + relu
=#
struct RConv3d;
    conv3d;
end

(c::RConv3d)(x) = relu.(c.conv3d(x))
RConv3d(in_channels::Int, out_channels::Int, kernel_size::Int, padding::Int, stride::Int) =
    RConv3d(Conv3d(in_channels, out_channels, kernel_size, padding, stride))

#=
Sequential
=#
struct Sequential;
    layers;
end

function (c::Sequential)(x)
    for layer in c.layers
        x = layer(x)
    end
    x
end

#=
Sequence of Conv3d + relu
=#
struct RConvBlock;
    sequential;
end

function (c::RConvBlock)(x)
    x = c.sequential(x)
end

function RConvBlock(n_stages::Int, n_filters_in::Int, n_filters_out::Int)
    seq = [RConv3d(n_filters_in, n_filters_out, 3, 1, 1)]
    for i=2:n_stages
        push!(seq, RConv3d(n_filters_out, n_filters_out, 3, 1, 1))
    end
    return RConvBlock(Sequential(seq))
end

#=
Sequence of Conv3d + batchnorm3d + relu
=#
struct BatchConvBlock;
    sequential;
end

function (c::BatchConvBlock)(x)
    x = c.sequential(x)
end

function BatchConvBlock(n_stages::Int, n_filters_in::Int, n_filters_out::Int)
    seq = [BatchConv3d(n_filters_in, n_filters_out, 3, 1, 1)]
    for i=2:n_stages
        push!(seq, BatchConv3d(n_filters_out, n_filters_out, 3, 1, 1))
    end
    return BatchConvBlock(Sequential(seq))
end

#=
Residual Sequence of Conv3d + batchnorm3d + relu
=#
struct ResidualBatchConvBlock;
    sequential;
    conv3d;
    batchnorm3d;
end

function (c::ResidualBatchConvBlock)(x)
    y = c.sequential(x)
    y = c.conv3d(y)
    y = c.batchnorm3d(y)
    x = x .+ y
    x = relu.(x)
    return x
end

function ResidualBatchConvBlock(n_stages::Int, n_filters_in::Int, n_filters_out::Int)
    seq = [BatchConv3d(n_filters_in, n_filters_out, 3, 1, 1)]
    for i=2:n_stages-1
        push!(seq, BatchConv3d(n_filters_out, n_filters_out, 3, 1, 1))
    end
    return ResidualBatchConvBlock(Sequential(seq), Conv3d(n_filters_out, n_filters_out, 3, 1, 1), BatchNorm3d(n_filters_out))
end

#=
Residual Sequence of Conv3d + relu
=#
struct ResidualRConvBlock;
    sequential;
    conv3d;
end

function (c::ResidualRConvBlock)(x)
    y = c.sequential(x)
    y = c.conv3d(y)
    x = x .+ y
    x = relu.(x)
    return x
end

function ResidualRConvBlock(n_stages::Int, n_filters_in::Int, n_filters_out::Int)
    seq = [RConv3d(n_filters_in, n_filters_out, 3, 1, 1)]
    for i=2:n_stages-1
        push!(seq, RConv3d(n_filters_out, n_filters_out, 3, 1, 1))
    end
    return ResidualRConvBlock(Sequential(seq), Conv3d(n_filters_out, n_filters_out, 3, 1, 1))
end

#=
Downsampling layer: Conv3d + relu
=#
struct DownsamplingRConv3d;
    conv3d;
end

(c::DownsamplingRConv3d)(x) = relu.(c.conv3d(x))
DownsamplingRConv3d(in_channels::Int, out_channels::Int, stride::Int) =
    DownsamplingRConv3d(Conv3d(in_channels, out_channels, stride, 0, stride))

#=
Downsampling layer: Conv3d + batchnorm3d + relu
=#
struct DownsamplingBatchConv3d;
    conv3d;
    batch3d;
end

(c::DownsamplingBatchConv3d)(x) = relu.(c.batch3d(c.conv3d(x)))
DownsamplingBatchConv3d(in_channels::Int, out_channels::Int, stride::Int) =
    DownsamplingBatchConv3d(Conv3d(in_channels, out_channels, stride, 0, stride), BatchNorm3d(out_channels))

#=
Transposed convolution
=#
struct ConvTranspose3d;
    w;
    b;
    padding;
    stride;
end

(c::ConvTranspose3d)(x) = deconv4(c.w, x; padding=c.padding, stride=c.stride, mode=1) .+ c.b # size(x) = (a, b, c, out_channels, batch)
ConvTranspose3d(in_channels::Int, out_channels::Int, kernel_size::Int, padding::Int, stride::Int) =
    ConvTranspose3d(param(kernel_size, kernel_size, kernel_size, out_channels, in_channels), param0(1,1,1,out_channels,1), padding, stride)

#=
Sequence of ConvTranspose3d + relu
=#
struct RConvTranspose3d;
    convtranspose3d;
end

(c::RConvTranspose3d)(x) = relu.(c.convtranspose3d(x))
RConvTranspose3d(in_channels::Int, out_channels::Int, kernel_size::Int, padding::Int, stride::Int) =
    RConvTranspose3d(ConvTranspose3d(in_channels, out_channels, kernel_size, padding, stride))

#=
Sequence of ConvTranspose3d + batchnorm3d + relu
=#
struct BatchConvTranspose3d;
    convtranspose3d;
    batch3d;
end

(c::BatchConvTranspose3d)(x) = relu.(c.batch3d(c.convtranspose3d(x)))
BatchConvTranspose3d(in_channels::Int, out_channels::Int, kernel_size::Int, padding::Int, stride::Int) =
    BatchConvTranspose3d(ConvTranspose3d(in_channels, out_channels, kernel_size, padding, stride), BatchNorm3d(out_channels))

#=

=#
struct RConvTranspose3d;
    convtranspose3d;
end

(c::RConvTranspose3d)(x) = relu.(c.convtranspose3d(x))
RConvTranspose3d(in_channels::Int, out_channels::Int, kernel_size::Int, padding::Int, stride::Int) =
    RConvTranspose3d(ConvTranspose3d(in_channels, out_channels, kernel_size, padding, stride))

#=

=#
struct UpsamplingRConvTranspose3d;
    convtranspose3d;
end

(c::UpsamplingRConvTranspose3d)(x) = relu.(c.convtranspose3d(x))
UpsamplingRConvTranspose3d(in_channels::Int, out_channels::Int, stride::Int) =
    UpsamplingRConvTranspose3d(ConvTranspose3d(in_channels, out_channels, stride, 0, stride))

ext: jl

max_stars_repo_name: skoc/julia-project

max_issues_repo_name: skoc/julia-project

max_forks_repo_name: skoc/julia-project


text:
import Test
import Random
import Statistics
import Zygote
import RestrictedBoltzmannMachines as RBMs
import CenteredRBMs
using CenteredRBMs: CenteredBinaryRBM, center, uncenter
using RestrictedBoltzmannMachines: visible, hidden, weights

@testset "center / uncenter" begin
    rbm = RBMs.BinaryRBM(randn(3), randn(2), randn(3,2))
    offset_v = randn(3)
    offset_h = randn(2)
    centered_rbm = center(rbm, offset_v, offset_h)
    @test centered_rbm.offset_v ≈ offset_v
    @test centered_rbm.offset_h ≈ offset_h
    @test visible(uncenter(centered_rbm)).θ ≈ visible(rbm).θ
    @test hidden(uncenter(centered_rbm)).θ ≈ hidden(rbm).θ
    @test weights(uncenter(centered_rbm)) ≈ weights(rbm)

    v = Random.bitrand(3,2)
    @test RBMs.mean_h_from_v(rbm, v) ≈ RBMs.mean_h_from_v(uncenter(rbm), v)
    h = Random.bitrand(2,2)
    @test RBMs.mean_v_from_h(rbm, h) ≈ RBMs.mean_v_from_h(uncenter(rbm), h)
end

@testset "rbm energy invariance" begin
    centered_rbm = CenteredBinaryRBM(randn(3), randn(2), randn(3,2), randn(3), randn(2))
    rbm = uncenter(centered_rbm)
    ΔE = RBMs.interaction_energy(rbm, centered_rbm.offset_v, centered_rbm.offset_h)::Number
    v = Random.bitrand(size(visible(rbm))..., 100)
    h = Random.bitrand(size(hidden(rbm))..., 100)
    @test RBMs.energy(centered_rbm, v, h) ≈ RBMs.energy(rbm, v, h) .+ ΔE
    @test RBMs.free_energy(centered_rbm, v) ≈ RBMs.free_energy(rbm, v) .+ ΔE
end

@testset "free energy" begin
    rbm = CenteredBinaryRBM(randn(3), randn(2), randn(3,2), randn(3), randn(2))
    v = Random.bitrand(size(visible(rbm))...)
    F = -log(sum(exp(-RBMs.energy(rbm, v, h)) for h in [[0,0], [0,1], [1,0], [1,1]]))
    @test RBMs.free_energy(rbm, v) ≈ F
end

@testset "∂free energy" begin
    rbm = CenteredBinaryRBM(randn(3), randn(2), randn(3,2), randn(3), randn(2))
    v = Random.bitrand(size(visible(rbm))...)
    gs = Zygote.gradient(rbm) do rbm
        Statistics.mean(RBMs.free_energy(rbm, v))
    end
    ∂ = RBMs.∂free_energy(rbm, v)
    @test ∂.visible.θ ≈ only(gs).rbm.visible.θ
    @test ∂.hidden.θ ≈ only(gs).rbm.hidden.θ
    @test ∂.w ≈ only(gs).rbm.w
end


text:
function DILocation(ctx, line, col, scope=nothing, inlined_at=nothing)
    # XXX: are null scopes valid? they crash LLVM:
    #      DILocation(Context(), 1, 2).scope
    DILocation(API.LLVMDIBuilderCreateDebugLocation(ctx, line, col,
                                                    something(scope, C_NULL),
                                                    something(inlined_at, C_NULL)))
end


text:
"""
# References
- Reference model (e.g., xd, vd, ad, ad_dot, ad_ddot)
[1] S. J. Su, Y. Y. Zhu, H. R. Wang, and C. Yun, “A Method to Construct a Reference Model for Model Reference Adaptive Control,” Adv. Mech. Eng., vol. 11, no. 11, pp. 1–9, 2019.
# Notes
d: degree of the reference model
auto_diff: auto_diff = true for tracking time-varying command; otherwise for set-point regulation
x_cmd_func(t): function of time
"""
struct ReferenceModel <: AbstractEnv
    d::Int  # degree
    Ks::AbstractArray
    auto_diff::Bool
    x_cmd_func::Union{Function, Nothing}
end
function ReferenceModel(d::Int; x_cmd_func=nothing)
    @assert d >= 0
    Ks = []
    if d == 4
        # [1]
        push!(Ks, Diagonal(1.0*ones(3)))
        push!(Ks, Diagonal(3.4*ones(3)))
        push!(Ks, Diagonal(5.4*ones(3)))
        push!(Ks, Diagonal(4.9*ones(3)))
        push!(Ks, Diagonal(2.7*ones(3)))
    else
        error("Assign values of matrix `Kx` manually")
    end
    auto_diff = x_cmd_func == nothing ? false : true
    ReferenceModel(d, Ks, auto_diff, x_cmd_func)
end

"""
`x_i` denotes `i`th (time) derivative of `x`
"""
function State(env::ReferenceModel)
    @unpack d = env
    return function (x0)
        xs_dict = Dict(:x_0 => x0)
        for i in 1:d
            xs_dict[Symbol(:x_, i)] = zeros(size(x0))  # e.g., xs_dict[:x1] = zeros(3)
        end
        ComponentArray((; zip(keys(xs_dict), values(xs_dict))...))
    end
end

function Dynamics!(env::ReferenceModel)
    @unpack d, Ks, auto_diff, x_cmd_func = env
    # derivatives for auto_diff
    funcs = nothing
    if auto_diff
        _funcs = Function[x_cmd_func]
        for i in 1:d+1
            push!(_funcs, (t) -> ForwardDiff.derivative(_funcs[i], t))
        end
    end
    funcs(t) = [_func(t) for _func in _funcs]
    @Loggable function dynamics!(dX, X, p, t; x_cmd=nothing)
        @onlylog state = X
        xs = [getproperty(X, Symbol(:x_, 0))]  # x_0, x_1, ...
        for i in 0:d-1
            _x_next = getproperty(X, Symbol(:x_, i+1))
            setproperty!(dX, Symbol(:x_, i), _x_next)  # e.g., dX.x_0 = x_1
            push!(xs, _x_next)
        end
        dx_d = nothing
        if auto_diff
            if x_cmd != nothing
                error("Do not provide a manual command for mode auto_diff=$(auto_diff)")
            end
            # dx_d = -sum(Ks .* xs) + sum([Ks..., I] .* funcs(t))
            desired_cmds = funcs(t)
            dx_d = -sum(Ks .* xs) + sum([Ks..., I] .* desired_cmds)
            @log x_cmd = desired_cmds[1]
        else
            if x_cmd == nothing
                error("Provide a manual command for mode auto_diff=$(auto_diff)")
            end
            dx_d = -sum(Ks .* xs) + Ks[1]*x_cmd
            @log x_cmd
        end
        setproperty!(dX, Symbol(:x_, d), dx_d)
    end
end


text:
using Test
using MianalyzerBackend
const HierachicalClustering = MianalyzerBackend.HierachicalClustering

@testset "HierachicalClustering" begin
  @testset ".call" begin
    actual = HierachicalClustering.fit([1 1; 1 1; 3 3; 3 3; 10 10],
                                       ["label1", "label2", "label3", "label4", "label5"],
                                       :euclidean,
                                       :ward)
    @test actual == [[["label1", "label2"], ["label3", "label4"]], "label5"]
  end
end


text:
using ITensors
using ITensorInfiniteMPS

##############################################################################
# VUMPS parameters
#

maxdim = 10 # Maximum bond dimension
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 20 # Maximum number of iterations of the VUMPS algorithm at a fixed bond dimension
tol = 1e-5 # Precision error tolerance for outer loop of VUMPS or TDVP
outer_iters = 5 # Number of times to increase the bond dimension
time_step = -Inf # -Inf corresponds to VUMPS
solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP)
multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]
conserve_qns = true
N = 2 # Number of sites in the unit cell (1 site unit cell is currently broken)

# Parameters of the transverse field Ising model
model_params = (J=1.0, h=0.9)

##############################################################################
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
#

model = Model"ising"()

function space_shifted(::Model"ising", q̃sz; conserve_qns=true)
  if conserve_qns
    return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
  else
    return [QN() => 2]
  end
end

space_ = fill(space_shifted(model, 0; conserve_qns=conserve_qns), N)
s = infsiteinds("S=1/2", N; space=space_)
initstate(n) = "↑"
ψ = InfMPS(s, initstate)

# Form the Hamiltonian
H = InfiniteSum{MPO}(model, s; model_params...)

# Check translational invariance
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

vumps_kwargs = (
  tol=tol,
  maxiter=max_vumps_iters,
  solver_tol=solver_tol,
  multisite_update_alg=multisite_update_alg,
)
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
#ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)

# Alternate steps of running VUMPS and increasing the bond dimension
@time for _ in 1:outer_iters
  println("\nIncrease bond dimension")
  global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
  println("Run VUMPS with new bond dimension")
  global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
end

# Check translational invariance
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

#
# Compare to DMRG
#

Nfinite = 100
sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
Hfinite = MPO(model, sfinite; model_params...)
ψfinite = randomMPS(sfinite, initstate)
@show flux(ψfinite)
sweeps = Sweeps(10)
setmaxdim!(sweeps, maxdim)
setcutoff!(sweeps, cutoff)
energy_finite_total, ψfinite = @time dmrg(Hfinite, ψfinite, sweeps)
@show energy_finite_total / Nfinite

function energy_local(ψ1, ψ2, h)
  ϕ = ψ1 * ψ2
  return (noprime(ϕ * h) * dag(ϕ))[]
end

function ITensors.expect(ψ, o)
  return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[]
end

# Exact energy at criticality: 4/pi = 1.2732395447351628

nfinite = Nfinite ÷ 2
orthogonalize!(ψfinite, nfinite)
hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_params...)
energy_finite = energy_local(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite)
energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], H[(1, 2)])
@show energy_finite, energy_infinite
@show abs(energy_finite - energy_infinite)

energy_exact = reference(model, Observable("energy"); model_params...)
@show energy_exact

Sz1_finite = expect(ψfinite[nfinite], "Sz")
orthogonalize!(ψfinite, nfinite + 1)
Sz2_finite = expect(ψfinite[nfinite + 1], "Sz")
Sz1_infinite = expect(ψ.AL[1] * ψ.C[1], "Sz")
Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz")

@show Sz1_finite, Sz2_finite
@show Sz1_infinite, Sz2_infinite

##############################################################################
# Compute eigenspace of the transfer matrix
#

using Arpack
using KrylovKit
using LinearAlgebra

T = TransferMatrix(ψ.AL)
Tᵀ = transpose(T)
vⁱᴿ = randomITensor(dag(input_inds(T)))
vⁱᴸ = randomITensor(dag(input_inds(Tᵀ)))

neigs = 10
tol = 1e-10
λ⃗ᴿ, v⃗ᴿ, right_info = eigsolve(T, vⁱᴿ, neigs, :LM; tol=tol)
λ⃗ᴸ, v⃗ᴸ, left_info = eigsolve(Tᵀ, vⁱᴸ, neigs, :LM; tol=tol)

@show norm(T(v⃗ᴿ[1]) - λ⃗ᴿ[1] * v⃗ᴿ[1])
@show norm(Tᵀ(v⃗ᴸ[1]) - λ⃗ᴸ[1] * v⃗ᴸ[1])

@show λ⃗ᴿ
@show λ⃗ᴸ
@show flux.(v⃗ᴿ)

neigs = length(v⃗ᴿ)

# Normalize the vectors
N⃗ = [(translatecell(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs]

v⃗ᴿ ./= sqrt.(N⃗)
v⃗ᴸ ./= sqrt.(N⃗)

# Form a second starting vector orthogonal to v⃗ᴿ[1]
# This doesn't work. TODO: project out v⃗ᴿ[1], v⃗ᴸ[1] from T
#λ⃗ᴿ², v⃗ᴿ², right_info_2 = eigsolve(T, vⁱᴿ², neigs, :LM; tol=tol)

# Projector onto the n-th eigenstate
function proj(v⃗ᴸ, v⃗ᴿ, n)
  Lⁿ = v⃗ᴸ[n]
  Rⁿ = v⃗ᴿ[n]
  return ITensorMap(
    [translatecell(Lⁿ, 1), translatecell(Rⁿ, -1)]; input_inds=inds(Rⁿ), output_inds=inds(Lⁿ)
  )
end

P⃗ = [proj(v⃗ᴸ, v⃗ᴿ, n) for n in 1:neigs]
T⁻P = T - sum(P⃗)

#vⁱᴿ² = vⁱᴿ - (translatecell(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1]
#@show norm(dag(vⁱᴿ²) * v⃗ᴿ[1])

λ⃗ᴾᴿ, v⃗ᴾᴿ, right_info = eigsolve(T⁻P, vⁱᴿ, neigs, :LM; tol=tol)
@show λ⃗ᴾᴿ

vⁱᴿ⁻ᵈᵃᵗᵃ = vec(array(vⁱᴿ))
λ⃗ᴿᴬ, v⃗ᴿ⁻ᵈᵃᵗᵃ = Arpack.eigs(T; v0=vⁱᴿ⁻ᵈᵃᵗᵃ, nev=neigs)

## XXX: this is giving an error about trying to set the element of the wrong QN block for:
## maxdim = 5
## cutoff = 1e-12
## max_vumps_iters = 10
## outer_iters = 10
## model_params = (J=1.0, h=0.8)
##
## v⃗ᴿᴬ = [itensor(v⃗ᴿ⁻ᵈᵃᵗᵃ[:, n], input_inds(T); tol=1e-4) for n in 1:length(λ⃗ᴿᴬ)]
## @show flux.(v⃗ᴿᴬ)

@show λ⃗ᴿᴬ

# Full eigendecomposition

Tfull = prod(T)
DV = eigen(Tfull, input_inds(T), output_inds(T))

@show norm(Tfull * DV.V - DV.Vt * DV.D)

d = diag(array(DV.D))

p = sortperm(d; by=abs, rev=true)
@show p[1:neigs]
@show d[p[1:neigs]]

println("Error if ED with Arpack")
@show d[p[1:neigs]] - λ⃗ᴿᴬ


text:
#####################
## model parameters
#####################

p = (
    γ = 0.5,    # labor share
    θ = 8.,     # trade elasticity
    damp = 0.1, # function iteration damping
    tol = 1e-3  # convergence tolerance
    )

#####################
## initial data
#####################

λ = [.3 .7;
     .4 .6]         # trade shares
Y = [100., 100.]    # output
X = λ .* p.γ .* Y   # bilateral expenditures
π = [.8 .2;
     .3 .7]         # migration shares
L = [1.5, 1.5]        # labor distribution

#####################
## counterfactual set up
#####################

T_hat = [5., 1.]  # counterfactual change in productivity

#####################
## initialize model
#####################

w_guess = ones(2) # initial guess

data = (λ = λ, π = π, L = L, Y = Y, X = X) # set up data tuple

# make sure shares sum to 1
sum(data.λ, dims = 2)
sum(data.π, dims = 2)

#####################
## solver functions
#####################

# solve P as function of wages and data
function solve_P_hat(p, data, w_hat, T_hat)
    ( data.λ * ( T_hat .* w_hat.^(-p.γ*p.θ) ) ).^(-1/p.θ)
end

# solve λ as function of wages and P
function solve_λ_hat(p, w_hat, T_hat, P_hat)
    T_hat .* ( w_hat.^p.γ ./ P_hat ).^-p.θ
end

# solve for migration shares as function of real wages and data
function solve_π_hat(p, data, w_hat, P_hat)
    (w_hat ./ P_hat) ./ (data.π * ((w_hat ./ P_hat)))
end

# solve for labor as function of migration shares and data
function solve_L_hat(p, data, π_hat)
    (1 ./ data.L) .* (( π_hat .* data.π ) * data.L)
end

# solve for change in wages (and other variables)
function solve_w_hat(p, data, w_hat)

    L_hat = similar(w_hat)
    π_hat = similar(w_hat)
    λ_hat = similar(w_hat)
    P_hat = similar(w_hat)
    error = Inf
    iters = 1

    while error > p.tol && iters < 1e3

        P_hat = solve_P_hat(p, data, w_hat, T_hat)  # given wage guess, solve for P
        λ_hat = solve_λ_hat(p, w_hat, T_hat, P_hat) # given P, solve for trade shares
        π_hat = solve_π_hat(p, data, w_hat, P_hat)  # given real wages, solve for migration shares
        L_hat = solve_L_hat(p, data, π_hat)         # given migration shares, solve for labor

        # wage fixed point we are iterating on until convergence
        w_hat_new = (1/p.γ) * (1 ./ L_hat) .* (1 ./ data.Y) .* ( ( data.X .* λ_hat ) * ( w_hat .* L_hat ) )

        # error checking
        error = maximum(abs.((w_hat_new .- w_hat)./w_hat))
        w_hat = p.damp .* w_hat_new .+ (1 - p.damp) .* w_hat
        iters += 1

    end

    return w_hat, L_hat, π_hat, λ_hat, P_hat

end

#####################
## SOLVE MODEL
#####################

w_hat, L_hat, π_hat, λ_hat, P_hat = solve_w_hat(p, data, w_guess)

max_stars_repo_name: rafserqui/class-repo

max_issues_repo_name: rafserqui/class-repo

max_forks_repo_name: rafserqui/class-repo


text:
"""
    slope_array(w; reduction = true)

Generate array to store spatial slopes of solutions
"""
slope_array(w::Number; kwargs...) = deepcopy(w)

function slope_array(w::AbstractArray; reduction = true)
    nd = ndims(w)
    ids = []
    for i = 1:nd
        push!(ids, [axes(w, i) |> first, axes(w, i) |> last])
    end

    sw = ifelse(
        reduction == true,
        cat(zero(w), zero(w), dims = ndims(w) + 1),
        cat(zero(w), zero(w), zero(w), dims = ndims(w) + 1),
    )

    if w isa MArray
        sw = static_array(sw)
    end

    if w isa OffsetArray && typeof(w).types[1] <: MArray
        sw = static_array(sw)
        if ndims(sw) == 2
            sw = OffsetArray(sw, ids[1][1]:ids[1][2], axes(sw)[end])
        elseif ndims(sw) == 3
            sw = OffsetArray(sw, ids[1][1]:ids[1][2], ids[2][1]:ids[2][2], axes(sw)[end])
        elseif ndims(sw) == 4
            sw = OffsetArray(
                sw,
                ids[1][1]:ids[1][2],
                ids[2][1]:ids[2][2],
                ids[3][1]:ids[3][2],
                axes(sw)[end],
            )
        end
    end

    return sw
end


text:
struct ScenarioBelief{S, O, RS<:DESPOTRandomSource} <: AbstractParticleBelief{S}
    scenarios::Vector{Pair{Int,S}}
    random_source::RS
    depth::Int
    _obs::O
end

rand(rng::AbstractRNG, b::ScenarioBelief) = last(b.scenarios[rand(rng, 1:length(b.scenarios))])
ParticleFilters.particles(b::ScenarioBelief) = (last(p) for p in b.scenarios)
ParticleFilters.n_particles(b::ScenarioBelief) = length(b.scenarios)
ParticleFilters.weight(b::ScenarioBelief, s) = 1/length(b.scenarios)
ParticleFilters.particle(b::ScenarioBelief, i) = last(b.scenarios[i])
ParticleFilters.weight_sum(b::ScenarioBelief) = 1.0
ParticleFilters.weights(b::ScenarioBelief) = (1/length(b.scenarios) for p in b.scenarios)
ParticleFilters.weighted_particles(b::ScenarioBelief) = (last(p)=>1/length(b.scenarios) for p in b.scenarios)
POMDPs.pdf(b::ScenarioBelief{S}, s::S) where S = sum(p==s for p in particles(b))/length(b.scenarios) # this is slow
POMDPs.mean(b::ScenarioBelief) = mean(last, b.scenarios)
POMDPs.currentobs(b::ScenarioBelief) = b._obs
@deprecate previous_obs POMDPs.currentobs
POMDPs.history(b::ScenarioBelief) = tuple((o=currentobs(b),))

initialize_belief(::PreviousObservationUpdater, b::ScenarioBelief) = previous_obs(b)


text:
"""
    DistributedOptimiser(optimiser)

Wrap the `optimiser` in a `DistributedOptimiser`. Before updating the
parameters, this adds the gradients across the processes using
non-blocking Allreduce

!!! note
    Remember to scale the loss function by `1 / total_workers()` to ensure
    that the gradients are correctly averaged
"""
struct DistributedOptimiser{O}
    optimiser::O
end

function apply!(o::DistributedOptimiser, state, x, y)
    y_ = allreduce!(y, +, MPI.COMM_WORLD)
    apply!(o.optimiser, state, x, y_)
end

init(o::DistributedOptimiser, x::AbstractArray) = init(o.optimiser, x)

"""
    allreduce_gradients(gs::NamedTuple; on_gpu::Bool=CUDA.functional())
"""
function allreduce_gradients(gs::NamedTuple; on_gpu::Bool=CUDA.functional())
    if on_gpu
        # Transfer data to CPU since OpenMPI Iallreduce! doesn't work for CUDA
        gs = fmap(cpu, gs)
    end
    requests = MPI.Request[]
    nonblocking_reduce_gradients(g) = g
    function nonblocking_reduce_gradients(g::AbstractArray)
        g, req = Iallreduce!(g, +, MPI.COMM_WORLD)
        push!(requests, req)
        return g
    end
    gs = fmap(nonblocking_reduce_gradients, gs)
    MPI.Waitall!(requests)
    if on_gpu
        # Transfer data back to GPU
        gs = fmap(gpu, gs)
    end
    return gs
end


text:
using StokesParticles, Test
import StokesParticles: do_intersect, get_orientation, get_nearest_point!, get_nearest_point, is_in_line

p1 = Point2D(0.0,0.0); q1 = Point2D(0.0,1.0)
p2 = Point2D(-0.5,0.5); q2 = Point2D(0.5,0.5)
@testset "do_intersect tests" begin
    @test do_intersect(p1,q1,p2,q2) == true
    @test do_intersect(p1,p2,q1,q2) == false
end

z = Point2D(0.5+Float64(pi),1.5+Float64(pi))
@testset "get_orientation tests" begin
    @test get_orientation(p2,q1,z) == 0
    @test get_orientation(p2,z,q1) == 0
    @test get_orientation(z,p2,q1) == 0
    @test get_orientation(p1,q1,Point2D(1.0,0.0)) == 1
    @test get_orientation(p1,q1,Point2D(-1.0,0.0)) == 2
end

wall1 = LineWall(Point2D(0.0,0.0),Point2D(2.0,0.0),0.2)
pointOnWall1 = Point2D(0.0,0.0)
pointOnWall2 = Point2D(0.0,0.0)
pointOnWallA = Point2D(0.1,0.1)
pointOnWallC = Point2D(0.0,0.0)
particle = Particle2D(Point2D(0.6,0.1),0.5,1)
wall2 = LineWall(Point2D(0.0,0.0),Point2D(0.0,2.0),0.2)
wallC = CircleWall(Point2D(0.0,0.0),1.0,0.2)
wallA = ArcWall(Point2D(0.0,-1.0),Point2D(0.0,0.0),Point2D(0.0,1.0),0.2)
@testset "get_nearest_point! tests" begin
    get_nearest_point!(pointOnWall1,particle,wall1)
    @test pointOnWall1.x == 0.6
    @test pointOnWall1.y == 0.0

    get_nearest_point!(pointOnWall2,particle,wall2)
    @test pointOnWall2.x == 0.0
    @test pointOnWall2.y == 0.1

    pointOnWall3 = get_nearest_point(particle,wall2)
    @test pointOnWall3.x == 0.0
    @test pointOnWall3.y == 0.1

    get_nearest_point!(pointOnWallA,particle,wallA)
    get_nearest_point!(pointOnWallC,particle,wallC)
end

@testset "is_in_line tests" begin
    @test is_in_line(wall1,pointOnWall1) == true
    @test is_in_line(wall2,pointOnWall2) == true
    @test is_in_line(wall2,particle.pos) == false
    @test is_in_line(wall2,Point2D(0.0,2.1)) == false
    @test is_in_line(wall2,Point2D(0.0,-0.1)) == false
    @test is_in_line(wallA,pointOnWallA) == true
    @test is_in_line(wallC,pointOnWallC) == true
    @test is_in_line(wallA,Point2D(sqrt(2)/2,sqrt(2)/2)) == true
    @test is_in_line(wallA,Point2D(-sqrt(2)/2,sqrt(2)/2)) == false
    @test is_in_line(wallA,Point2D(sqrt(2)/2,-sqrt(2)/2)) == true
    @test is_in_line(wallA,Point2D(-sqrt(2)/2,-sqrt(2)/2)) == false
    @test is_in_line(wallA,Point2D(0.0,-0.1)) == false
    @test is_in_line(wallC,Point2D(sqrt(2)/2,sqrt(2)/2)) == true
    @test is_in_line(wallC,Point2D(0.0,-0.1)) == false
end


text:
# ===================================
# Approximations in the infinity norm
# ===================================

"""
    box_approximation(S::LazySet{N})::Union{Hyperrectangle{N}, EmptySet{N}}
        where {N<:Real}

Overapproximate a convex set by a tight hyperrectangle.

### Input

- `S` -- convex set

### Output

A tight hyperrectangle.

### Algorithm

The center of the hyperrectangle is obtained by averaging the support function
of the given set in the canonical directions, and the lengths of the sides can
be recovered from the distance among support functions in the same directions.
"""
function box_approximation(S::LazySet{N}
                          )::Union{Hyperrectangle{N}, EmptySet{N}} where {N<:Real}
    (c, r) = box_approximation_helper(S)
    if r[1] < 0
        return EmptySet{N}()
    end
    return Hyperrectangle(c, r)
end

# special case: Hyperrectangle
box_approximation(S::Hyperrectangle) = S

# special case: other rectangle
box_approximation(S::AbstractHyperrectangle) =
    Hyperrectangle(center(S), radius_hyperrectangle(S))

# special case: empty set
box_approximation(∅::EmptySet) = ∅

"""
    box_approximation(r::Rectification{N}
                     )::Union{Hyperrectangle{N}, EmptySet{N}} where {N<:Real}

Overapproximate the rectification of a convex set by a tight hyperrectangle.

### Input

- `S` -- rectification of a convex set

### Output

A hyperrectangle.

### Algorithm

Box approximation and rectification distribute.
Hence we first check whether the wrapped set is empty.
If so, we return the empty set.
Otherwise, we compute the box approximation of the wrapped set, rectify the
resulting box (which is simple), and finally convert the resulting set to a box.
"""
function box_approximation(r::Rectification{N}
                          )::Union{Hyperrectangle{N}, EmptySet{N}} where {N<:Real}
    if isempty(r.X)
        return EmptySet{N}()
    end
    return convert(Hyperrectangle, Rectification(box_approximation(r.X)))
end

"""
    interval_hull

Alias for `box_approximation`.
"""
interval_hull = box_approximation

"""
    box_approximation_symmetric(S::LazySet{N}
                               )::Union{Hyperrectangle{N}, EmptySet{N}}
                                where {N<:Real}

Overapproximate a convex set by a tight hyperrectangle centered in the origin.

### Input

- `S` -- convex set

### Output

A tight hyperrectangle centered in the origin.

### Algorithm

The center of the box is the origin, and the radius is obtained by computing the
maximum value of the support function evaluated at the canonical directions.
"""
function box_approximation_symmetric(S::LazySet{N};
                                    )::Union{Hyperrectangle{N},
                                             EmptySet{N}} where {N<:Real}
    (c, r) = box_approximation_helper(S)
    if r[1] < 0
        return EmptySet{N}()
    end
    return Hyperrectangle(zeros(N, length(c)), abs.(c) .+ r)
end

# special case: empty set
box_approximation_symmetric(∅::EmptySet) = ∅

"""
    symmetric_interval_hull

Alias for `box_approximation_symmetric`.
"""
symmetric_interval_hull = box_approximation_symmetric

"""
    box_approximation_helper(S::LazySet{N};
                            ) where {N<:Real}

Common code of `box_approximation` and `box_approximation_symmetric`.

### Input

- `S` -- convex set

### Output

A tuple containing the data that is needed to construct a tightly
overapproximating hyperrectangle.

- `c` -- center
- `r` -- radius

### Algorithm

The center of the hyperrectangle is obtained by averaging the support function
of the given convex set in the canonical directions.
The lengths of the sides can be recovered from the distance among support
functions in the same directions.
"""
@inline function box_approximation_helper(S::LazySet{N};
                                         ) where {N<:Real}
    zero_N = zero(N)
    one_N = one(N)
    n = dim(S)
    c = Vector{N}(undef, n)
    r = Vector{N}(undef, n)
    d = zeros(N, n)
    @inbounds for i in 1:n
        d[i] = one_N
        htop = ρ(d, S)
        d[i] = -one_N
        hbottom = -ρ(d, S)
        d[i] = zero_N
        c[i] = (htop + hbottom) / 2
        r[i] = (htop - hbottom) / 2
        if r[i] < 0
            # contradicting bounds => set is empty
            # terminate with first radius entry being negative
            r[1] = r[i]
            break
        end
    end
    return c, r
end

"""
    ballinf_approximation(S::LazySet{N};
                         )::BallInf{N} where {N<:Real}

Overapproximate a convex set by a tight ball in the infinity norm.

### Input

- `S` -- convex set

### Output

A tight ball in the infinity norm.

### Algorithm

The center and radius of the box are obtained by evaluating the support function
of the given convex set along the canonical directions.
"""
function ballinf_approximation(S::LazySet{N};
                              )::Union{BallInf{N}, EmptySet{N}} where {N<:Real}
    zero_N = zero(N)
    one_N = one(N)
    n = dim(S)
    c = Vector{N}(undef, n)
    r = zero_N
    d = zeros(N, n)

    @inbounds for i in 1:n
        d[i] = one_N
        htop = ρ(d, S)
        d[i] = -one_N
        hbottom = -ρ(d, S)
        d[i] = zero_N
        c[i] = (htop + hbottom) / 2
        rcur = (htop - hbottom) / 2
        if (rcur > r)
            r = rcur
        elseif rcur < 0
            # contradicting bounds => set is empty
            return EmptySet{N}()
        end
    end
    return BallInf(c, r)
end

# special case: empty set
ballinf_approximation(∅::EmptySet) = ∅


text:
disthisttrans(xy, step, N, window, pcf=false) =
    disthisttransf(xy, (_,_) -> 1.0, step, N, window, pcf)
disthisttransf(xy, f, step, N, window, pcf=false) =
    disthisttrans12f(xy, xy, f, step, N, window, pcf)
disthisttrans12(xy1, xy2, step, N, window, pcf=false) =
    disthisttrans12f(xy1, xy2, (_,_) -> 1.0, step, N, window, pcf)
function disthisttrans12f(xy1, xy2, f, step, N, window, pcf=false)
    hist = zeros(N)
    if length(xy1) == 0 || length(xy2) == 0
        return hist
    end
    Rmax = (N-1)*step
    T = promote_type(typeof(getx(xy1[1])), typeof(step))
    R2max = T(Rmax^2)
    mult = xy1 === xy2 ? 2.0 : 1.0

    for i1 in 1:length(xy1)
        i2start = xy1 === xy2 ? i1+1 : 1
        for i2 in i2start:length(xy2)
            #dxy = xy1[i1] .- xy2[i2]
            dx = getx(xy1[i1]) - getx(xy2[i2])
            dy = gety(xy1[i1]) - gety(xy2[i2])
            dd = T(dx^2 + dy^2)
            # dd = d2(dxy)
            if dd <= R2max
                d = sqrt(dd)/step
                k = ceil(Int, d)+1
                frac = k - 1 - d
                @assert k <= N
                val = mult*edgetrans(window, (dx, dy))*f(xy1[i1], xy2[i2])
                if pcf && k > 1
                    hist[k-1] += frac*val
                    hist[k] += (1-frac)*val
                else
                    hist[k] += val
                end
            end
        end
    end
    hist
end


text:
module AlgebraicDistance

using LinearAlgebra
using Random

function sumzero!(X)
    for j in axes(X, 2)
        acc = 0
        @simd for i in axes(X, 1)
            acc += X[i, j]
        end
        acc /= size(X, 1)
        @simd for i in axes(X, 1)
            X[i, j] -= acc
        end
    end
    return X
end

gauss_seidel(L) = LowerTriangular(L)
jacobi(L) = Diagonal(diag(L))

function algdist!(X, R, L, P; numiters=10)
    sumzero!(X)
    for i in 1:numiters
        mul!(R, L, X)
        ldiv!(P, R)
        X += R
        sumzero!(X)
    end
    return X
end

function algdist(B, f, p; numvecs=1, numiters=10, rng=Xoshiro())
    L = B * B'
    P = f(L)
    X = randn(rng, size(B, 1), numvecs)
    R = similar(X)
    algdist!(X, R, L, P; numiters=numiters)
    grads = B' * X
    dists = p == Inf ? maximum(abs, grads; dims=2) : sum(x -> abs(x) ^ p, grads; dims=2) .^ (1 / p)
    return X, reshape(dists, :)
end

export algdist, gauss_seidel, jacobi

end # module

max_stars_repo_name: cameton/AlgebraicDistance

max_issues_repo_name: cameton/AlgebraicDistance

max_forks_repo_name: cameton/AlgebraicDistance


text:
using FilePathsBase: FileBuffer

@testset "FileBuffer Tests" begin
    @testset "read" begin
        p = p"../README.md"
        io = FileBuffer(p)
        try
            @test isreadable(io)
            @test !iswritable(io)
            @test !eof(io)
            @test position(io) == 0
            @test read(p) == read(io)
            @test eof(io)
            @test position(io) == length(read(p))
            seekstart(io)
            @test position(io) == 0
            @test !eof(io)
            @test read(p, String) == read(io, String)
            @test eof(io)
        finally
            close(io)
        end

        io = FileBuffer(p)
        try
            for b in read(p)
                @test read(io, UInt8) == b
            end
            @test eof(io)
        finally
            close(io)
        end

        # issue #126: data on first read
        mktemp(SystemPath) do p, _
            write(p, "testing")
            io = FilePathsBase.FileBuffer(p)
            @test read(io, 4) == UInt8['t', 'e', 's', 't']
        end
    end

    @testset "write" begin
        mktmpdir() do d
            p1 = absolute(p"../README.md")

            cd(d) do
                p2 = p"README.md"
                cp(p1, p2)

                io = FileBuffer(p2; read=true,write=true)
                try
                    @test isreadable(io)
                    @test iswritable(io)
                    @test !eof(io)
                    @test position(io) == 0
                    @test read(p1) == read(io)
                    @test eof(io)
                    write(io, "\nHello")
                    write(io, " World!\n")
                    flush(io)

                    @test position(io) == length(read(p1)) + 14
                    txt1 = read(p1, String)
                    txt2 = read(p2, String)
                    @test txt1 != txt2
                    @test occursin(txt1, txt2)
                    @test occursin("Hello World!", txt2)
                finally
                    close(io)
                end

                rm(p2)

                io = FileBuffer(p2; read=true,write=true)
                try
                    write(io, read(p1))
                    flush(io)

                    @test position(io) == length(read(p1))
                    seekstart(io)
                    @test position(io) == 0
                    for b in read(p1)
                        write(io, b)
                    end
                    flush(io)
                    @test read(p1) == read(p2)
                finally
                    close(io)
                end
            end
        end
    end

    @testset "eof" begin
        mktemp(SystemPath) do p, _
            io = FileBuffer(p; write=true)
            @test eof(io)
            write(io, "Hey")
            flush(io)
            @test eof(io)
            seekstart(io)
            @test !eof(io)
            read(io)
            @test eof(io)
        end
    end

    @testset "Custom Types" begin
        jlso = JLSOFile(:msg => "Hello World!")
        mktmpdir() do d
            cd(d) do
                write(p"hello_world.jlso", jlso)
                new_jlso = read(p"hello_world.jlso", JLSOFile)
                @test new_jlso[:msg] == "Hello World!"

                rm(p"hello_world.jlso")
                data = IOBuffer()
                write(data, jlso)
                open(p"hello_world.jlso", "w") do io
                    for x in take!(data)
                        write(io, x)
                    end
                end

                open(p"hello_world.jlso") do io
                    data = UInt8[]
                    push!(data, read(io, UInt8))

                    while !eof(io)
                        push!(data, read(io, UInt8))
                    end

                    new_jlso = read(IOBuffer(data), JLSOFile)
                    @test new_jlso[:msg] == "Hello World!"
                end
            end
        end
    end
end


text:
import Test: @testset, @test

using ScottishTaxBenefitModel

using .FRSHouseholdGetter:
    get_household,
    get_num_households,
    initialise

using .ModelHousehold:
    Household,
    Person

using .Weighting:
    DEFAULT_TARGETS,
    generate_weights,
    initialise_target_dataframe,
    make_target_dataset
using .ExampleHelpers

using CSV
using Tables

start_year=2015

println( "num_households=$(get_num_households())")

#if get_num_households() == 0
@time nhh,total_num_people,nhh2 = initialise( Settings() )

#end

function lsum(v::Vector)
    reduce(+,v)
end

@testset "weighting tests" begin
    nhh = get_num_households()
    hhlds_in_popn = lsum( DEFAULT_TARGETS[42:48]) # sum of all hhld types
    data = make_target_dataset( nhh )
    nr,nc = size(data)

    @test nr == nhh
    @test nc == size( DEFAULT_TARGETS )[1]

    initial_weights = ones(nhh)*hhlds_in_popn/nr

    initial_weighted_popn = (initial_weights' * data)'

    println( "initial-weighted_popn vs targets" )
    println( "target\tinitial\t%diff" )
    for c in 1:nc
        diffpc = 100*(initial_weighted_popn[c]-DEFAULT_TARGETS[c])/DEFAULT_TARGETS[c]
        println( "$c $(DEFAULT_TARGETS[c])\t$(initial_weighted_popn[c])\t$diffpc%")
    end

    open( "scotmat.csv", "w" ) do io
         CSV.write(io, Tables.table(data))
    end

    for c in 1:nc
        @test( sum(data[:,c]) > 0 )
    end

    @time weights = generate_weights(nhh)

    weighted_popn = (weights' * data)'
    @test weighted_popn ≈ DEFAULT_TARGETS
    for c in 1:nc
        diffpc = 100*(weighted_popn[c]-DEFAULT_TARGETS[c])/DEFAULT_TARGETS[c]
        println( "$c $(DEFAULT_TARGETS[c])\t$(weighted_popn[c])\t$diffpc%")
    end

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const DOWNLOAD_HOOKS = Callable[]

function download_url(url::AbstractString)
    for hook in DOWNLOAD_HOOKS
        url = String(hook(url)::AbstractString)
    end
    return url
end

Downloads() = require(PkgId(
        UUID((0xf43a241f_c20a_4ad4, 0x852c_f6b1247861c6)),
        "Downloads",
    ))

"""
    download(url::AbstractString, [path::AbstractString = tempname()]) -> path

Download a file from the given url, saving it to the location `path`, or if not
specified, a temporary path. Returns the path of the downloaded file.

!!! note
    Since Julia 1.6, this function is deprecated and is just a thin wrapper
    around `Downloads.download`. In new code, you should use that function
    directly instead of calling this.
"""
function download(url::AbstractString, path::AbstractString)
    depwarn("Base.download is deprecated; use Downloads.download instead", :download)
    invokelatest(Downloads().download, download_url(url), path)
end
function download(url::AbstractString)
    depwarn("Base.download is deprecated; use Downloads.download instead", :download)
    invokelatest(Downloads().download, download_url(url))
end

max_stars_repo_name: shikhargoswami/julia

max_issues_repo_name: PallHaraldsson/julia

max_forks_repo_name: PallHaraldsson/julia


text:
function asap2(style = Rectangular(2,1))
    arch = TopLevel{2}("asap2")
    arch.metadata["style"] = style

    processor = build_processor_tile(style)
    for r in 0:11, c in 0:12
        add_child(arch, processor, Address(r,c))
    end
    # Add the two processors on the bottom that are not memory processors.
    # Do not include "viterbi", "fft", or "motion estimation" for now.
    add_child(arch, processor, Address(12, 4))
    add_child(arch, processor, Address(12, 9))

    # Add the memory processors.
    memory_processor = build_processor_tile(style, include_memory = true)
    for r in 12, c in (2, 3, 5, 6, 7, 8)
        add_child(arch, memory_processor, Address(r, c))
    end

    # Add the memories. Use left corner as address.
    memory = build_memory(2)
    for r in 13, c in (2, 5, 7)
        add_child(arch, memory, Address(r, c))
    end

    input_handler = build_input_handler(12)
    add_child(arch, input_handler, Address(0,-1))

    output_handler = build_output_handler(12)
    add_child(arch, output_handler, Address(0,13))

    connect_processors(arch, style)
    connect_memories(arch, style)
    connect_io_asap2(arch)

    return arch
end

function connect_io_asap2(arch)

    # Build metadata dictionary for capacity and cost
    metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "link_class"
    )

    # Connect input manually
    for i in 0:11
        src_address = Address(0,-1)
        dst_address = Address(i,0)
        src_name = getname(arch, src_address)
        dst_name = getname(arch, dst_address)

        src_port = Path{Port}(src_name, "out[$i]")
        dst_port = Path{Port}(dst_name, "west_in[0]")
        add_link(arch, src_port, dst_port, metadata = metadata)
    end

    # Connect output manually
    for i in 0:11
        src_address = Address(i,12)
        dst_address = Address(0,13)
        src_name = getname(arch, src_address)
        dst_name = getname(arch, dst_address)

        src_port = Path{Port}(src_name, "east_out[0]")
        dst_port = Path{Port}(dst_name, "in[$i]")
        add_link(arch, src_port, dst_port, metadata = metadata)
    end
end


text:
using Test
using ApproximateComputations

using ImportAll
@importall Base

@testset "Approximation Type" begin
## Approximation Tests
 a = Approximation(6.0)
 b = 4.0

 @test typeof(a - b) == Approximation{typeof(b)}
 @test typeof(a + b) == Approximation{typeof(b)}
 @test typeof(a * b) == Approximation{typeof(b)}
 @test typeof(a / b) == Approximation{typeof(b)}

 @test typeof(Get(a)) == Float64

end

@testset "AST Approximation" begin
	little(x)  = (x * 2) + 5

	@test UpdateEnvironmentForFunction(little) == nothing

	littletree = little(Variable(123))

	@test EmulateTree(littletree) == 251

	ReplaceSubTree(littletree, Variable(0), 1)

	@test EmulateTree(littletree) == 5
end

@testset "Loop Perforation" begin
	expr =	quote
				function newfunc()
					aa = 0
					for i in 1:10
						aa = aa + 1
					end
					aa
				end
			end

	LoopPerforation(expr, UnitRange, ClipFrontAndBack)
	eval(expr)
	@test newfunc() == 8

end

@testset "Memoisation" begin

	# Test custom overwrite hashing
	sinhash(fn, val) = hash(val)
	memoDict = Dict()
	ApproximateHashingMemoise(sin, memoDict, sinhash, 0.4)
	ApproximateHashingMemoise(sin, memoDict, sinhash, 0.1)

	@test ApproximateHashingMemoise(sin, memoDict, sinhash, 0.4) == 0.3894183423086505
	@test ApproximateHashingMemoise(sin, memoDict, sinhash, 0.1) == 0.09983341664682815

	# Test trending batched hash values
	trendingArray = []
	for i in 1:5000
	   push!(trendingArray,[0,0.0,0.0])
	end

	trendingsinhash(fn, val) = 1+Int64.(round(val*10.0))

	@test TrendingMemoisation(sin, trendingsinhash, trendingArray, 3, 0.5)  == 0.479425538604203
	@test TrendingMemoisation(sin, trendingsinhash, trendingArray, 3, 0.55) == 0.5226872289306592
	@test TrendingMemoisation(sin, trendingsinhash, trendingArray, 3, 0.51) == 0.48817724688290753
	@test TrendingMemoisation(sin, trendingsinhash, trendingArray, 3, 0.59) == 0.5563610229127838
	@test TrendingMemoisation(sin, trendingsinhash, trendingArray, 3, 0.59) == 0.5563610229127838
end


text:
@testset "send_extra.jl" begin
    @testset "If already running, timer stays running when new segment sent" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abc")
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        tick!(sender, rto - 5)
        expect_no_seg(sender)
        write_bytes!(sender, "def")
        expect_seg(sender, data="def", payload_size=3)
        tick!(sender, 6)
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        expect_no_seg(sender)
    end

    @testset "Retransmission still happens when expiration time not hit exactly" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abc")
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        tick!(sender, rto - 5)
        expect_no_seg(sender)
        write_bytes!(sender, "def")
        expect_seg(sender, data="def", payload_size=3)
        tick!(sender, 200)
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        expect_no_seg(sender)
    end

    @testset "Timer restarts on ACK of new data" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abc")
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        tick!(sender, rto - 1)
        expect_no_seg(sender)
        write_bytes!(sender, "def")
        expect_seg(sender, data="def", payload_size=3)
        ack_received_test(sender, WrappingInt32(isn + 4), 1000)
        tick!(sender, rto - 1)
        expect_no_seg(sender)
        tick!(sender, 2)
        expect_seg(sender, data="def", payload_size=3, seqno=WrappingInt32(isn + 4))
    end

    @testset "Timer doesn't restart without ACK of new data" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abc")
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        tick!(sender, rto - 5)
        expect_no_seg(sender)
        write_bytes!(sender, "def")
        expect_seg(sender, data="def", payload_size=3, seqno=WrappingInt32(isn + 4))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        tick!(sender, 6)
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        expect_no_seg(sender)
        tick!(sender, 2rto - 5)
        expect_no_seg(sender)
        tick!(sender, 8)
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        expect_no_seg(sender)
    end

    @testset "RTO resets on ACK of new data" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abc")
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        tick!(sender, rto - 5)
        write_bytes!(sender, "def")
        expect_seg(sender, data="def", payload_size=3, seqno=WrappingInt32(isn + 4))
        write_bytes!(sender, "ghi")
        expect_seg(sender, data="ghi", payload_size=3, seqno=WrappingInt32(isn + 7))

        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        tick!(sender, 6)
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        expect_no_seg(sender)
        tick!(sender, 2rto - 5)
        expect_no_seg(sender)
        tick!(sender, 5)
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        expect_no_seg(sender)

        tick!(sender, 4rto - 5)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 4), 1000)
        tick!(sender, rto - 1)
        expect_no_seg(sender)
        tick!(sender, 2)
        expect_seg(sender, data="def", payload_size=3, seqno=WrappingInt32(isn + 4))
        expect_no_seg(sender)
    end

    @testset "fill_window() correctly fills a big window" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))
        MAX_PAYLOAD_SIZE = 1460
        DEFAULT_CAPACITY = 64000
        nicechars = "abcdefghijklmnopqrstuvwxyz"
        bigstring = join(rand(nicechars, DEFAULT_CAPACITY))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        write_bytes!(sender, bigstring)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), window_size)
        expect_state(sender, JLSponge.SYN_ACKED)
        i = UInt(0)
        while i + MAX_PAYLOAD_SIZE < min(length(bigstring), window_size)
            expected_size = min(MAX_PAYLOAD_SIZE, min(length(bigstring), window_size) - i)
            expect_seg(sender, no_flag=true,
                payload_size=expected_size,
                data=@view(bigstring[i+1:i+expected_size]),
                seqno=WrappingInt32(isn + 1 + i))
            i += MAX_PAYLOAD_SIZE
        end
    end

    @testset "Retransmit a FIN-containing segment same as any other" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))

        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abc", true)
        expect_seg(sender, data="abc", payload_size=3, fin=true, seqno=WrappingInt32(isn + 1))
        tick!(sender, rto - 1)
        expect_no_seg(sender)
        tick!(sender, 2)
        expect_seg(sender, data="abc", payload_size=3, fin=true, seqno=WrappingInt32(isn + 1))
    end

    @testset "Retransmit a FIN-only segment same as any other" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abc")
        expect_seg(sender, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        sender_close(sender)
        expect_seg(sender, fin=true, payload_size=0, seqno=WrappingInt32(isn + 4))
        tick!(sender, rto - 1)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 4), 1000)
        tick!(sender, rto - 1)
        expect_no_seg(sender)
        tick!(sender, 2)
        expect_seg(sender, fin=true, payload_size=0, seqno=WrappingInt32(isn + 4))
        tick!(sender, 2rto - 5)
        expect_no_seg(sender)
        tick!(sender, 10)
        expect_seg(sender, fin=true, payload_size=0, seqno=WrappingInt32(isn + 4))
    end

    @testset "Don't add FIN if this would make the segment exceed the receiver's window" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        write_bytes!(sender, "abc", true)
        ack_received_test(sender, WrappingInt32(isn + 1), 3)
        expect_state(sender, JLSponge.SYN_ACKED)

        expect_seg(sender, no_flag=true, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        ack_received_test(sender, WrappingInt32(isn + 2), 2)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 3), 1)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 4), 1)
        expect_seg(sender, fin=true, payload_size=0, seqno=WrappingInt32(isn + 4))

    end

    @testset "Don't send FIN by itself if the window is full" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        write_bytes!(sender, "abc")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1), 3)
        expect_state(sender, JLSponge.SYN_ACKED)

        expect_seg(sender, no_flag=true, data="abc", payload_size=3, seqno=WrappingInt32(isn + 1))
        sender_close(sender)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 2), 2)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 3), 1)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 4), 1)
        expect_seg(sender, fin=true, payload_size=0, seqno=WrappingInt32(isn + 4))

    end

    @testset "MAX_PAYLOAD_SIZE limits payload only" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))
        MAX_PAYLOAD_SIZE = 1460
        nicechars = "abcdefghijklmnopqrstuvwxyz"
        bigstring = join(rand(nicechars, MAX_PAYLOAD_SIZE))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        write_bytes!(sender, bigstring, true)
        ack_received_test(sender, WrappingInt32(isn + 1), 40000)
        expect_seg(sender, data=bigstring, fin=true, payload_size=MAX_PAYLOAD_SIZE, seqno=WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.FIN_SENT)
        ack_received_test(sender, WrappingInt32(isn + 2 + MAX_PAYLOAD_SIZE))
        expect_state(sender, JLSponge.FIN_ACKED)

    end

    @testset "When filling window, treat a '0' window size as equal to '1' but don't back off RTO" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))

        write_bytes!(sender, "abc")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1), 0)
        expect_state(sender, JLSponge.SYN_ACKED)
        expect_seg(sender, data="a", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 1))
        sender_close(sender)
        expect_no_seg(sender)

        for i in 1:5
            tick!(sender, rto - 1)
            expect_no_seg(sender)
            tick!(sender, 1)
            expect_seg(sender, data="a", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 1))
        end
        ack_received_test(sender, WrappingInt32(isn + 2), 0)
        expect_seg(sender, data="b", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 2))

        for i in 1:5
            tick!(sender, rto - 1)
            expect_no_seg(sender)
            tick!(sender, 1)
            expect_seg(sender, data="b", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 2))
        end

        ack_received_test(sender, WrappingInt32(isn + 3), 0)
        expect_seg(sender, data="c", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 3))
        for i in 1:5
            tick!(sender, rto - 1)
            expect_no_seg(sender)
            tick!(sender, 1)
            expect_seg(sender, data="c", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 3))
        end

        ack_received_test(sender, WrappingInt32(isn + 4), 0)
        expect_seg(sender, data="", fin=true, payload_size=0, seqno=WrappingInt32(isn + 4))
        for i in 1:5
            tick!(sender, rto - 1)
            expect_no_seg(sender)
            tick!(sender, 1)
            expect_seg(sender, data="", fin=true, payload_size=0, seqno=WrappingInt32(isn + 4))
        end
    end

    @testset "Unlike a zero-size window, a full window of nonzero size should be respected" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))

        write_bytes!(sender, "abc")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1), 1)
        expect_state(sender, JLSponge.SYN_ACKED)
        expect_seg(sender, data="a", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 1))

        tick!(sender, rto - 1)
        expect_no_seg(sender)
        tick!(sender, 1)
        expect_seg(sender, data="a", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 1))

        sender_close(sender)
        tick!(sender, 2rto - 1)
        expect_no_seg(sender)
        tick!(sender, 1)
        expect_seg(sender, data="a", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 1))

        tick!(sender, 4rto - 1)
        expect_no_seg(sender)
        tick!(sender, 1)
        expect_seg(sender, data="a", no_flag=true, payload_size=1, seqno=WrappingInt32(isn + 1))

        ack_received_test(sender, WrappingInt32(isn + 2), 3)
        expect_seg(sender, data="bc", fin=true, payload_size=2, seqno=WrappingInt32(isn + 2))
    end

    @testset "Repeated ACKs and outdated ACKs are harmless" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        expect_state(sender, JLSponge.SYN_ACKED)

        write_bytes!(sender, "abcdefg")
        expect_seg(sender, data="abcdefg", no_flag=true, payload_size=7, seqno=WrappingInt32(isn + 1))
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        expect_no_seg(sender)
        write_bytes!(sender, "ijkl", true)
        expect_seg(sender, data="ijkl", fin=true, payload_size=4, seqno=WrappingInt32(isn + 8))
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        ack_received_test(sender, WrappingInt32(isn + 8), 1000)
        ack_received_test(sender, WrappingInt32(isn + 12), 1000)
        ack_received_test(sender, WrappingInt32(isn + 12), 1000)
        ack_received_test(sender, WrappingInt32(isn + 12), 1000)
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        tick!(sender, 5rto)
        expect_seg(sender, data="ijkl", fin=true, payload_size=4, seqno=WrappingInt32(isn + 8))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 13), 1000)
        ack_received_test(sender, WrappingInt32(isn + 1), 1000)
        tick!(sender, 5rto)
        expect_no_seg(sender)

        expect_state(sender, JLSponge.FIN_ACKED)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using Test
using IntervalArithmetic
using TaylorModels
# using Revise
using ImplicitDomainQuadrature

set_variables(Float64, "x", order = 2, numvars = 2)
box = IntervalBox(-1..1,2)
x0 = IntervalBox(mid(box))
tm = TaylorModelN(1,1,x0,box)
a = 5.0
b = 2.0
@test muladd(tm,a,b) == a*tm + b
@test muladd(a,tm,b) == a*tm + b
@test muladd(tm,tm,b) == tm*tm + b
@test muladd(a,b,tm) == a*b+tm
tm1 = TaylorModelN(1,1,x0,box)
tm2 = TaylorModelN(1,1,x0,box)
tm3 = TaylorModelN(1,1,x0,box)
@test muladd(tm1,a,tm2) == a*tm1 + tm2
@test muladd(tm1,tm2,tm3) == tm1*tm2 + tm3

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -37.92418043260588 <= q <= 19.457860872935036)

                     add_NL_constraint(m, :(gelu(-0.5948345021955626 + 0.7909419235855366*gelu(0.16572281076834727 + 0.7701313689304001*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.9350053029495258*gelu(0.6292157071886879 + 0.7447123643459479*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.6696024550400765*gelu(-0.5341808090129341 + 0.4032770914369115*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.9733717899660395*gelu(0.43846326345255404 + -0.02862602361460853*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.07261217388131413*gelu(0.02496868545483233 + -0.04485855107115411*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + gelu(-0.7044745597894595 + -0.85619817917634*gelu(0.16572281076834727 + 0.7701313689304001*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.37830231556965854*gelu(0.6292157071886879 + 0.7447123643459479*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.12156626196118481*gelu(-0.5341808090129341 + 0.4032770914369115*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.3744591738356289*gelu(0.43846326345255404 + -0.02862602361460853*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.4985780257939494*gelu(0.02496868545483233 + -0.04485855107115411*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + gelu(-0.9683275825894269 + -0.7637278723892034*gelu(0.16572281076834727 + 0.7701313689304001*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.9143661573707975*gelu(0.6292157071886879 + 0.7447123643459479*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.7824888723979511*gelu(-0.5341808090129341 + 0.4032770914369115*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.6866950641423624*gelu(0.43846326345255404 + -0.02862602361460853*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.09633196432882363*gelu(0.02496868545483233 + -0.04485855107115411*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + gelu(-0.5143076265809325 + -0.11312975790997548*gelu(0.16572281076834727 + 0.7701313689304001*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7215302240627146*gelu(0.6292157071886879 + 0.7447123643459479*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.8436266175918532*gelu(-0.5341808090129341 + 0.4032770914369115*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.8029468027318871*gelu(0.43846326345255404 + -0.02862602361460853*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.07363029133667531*gelu(0.02496868545483233 + -0.04485855107115411*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + gelu(-0.8626839052591899 + -0.6726505975051698*gelu(0.16572281076834727 + 0.7701313689304001*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.7431517048783043*gelu(0.6292157071886879 + 0.7447123643459479*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.03587478002549904*gelu(-0.5341808090129341 + 0.4032770914369115*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7126331133171049*gelu(0.43846326345255404 + -0.02862602361460853*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7945114361960401*gelu(0.02496868545483233 + -0.04485855107115411*gelu(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*gelu(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*gelu(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*gelu(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*gelu(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
" Implementation of a separating axis test for the intersection of two convex shapes.
# Arguments
- `shape1::Dict` - Convex shape with entry `vertices`
- `shape2::Dict` - Convex shape with entry `vertices`
"
function do_shapes_intersect(shape1::Dict, shape2::Dict; axes=nothing)
    int_result = true
    min_distances = [-1.]

    if isnothing(axes)
        axes = get_axes([shape1, shape2])
    end

    num_axes = size(axes)[2]
    for i=1:num_axes
        axis = axes[:, i:i]
        p1 = project(shape1, axis)
        p2 = project(shape2, axis)
        overlap = (p2[1] <= p1[1] <= p2[2] || p2[1] <= p1[2] <= p2[2] ||
                   p1[1] <= p2[1] <= p1[2] || p1[1] <= p2[2] <= p1[2] )
        push!(min_distances, minimum([abs(p1[2] - p2[1]), abs(p2[2] - p1[1])]))
        if !overlap
            # push!(min_distances, minimum([abs(p1[2] - p2[1]), abs(p2[2] - p1[1])]))
            int_result = false
        end
    end

    return int_result, maximum(min_distances)
end

" Checks if shape2 contains shape1.
# Arguments
- `shape1::Dict` - Convex shape with entry `vertices`
- `shape2::Dict` - Convex shape with entry `vertices`
"
function does_shape2_contain_shape1(shape1::Dict, shape2::Dict; axes=nothing)
    int_result = false
    min_distance = Inf

    if isnothing(axes)
        axes = get_axes([shape1, shape2])
    end

    num_axes = size(axes)[2]
    for i=1:num_axes
        axis = axes[:, i:i]
        p1 = project(shape1, axis)
        p2 = project(shape2, axis)
        # overlap = (p2[1] <= p1[1] <= p2[2] || p2[1] <= p1[2] <= p2[2] ||
        #            p1[1] <= p2[1] <= p1[2] || p1[1] <= p2[2] <= p1[2] )
        contain = p2[1] <= p1[1] <= p2[2] && p2[1] <= p1[2] <= p2[2]
        if contain
            min_distance = minimum([abs(p1[1] - p2[1]), abs(p1[2] - p2[2]), min_distance])
            int_result = true
        else
            int_result = false
            break
        end
    end

    return int_result, min_distance
end

" Projects the shape vertices onto the axis.
# Arguments
- `shape::Dict` - Shape with entry `vertices`
- `axis` - Axis of projection
"
function project(shape::Dict, axis)
    dots = [(vertex*axis)[1] for vertex in shape["vertices"]]
    min = minimum(dots)
    max = maximum(dots)
    return [min, max]
end

" Gets the principal normals of each shape.
# Arguments
- `shapes::Vector` - Vector with shape definitions
"
function get_axes(shapes::Vector)
    axes = []
    for shape in shapes
        num_v = length(shape["vertices"])
        for i=1:num_v
            normals = nullspace(shape["vertices"][i]-(i==num_v ? shape["vertices"][1] : shape["vertices"][i+1]))
            if isempty(axes)
                axes = normals
            else
                # TODO: remove redundant axes
                axes = hcat([axes normals])
            end
        end
    end
    return axes
end