text:
## File contains the fitting funcs with a parametric expansion...

using JuMP
using Gurobi
using Distributions
using Roots

#represent the desired funciton as a combination of features
#featureFuncs is itself a collection of function pointers
# phi : (prices, gdp) -> Real

#Cheapest way possible
function genFeaturesPoly(deg)
	features = Array(Any, binomial(deg + 3, deg))
	indx = 1
	for p1deg = 0:deg
		for p2deg = 0:deg
			for gdpdeg = 0:deg
				if p1deg + p2deg + gdpdeg <= deg
					features[indx] = (p1, p2, gdp) -> p1^p1deg * p2^p2deg * gdp^gdpdeg
					indx+=1
				end
			end
		end
	end
	features
end

function genFeaturesExp( rng )
	features = [(p1, p2, gdp)-> exp(r * p1) for r in rng]
	append!(features, [(p1, p2, gdp) -> exp(r * p2) for r in rng])
	append!(features, [(p1, p2, gdp) -> exp(r * gdp) for r in rng])
	append!(features, [(p1, p2, gdp) -> p1, (p1, p2, gdp) -> p2, (p1, p2, gdp) -> gdp, (p1, p2, gdp)-> 1])
	features
end

createfun(flambdas, features) =
	(prices, gdp) -> sum([flambdas[i] * features[i](prices..., gdp) for i=1:length(flambdas) ])

function setUpFitting(features)
	m = Model(solver=GurobiSolver(OutputFlag=false))
#	m = Model()
	N=length(features)
	@defVar(m, vlambdas1[1:N])
	@defVar(m, vlambdas2[1:N])
	return m, vlambdas1, vlambdas2
end

#Gens the exp features and coefficients imultaneously to constrain increasingness etc.
function setUpFitting( rng::Range{Float64} )
	N = 3*length(rng) + 4 # add linears and the constants
	features = genFeaturesExp(rng)
	m = Model(solver=GurobiSolver(OutputFlag=false))
	@defVar(m, vlambdas1[1:N])
	@defVar(m, vlambdas2[1:N])

	#require that MR(infty) -> -infty
	# for (ix, r) = enumerate(rng)
	# 	if r >= 0
	# 		@addConstraint(m, vlambdas1[ix] <= 0)
	# 		@addConstraint(m, vlambdas2[ix + length(rng)] <= 0)
	# 	end
	# 	# if r >= 0  #increasing functions
	# 	# 	@addConstraint(m, vlambdas1[ix] == 0)  #MR1 should be decreasing in p1
	# 	# 	@addConstraint(m, vlambdas1[ix + length(rng)] >=0 )  #MR1 should be increasing in p2
	# 	# 	@addConstraint(m, vlambdas2[ix] >= 0)  #MR2 should be increasing in p1
	# 	# 	@addConstraint(m, vlambdas2[ix + length(rng)] == 0)  #MR2 should be decreasing in p2
	# 	# elseif r <= 0  #decreasing functions
	# 	# 	@addConstraint(m, vlambdas1[ix] >= 0)   #MR1 should be decreasing in p1
	# 	# 	@addConstraint(m, vlambdas1[ix + length(rng)] == 0)   #MR1 should be increasing in p2
	# 	# 	@addConstraint(m, vlambdas2[ix] == 0)  #MR2 should be increasing in p1
	# 	# 	@addConstraint(m, vlambdas2[ix + length(rng)] >=0)  #MR2 should be decreasing in p2
	# 	# end
	# end

	#same increasing/decreasing for the linear terms
	# @addConstraint(m, vlambdas1[3 *length(rng) + 1] <= 0)
	# @addConstraint(m, vlambdas2[3 *length(rng) + 2] <= 0)
	# @addConstraint(m, vlambdas1[3 *length(rng) + 2] >= 0)
	# @addConstraint(m, vlambdas2[3 *length(rng) + 1] >= 0)

	#disallow the weirdness in the functional form
	for ix = 1:length(rng)
		@addConstraint(m, vlambdas1[ix + 2 *length(rng)] == 0)
		@addConstraint(m, vlambdas2[ix + 2 *length(rng)] == 0)
	end
	return m, vlambdas1, vlambdas2, features
end

#insist that the fitted function is decreasing on the data points for mean values of other stuff
function addDecreasingCnsts(m, vlambdas, features, otherP, meanGDP, price_data, pmax, ID)
	price_data = sort(price_data)
	sort_prices(myP, yourP) = ID == :1 ? (myP, yourP) : (yourP, myP)
	for j = 1:(length(price_data)- 1)
		prices_j = sort_prices(price_data[j], otherP)
		prices_jp = sort_prices(price_data[j+1], otherP)
		@addConstraint(m, sum{features[i](prices_j..., meanGDP) * vlambdas[i], i=1:length(vlambdas)} >=
		 					sum{features[i](prices_jp..., meanGDP) * vlambdas[i], i=1:length(vlambdas)} )
	end

	#add one last one for pmax in case its not achieved...
	prices_j = sort_prices(price_data[length(price_data)], otherP)
	prices_jp = sort_prices(pmax, otherP)
	@addConstraint(m, sum{features[i](prices_j..., meanGDP) * vlambdas[i], i=1:length(vlambdas)} >=
	 					sum{features[i](prices_jp..., meanGDP) * vlambdas[i], i=1:length(vlambdas)} )

end

#insist that the fitted function is positive at zero for mean values
function addPositivity(m, vlambdas, features, otherP, myP, meanGDP, ID; TOL=0)
	if ID == :1
		@addConstraint(m, sum{features[i](myP, otherP, meanGDP) * vlambdas[i], i=1:length(vlambdas)} >= TOL)
	else
		@addConstraint(m, sum{features[i](otherP, myP, meanGDP) * vlambdas[i], i=1:length(vlambdas)} >= TOL)
	end
end

#insist that the MR at the mean point should be equal to given  quantity for each
function normalize(m, vlambdas1, vlambdas2, features, meanPs, meanGDP, MR1, MR2)
	N = length(vlambdas1)
	@addConstraint(m, MR1 == sum{features[i](meanPs..., meanGDP)*vlambdas1[i], i=1:N})
	@addConstraint(m, MR2 == sum{features[i](meanPs..., meanGDP)*vlambdas2[i], i=1:N})
end

function normalize(m, vlambdas, features, prices, error, MR)
	N = length(vlambdas)
	@addConstraint(m, MR == sum{features[i](prices..., error) * vlambdas[i], i=1:N})
end

#adds the constraints to the model for the residual for this observation
function addResids(m, vlambdas1, vlambdas2, prices_t, GDP_t, pmax, features)
	#add some supplementary variables for ease:
	x = [prices_t GDP_t]; N = length(vlambdas1)
	@defVar(m, MR[1:2])
	@addConstraint(m, MR[1] == sum{features[i](prices_t..., GDP_t)*vlambdas1[i], i=1:N})
	@addConstraint(m, MR[2] == sum{features[i](prices_t..., GDP_t)*vlambdas2[i], i=1:N})

	@defVar(m, resid >= 0)
	@defVar(m, y[1:2] >= 0)
	@addConstraint(m, y[1] >= MR[1])
	@addConstraint(m, y[2] >= MR[2])
	@addConstraint(m, pmax * (y[1] + y[2]) - prices_t[1] * MR[1] - prices_t[2] * MR[2] <= resid)
	return resid
end

#calc an out of sample residual.  may be a better way to do this.
function calcResid(flambdas1, flambdas2, prices, GDP, pmax, features)
	N = length(flambdas1); MR = Array(Float64, 2)
	MR[1] = sum([flambdas1[i]*features[i](prices..., GDP) for i = 1:N])
	MR[2] = sum([flambdas2[i]*features[i](prices..., GDP) for i = 1:N])

	y = Array(Float64, 2)
	y[1] = max(0, MR[1])
	y[2] = max(0, MR[2])

	return pmax * (y[1] + y[2]) - prices[1] * MR[1] - prices[2] * MR[2]
end

function bestResponsePrice(flambdas, features, prices, error, pmax, ID)
	fMR = createfun(flambdas, features)
	if ID == :1
		# pstar = newton(p-> fMR([p, prices[2]], error), prices[1])
		pstar = fzero(p-> fMR([p, prices[2]], error), prices[1])
		#pstar = customRoot(p-> fMR([p, prices[2]], error), pmax)  #works really well with exp features
	else
		# pstar = newton(p-> fMR([prices[1], p], error), prices[2])
		pstar = fzero(p-> fMR([prices[1], p], error), prices[2])
		# pstar = customRoot(p-> fMR([prices[1], p], error), pmax)   #works really well with exp features
	end
	# assert( pstar >= 0 )
	pstar = max(0, pstar)
	return min(pstar, pmax)
end

#assumes f(0) > 0 and we have decreasingness
function customRoot(f, pmax; MAX_BRACKET = 100)
	assert( f(1e-3) > 0 )  # this should be because of the positivity issue
	if f(pmax) > 0
		return pmax
	end
	return fzero(f, [1e-3, pmax])
end

#solve for equilibrium prices.  True model is log with a random, firm dependent effect
function solveNashPrices(flambdas1, flambdas2, features, error, pmax; MAX_ITERS = 20, TOL=1e-6, trace=false,
		prices=fill(convert(Float64, pmax / 2), 2))
	dist = 0
	for i = 1:MAX_ITERS
		prices_old = copy(prices)
		prices[1] = bestResponsePrice(flambdas1, features, prices_old, error, pmax, :1)
		prices[2] = bestResponsePrice(flambdas2, features, prices_old, error, pmax, :2)
		dist = norm(prices - prices_old)
		if trace
			println([prices; dist]')
		end
		if dist < TOL
			return prices
		end
	end
	println("Warning: Max Iters Reached \t $dist")
	return prices
end

ext: jl

max_stars_repo_name: nlperic/InverseVIs

max_issues_repo_name: nlperic/InverseVIs


text:
"""
    SpectralNeuralOperator{F,L,F1,F2,P}

A type representing a Neural Operator whose forward pass is of the form

    y(t) = σ((B*x)(t) - λ(t)*(B*x)(t) + b(t))

where `x`, `λ` and `b` are functions, and `B` is a linear operator.

By default, `y` is projected in the Chebyshev polynomial basis before outputing. This is done to avoid computing multiple nested integrals when stacking multiple layers.

## Fields
 * `B::L` A linear operator on functions ``[0,1] → R``, generally a kernel operator.
 * `λ::F1` a function ``[0,1] → R``
 * `b::F2` a function ``[0,1] → R`` used as bias
 * `σ::F` the activation function
 * `project::P` a projection operator used on the output of the layer. Since by default, the output is a lazy function, this is meant to
"""
struct SpectralNeuralOperator{F,L,F1,F2,P}
    B::L # Linear Operator
    λ::F1 # Spectral Scaling Function
    b::F2 # Bias function
    σ::F # activation function
    project::P # Projection operator
end

SpectralNeuralOperator(B, λ, b) = SpectralNeuralOperator(B, λ, b, identity, f -> Fun(f, 0..1))

function Base.show(io::IO, ::MIME"text/plain", x::SpectralNeuralOperator{F,L,F1,F2,P}) where {F,L,F1,F2,P}
    println(io, "SpectralNeuralOperator with\n kernel type $L\n scaling type $F1\n bias type $F2\n activation function $F\n projection operator $P")
end

Flux.@functor SpectralNeuralOperator

function kernel_fun(sno::SpectralNeuralOperator{F,<:KernelOperator,F1,F2,P}) where {F1,F2,F,P}
    return function k_fun(t,s)
        (1 - sno.λ(t))*sno.B.kernel(t,s)
    end
end

function (sno::SpectralNeuralOperator)(x)
    B, λ, b, σ = sno.B, sno.λ, sno.b, sno.σ
    y = B(x)
    z = t -> σ((1- λ(t))*y(t) + b(t))
    sno.project(z)
end

function b_pullback(sno::SpectralNeuralOperator{F,L,F1,F2,P}, δz, y) where {F,L,F1,F2<:ParametricFun,P}
    dσ(u) = ForwardDiff.derivative(sno.σ, u)
    function fun(u)
        dσ(y(u)) .* to_grad(sno.b, u)
    end
    δb = product_integral(fun, δz)
    return from_grad(sno.b, δb)
end

# default behaviour, no tangent
b_pullback(sno, δz, δy) = ZeroTangent()

function λ_pullback(sno::SpectralNeuralOperator{F,L,F1,F2,P}, δz, y) where {F,L,F1<:ParametricFun,F2,P}
    dσ(u) = ForwardDiff.derivative(sno.σ, u)
    function fun(u)
        dσ(y(u)) .* to_grad(sno.λ, u)
    end
    δλ =  kernel_double_integral(kernel_fun(sno), fun, δz)
    return from_grad(sno.λ, δλ)
end

λ_pullback(sno, δz, δy) = ZeroTangent()

function x_pullback(sno::SpectralNeuralOperator, δz, y)
    dσ(u) = ForwardDiff.derivative(sno.σ, u)
    function δx(s)
        dσ(y(s)) * kernel_transpose_integral(kernel_fun(sno), δz, s)
    end
    return sno.project(δx)
end

function ChainRulesCore.rrule(sno::SpectralNeuralOperator, x)
    Bx = sno.B(x)
    y = sno.project(t -> (1 - sno.λ(t))*Bx(t) + sno.b(t))
    z = t -> sno.σ( (1 - sno.λ(t)) * y(t) + sno.b(t) )
    function sno_pullback(δz)
        δB = ZeroTangent()
        δλ = λ_pullback(sno, δz, y)
        δb = b_pullback(sno, δz, y)
        δσ = NoTangent()
        δproject = NoTangent()
        δx = x_pullback(sno, δz, y)
        return Tangent{SpectralNeuralOperator}(;B=δB, λ=δλ, b=δb, σ=δσ, project=δproject), δx
    end
    return sno.project(z), sno_pullback
end


text:
struct NotSet end
struct CASing end

Reagents.Ref{T}() where {T} = GenericRef{T}()
Reagents.Ref{T}(x) where {T} = GenericRef{T}(x)
Reagents.Ref(x::T) where {T} = GenericRef{T}(x)

mutable struct GenericRef{T} <: Reagents.Ref{T}
    @atomic value::Any
end

GenericRef{T}() where {T} = GenericRef{T}(NotSet())


text:
x="println(\"x=\$(repr(x))\\n\$x\")"
println("x=$(repr(x))\n$x")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file includes the plugin for the calculation variance of data

using Statistics

"""
    Variance(dims::Int)

Constructs a `Variance` plugin. The [`process(plg::Variance, x)`](@ref) function takes the variance of the input data `x` along dimension `dims`.
"""
struct Variance <: AbstractPlugin
    dims::Int
end
Variance(;dims::Int=1) = Variance(dims)

show(io::IO, plg::Variance) = print(io, "Mean(dims:$(plg.dims))")

"""
    process(plg::Std, x)

Returns the standard deviation of `x` along the dimension `plg.dims`.

# Example
```julia
julia> x = collect(reshape(1:16, 4,4))
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> plg = Plugins.Variance(dims=1)
Mean(dims:1)

julia> process(plg, x)
1×4 Array{Float64,2}:
 1.66667  1.66667  1.66667  1.66667
```
"""
process(plg::Variance, x) = var(x, dims=plg.dims)


text:
@testset "Conversions" begin
    x = 1.25
    i = 4

    # Floats test conversion back and forth
    @test x == Float64(Posit8(x))
    @test x == Float64(Posit16(x))
    @test x == Float64(Posit32(x))

    @test x == Float32(Posit8_2(x))
    @test x == Float32(Posit16_2(x))
    @test x == Float32(Posit24_2(x))

    # Posit to Posit back and forth
    if Sys.isapple() || Sys.iswindows()
        @test_skip Posit8(x) == Posit8(Posit16(x))
        @test_skip Posit32(x) == Posit32(Posit16(Posit8_2(x)))
    else
        @test Posit8(x) == Posit8(Posit16(x))
        @test Posit32(x) == Posit32(Posit16(Posit8_2(x)))
    end

    # Int to posit back and forth
    @test i == Int64(Posit8(i))
    @test i == Int64(Posit16(i))
    @test i == Int64(Posit32(i))

    # Boolean tests
    @testset for T in (Posit8,Posit16,Posit32)
        T(x) == true*T(x)
        zero(T) == false*T(x)
        T(x)+one(T) == true + T(x)
        T(x) == false + T(x)
        T(true) == one(T)
        T(false) == zero(T)
    end
end

@testset "0,±1,±2,±4" begin
	for f in Float32[-4,-2,-1,0,1,2,4]
		@test Posit16(f) == Posit16_old(f)
		@test f == Float32(Posit16(f))
	end
end

@testset "0,±1/2,±1/4,±1/8" begin
	for f in Float32[-1/8,-1/4,-1/2,0,1/2,1/4,1/8]
		@test Posit16(f) == Posit16_old(f)
		@test f == Float32(Posit16(f))
	end
end

@testset "With zero mantissa bits" begin
	for f in Float32[	0.000244140625,
						0.0009765625,
						0.001953125,
						0.00390625,
						0.0078125,
						0.015625,
						0.03125,
						0.0625,
						0.125,
						0.25,
						0.5]
		@test Posit16(f) == Posit16_old(f)
		@test Posit16(-f) == Posit16_old(-f)
		@test Posit16(1/f) == Posit16_old(1/f)
		@test Posit16(-1/f) == Posit16_old(-1/f)

		@test f == Float32(Posit16(f))
		@test -f == Float32(Posit16(-f))
		@test 1/f == Float32(Posit16(1/f))
		@test -1/f == Float32(Posit16(-1/f))
	end
end

@testset "NaN, Inf" begin
	for f in [NaN32,Inf32,-Inf32]
		@test Posit16(f) == Posit16_old(f)
		@test isnan(Float32(Posit16(f)))
	end

	for _ in 1:10
		# create various NaNs
		f = reinterpret(Float32,reinterpret(UInt32,NaN32)+(rand(UInt32)>>10))
		@test Posit16(f) == Posit16_old(f)
		@test isnan(Posit16(f))
	end
end

@testset "U(0,1)" begin
	for f in rand(Float32,100)
		@test Posit16(f) == Posit16_old(f)
		@test Float32(Posit16(f)) == Float32_old(Posit16_old(f))

		# idempotence
		f32 = Float32(Posit16(f))
		@test f32 == Float32(Posit16(f32))
	end
end

@testset "U(1,21)" begin
	for f in 1 .+ 20*rand(Float32,100)
		@test Posit16(-f) == Posit16_old(-f)
		@test Posit16(f) == Posit16_old(f)
		@test Float32(Posit16(f)) == Float32_old(Posit16_old(f))

		# idempotence
		f32 = Float32(Posit16(f))
		@test f32 == Float32_old(Posit16_old(f32))
	end
end

@testset "N(0,1)" begin
	for f in randn(Float32,100)
		@test Posit16(f) == Posit16_old(f)
		@test Float32(Posit16(f)) == Float32_old(Posit16_old(f))

		# idempotence
		f32 = Float32(Posit16(f))
		@test f32 == Float32(Posit16(f32))
	end
end

@testset "Overflow and underflow" begin
    for f in Float32[2f8,3f8,4f8,5f8,1f9,1f10,1f15,1f20,1f25,1f30,1f35,floatmax(Float32)/4]
        @test floatmax(Posit16) == Posit16(f)
        @test -floatmax(Posit16) == Posit16(-f)
    end

    @test isnan(Posit16(NaN32))
    @test isnan(Posit16(Inf32))
    @test isnan(Posit16(-Inf32))
end


text:
module ThreedWrapper

using JavaCall

const AbstractRegion = @jimport org.hipparchus.geometry.partitioning.AbstractRegion
const AbstractSubHyperplane = @jimport org.hipparchus.geometry.partitioning.AbstractSubHyperplane
const BSPTree = @jimport org.hipparchus.geometry.partitioning.BSPTree
const BoundaryProjection = @jimport org.hipparchus.geometry.partitioning.BoundaryProjection
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const EnclosingBall = @jimport org.hipparchus.geometry.enclosing.EnclosingBall
const Enum = @jimport java.lang.Enum
const Euclidean2D = @jimport org.hipparchus.geometry.euclidean.twod.Euclidean2D
const Euclidean3D = @jimport org.hipparchus.geometry.euclidean.threed.Euclidean3D
const Field = @jimport org.hipparchus.Field
const FieldLine = @jimport org.hipparchus.geometry.euclidean.threed.FieldLine
const FieldPosition = @jimport java.text.FieldPosition
const FieldRotation = @jimport org.hipparchus.geometry.euclidean.threed.FieldRotation
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const Hyperplane = @jimport org.hipparchus.geometry.partitioning.Hyperplane
const IntervalsSet = @jimport org.hipparchus.geometry.euclidean.oned.IntervalsSet
const JString = @jimport java.lang.String
const JVector = @jimport org.hipparchus.geometry.Vector
const Line = @jimport org.hipparchus.geometry.euclidean.threed.Line
const List = @jimport java.util.List
const Locale = @jimport java.util.Locale
const NumberFormat = @jimport java.text.NumberFormat
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const OutlineExtractor = @jimport org.hipparchus.geometry.euclidean.threed.OutlineExtractor
const ParsePosition = @jimport java.text.ParsePosition
const Plane = @jimport org.hipparchus.geometry.euclidean.threed.Plane
const Point = @jimport org.hipparchus.geometry.Point
const PolyhedronsSet = @jimport org.hipparchus.geometry.euclidean.threed.PolyhedronsSet
const PolyhedronsSet_BRep = @jimport org.hipparchus.geometry.euclidean.threed.PolyhedronsSet$BRep
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const Region = @jimport org.hipparchus.geometry.partitioning.Region
const Region_Location = @jimport org.hipparchus.geometry.partitioning.Region$Location
const Rotation = @jimport org.hipparchus.geometry.euclidean.threed.Rotation
const RotationConvention = @jimport org.hipparchus.geometry.euclidean.threed.RotationConvention
const RotationOrder = @jimport org.hipparchus.geometry.euclidean.threed.RotationOrder
const Segment = @jimport org.hipparchus.geometry.euclidean.threed.Segment
const Space = @jimport org.hipparchus.geometry.Space
const SphereGenerator = @jimport org.hipparchus.geometry.euclidean.threed.SphereGenerator
const SphericalCoordinates = @jimport org.hipparchus.geometry.euclidean.threed.SphericalCoordinates
const StringBuffer = @jimport java.lang.StringBuffer
const SubHyperplane = @jimport org.hipparchus.geometry.partitioning.SubHyperplane
const SubHyperplane_SplitSubHyperplane = @jimport org.hipparchus.geometry.partitioning.SubHyperplane$SplitSubHyperplane
const SubLine = @jimport org.hipparchus.geometry.euclidean.threed.SubLine
const SubPlane = @jimport org.hipparchus.geometry.euclidean.threed.SubPlane
const Transform = @jimport org.hipparchus.geometry.partitioning.Transform
const Vector1D = @jimport org.hipparchus.geometry.euclidean.oned.Vector1D
const Vector2D = @jimport org.hipparchus.geometry.euclidean.twod.Vector2D
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D
const Vector3DFormat = @jimport org.hipparchus.geometry.euclidean.threed.Vector3DFormat
const VectorFormat = @jimport org.hipparchus.geometry.VectorFormat

export Euclidean3D
export FieldLine
export FieldRotation
export FieldVector3D
export Line
export OutlineExtractor
export Plane
export PolyhedronsSet
export Rotation
export RotationConvention
export RotationOrder
export Segment
export SphereGenerator
export SphericalCoordinates
export SubLine
export SubPlane
export Vector3D
export Vector3DFormat
export add
export angle
export apply_inverse_to
export apply_to
export apply_transform
export ball_on_support
export build_new
export check_point
export closest_point
export compare_to
export compose
export compose_inverse
export contains
export copy_self
export cross_product
export describe_constable
export distance
export distance1
export distance_inf
export distance_sq
export dot_product
export empty_hyperplane
export equals
export first_intersection
export format
export from_direction
export get_a1
export get_a2
export get_a3
export get_abscissa
export get_alpha
export get_angle
export get_angles
export get_available_locales
export get_axis
export get_b_rep
export get_barycenter
export get_boundary_size
export get_cartesian
export get_class
export get_declaring_class
export get_delta
export get_dimension
export get_direction
export get_end
export get_format
export get_hyperplane
export get_identity
export get_instance
export get_line
export get_matrix
export get_minus_i
export get_minus_j
export get_minus_k
export get_na_n
export get_negative_infinity
export get_norm
export get_norm1
export get_norm_inf
export get_norm_sq
export get_normal
export get_offset
export get_origin
export get_outline
export get_phi
export get_plus_i
export get_plus_j
export get_plus_k
export get_point_at
export get_positive_infinity
export get_prefix
export get_q0
export get_q1
export get_q2
export get_q3
export get_r
export get_remaining_region
export get_rotation_order
export get_segments
export get_separator
export get_size
export get_space
export get_start
export get_sub_space
export get_suffix
export get_theta
export get_tolerance
export get_tree
export get_u
export get_v
export get_vector3d_format
export get_x
export get_y
export get_z
export get_zero
export hash_code
export intersection
export is_empty
export is_full
export is_infinite
export is_na_n
export is_similar_to
export name
export negate
export normalize
export notify
export notify_all
export ordinal
export orthogonal
export parse
export point_at
export project
export project_to_boundary
export reset
export reunite
export revert
export revert_self
export rotate
export same_orientation_as
export scalar_multiply
export split
export subtract
export to_array
export to_cartesian_gradient
export to_cartesian_hessian
export to_rotation
export to_space
export to_string
export to_sub_space
export to_vector3d
export translate
export value_of
export values
export wait
export whole_hyperplane
export whole_line
export whole_space

include("euclidean3_d.jl")
include("field_line.jl")
include("field_rotation.jl")
include("field_vector3_d.jl")
include("line.jl")
include("outline_extractor.jl")
include("plane.jl")
include("polyhedrons_set.jl")
include("rotation.jl")
include("rotation_convention.jl")
include("rotation_order.jl")
include("segment.jl")
include("sphere_generator.jl")
include("spherical_coordinates.jl")
include("sub_line.jl")
include("sub_plane.jl")
include("vector3_d.jl")
include("vector3_d_format.jl")

end


text:
import VecMathOptInterface
const MOI = VecMathOptInterface
const MOIU = MOI.Utilities

abstract type AbstractDummyModel <: MOI.ModelLike end

function MOI.empty!(::AbstractDummyModel) end

function MOI.copy_to(dest::AbstractDummyModel, src::MOI.ModelLike)
    return MOIU.default_copy_to(dest, src)
end

MOI.supports_incremental_interface(::AbstractDummyModel) = true

MOI.supports(::AbstractDummyModel, ::MOI.ObjectiveSense) = true

function MOI.supports(
    ::AbstractDummyModel,
    ::MOI.ConstraintPrimalStart,
    ::Type{<:MOI.ConstraintIndex},
)
    return true
end
function MOI.supports_constraint(
    ::AbstractDummyModel,
    ::Type{MOI.VariableIndex},
    ::Type{MOI.EqualTo{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::AbstractDummyModel,
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.Zeros},
)
    return true
end

struct DummyModel <: AbstractDummyModel end

# Implements add_variable and add_constraint
struct DummyModelWithAdd <: AbstractDummyModel end
MOI.add_variable(::DummyModelWithAdd) = MOI.VariableIndex(0)
MOI.add_variables(::DummyModelWithAdd, n) = fill(MOI.VariableIndex(0), n)
function MOI.add_constraint(
    ::DummyModelWithAdd,
    ::MOI.VariableIndex,
    ::MOI.EqualTo{Float64},
)
    return MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(0)
end


text:
using KernelAbstractions

@kernel function darcyKernel!(uⁿ⁺¹, vⁿ⁺¹, wⁿ⁺¹, h, Δx, Δy, Δz, K)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        uⁿ⁺¹[i,j,k] = -K[i,j,k] * ∇3dᶜ(h, i, j, k, (Δx, Δy, Δz))[1]    # u
        vⁿ⁺¹[i,j,k] = -K[i,j,k] * ∇3dᶜ(h, i, j, k, (Δx, Δy, Δz))[2]    # v
        wⁿ⁺¹[i,j,k] = -K[i,j,k] * ∇3dᶜ(h, i, j, k, (Δx, Δy, Δz))[3]    # w
    end
end


text:
using Test
using Distances
using Distributions
using TripletEmbeddings
using EuclideanDistanceMatrices

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

@testset "EDMs.jl" begin
    # EDM
    d = 20
    n = 10
    @test_throws AssertionError EDM(d, n)
    @test_throws DomainError EDM(-1, n)
    @test_throws DomainError EDM(d, -1)

    d = 2
    n = 10
    D = pairwise(SqEuclidean(), rand(d, n), dims=2)
    @test_throws DomainError EDM(D + rand(n, n))

    D = Matrix{Float64}(undef, 1, 1)
    D[1,1] = 1
    @test_throws DomainError EDM(D)

    # NoisyEDM
    D = EDM(2,10)
    N = rand(2,2)
    X = Embedding(2,10)

    @test_throws DimensionMismatch NoisyEDM(D, N)
    @test_throws DimensionMismatch NoisyEDM(X, N)

    # NoisyMaskedEDM
    X = Embedding(2,10)
    N = rand(10,10)
    mask = rand(3,4)
    @test_throws DimensionMismatch NoisyMaskedEDM(X, N, mask)

    D = EDM(2,9)
    N = rand(10,10)
    mask = rand(3,4)
    @test_throws DimensionMismatch NoisyMaskedEDM(X, N, mask)

    # isedm
    X = rand(3,10)
    D = pairwise(SqEuclidean(), X, dims=2)
    @test isedm(D) == true
    @test isedm(X) == false

    X = rand(1,2)
    D = pairwise(SqEuclidean(), X, dims=2)
    @test isedm(D) == true

    D = pairwise(SqEuclidean(), X, dims=1)
    @test_throws DomainError isedm(D)
end

@testset "masks.jl" begin
    n = 20
    deletions = 20
    mask = randommask(n, deletions)
    @test count(x -> x == 0, mask) == 2 * deletions

    @test_throws DomainError randommask(-1, 10)
    @test_throws DomainError randommask(-1, Bernoulli(0.5))
    @test_throws DomainError randommask(2, 0)

    n = 20
    fraction = 0.2
    mask = randommask(n, fraction)
    @test count(x -> x == 0, mask)/(prod(size(mask)) - n) == fraction

    # idxmask
    D = EDM(2, 10)
    @test_throws BoundsError idxmask(D, CartesianIndices((9:11, 9:11)))

    D = EDM(2, 10)
    mask = [1  1  1  1  1  1  1  1  1  1;
            1  1  0  0  1  1  1  1  1  1;
            1  0  1  0  1  1  1  1  1  1;
            1  0  0  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1]
    @test all(mask .== idxmask(D, CartesianIndices((2:4, 2:4))))

    mask = [1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  0  0  0  0;
            1  1  1  1  1  1  0  0  0  0;
            1  1  1  1  1  1  0  0  0  0;
            1  1  1  1  1  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1]
    @test all(mask .== idxmask(D, CartesianIndices((7:10, 3:5)), symmetric=true))
    @test all(mask .== idxmask(D, CartesianIndices((7:10, 3:5))))

    mask = [1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  1  1  1  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1;
            1  1  0  0  0  1  1  1  1  1]
    @test all(mask .== idxmask(D, CartesianIndices((7:10, 3:5)), symmetric=false))
end


text:
"""A full basin scenario example replicating the Lower Campaspe subcatchment."""

using DataFrames, CSV, FileIO
using Statistics, OrderedCollections
using Distributed
using Surrogates
using Agtor, Dates

using Interact: @manipulate

import Base.Filesystem: rm

"""Reallocate water on given date."""
function allocation_precall!(zone, dt_i)
    # Resetting allocations for each growing season
    if monthday(dt_i) == gs_start
        # Full annual water allocations each year
        tmp = LittleDict(ws.name=>ws.entitlement for ws in zone.water_sources)
        t_names = Tuple(Symbol.(collect(keys(tmp))))
        allocs = NamedTuple{t_names}(collect(values(tmp)))

        update_available_water!(zone, allocs)
    end
end

function scenario_run(samples, basin; pre=nothing, post=nothing)::Dict
    # Need to change this to be a Dict so we have scenario/
    scenario_results = Dict()
    sizehint!(scenario_results, nrow(samples))

    for (row_id, r) in enumerate(eachrow(samples))
        tmp_b = deepcopy(basin)

        update_model!(tmp_b, r)
        scenario_results[row_id] = run_model(tmp_b, run_timestep!; pre=pre, post=post)
    end

    return scenario_results
end

# Start month/day of growing season
# Represents when water allocations are announced.
const gs_start = (5, 1)

data_dir = "test/data/"

# Basin setup
basin_spec_dir = "examples/campaspe/basin/"
basin_spec = load_spec(basin_spec_dir)[:Campaspe]

basin_name = basin_spec[:name]
zone_specs = basin_spec[:zone_spec]

OptimizingManager = BaseManager("optimizing")

climate_data = "examples/campaspe/climate/basin_historic_climate_data.csv"

# Set the optimizing manager as the farmer for all zones
manage_zones = ((OptimizingManager, Tuple(collect(keys(zone_specs)))), )
campaspe_basin = Basin(name=basin_name, zone_spec=zone_specs,
                       climate_data=climate_data, managers=manage_zones)

agparams = collect_agparams!(campaspe_basin)

samples = sample(50, agparams[:, :min_val], agparams[:, :max_val], SobolSample())

# Match sampled values with parameter names
df = names!(DataFrame(samples), map(Symbol, agparams[:, :name]))

res = scenario_run(df, campaspe_basin; pre=allocation_precall!)

save_results!("campaspe_example.jld2", res)

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
module AdventOfCodeSolutions

struct Puzzle{year, day, part} end
Puzzle(day, part) = Puzzle{2020, day, part}()
Puzzle(year, day, part) = Puzzle{year, day, part}()

function openInput(::Puzzle{n, d, p}) where {n, d, p}
    return open(joinpath(@__DIR__, "../$n/data/day_$d.txt"), "r")
end

export Puzzle, openInput
end


text:
export image, preimage, has_preimage, isinjective,
issurjective, isbijective, automorphisms, morphisms, find_small_group,
inv, *, id_hom, domain, codomain, divisors, multiples

################################################################################
#
#  Morphism functionality
#
################################################################################

@doc Markdown.doc"""
    image(f::GrpGenToGrpGenMor, g::GrpGenElem) -> h::GrpGenElem

Returns the image of $g$ under $f$.
"""
image(f::GrpGenToGrpGenMor, g::GrpGenElem) = f.img[g.i]

@doc Markdown.doc"""
    preimage(f::GrpGenToGrpGenMor, g::GrpGenElem) -> h::GrpGenElem

Returns one element of the preimage of $g$ under $f$.
"""
function preimage(f::GrpGenToGrpGenMor, g::GrpGenElem)
  h = findfirst(x -> f(x) == g, collect(f.domain))
   if h == nothing
     error("$g has no preimage under $f")
   end
   return f.domain[h]
end

@doc Markdown.doc"""
    has_preimage(f::GrpGenToGrpGenMor, g::GrpGenElem) -> (b::Bool, h::GrpGenElem)

Returns whether $g$ has a preimage under $f$. If so, the second return value is an
element $h$ with $f(h) = g$.
"""
function has_preimage(f::GrpGenToGrpGenMor, g::GrpGenElem)
  h = findfirst(x -> f(x) == g, collect(f.domain))
   if h == nothing
     return false
   end
   return (true, f.domain[h])
end

@doc Markdown.doc"""
    *(f::GrpGenToGrpGenMor, g::GrpGenToGrpGenMor) -> h::GrpGenToGrpGenMor

Returns the composition $(f * g) = g(f)$.
"""
function *(f::GrpGenToGrpGenMor, g::GrpGenToGrpGenMor)
  return GrpGenToGrpGenMor(f.domain, g.codomain, [g(f(x)) for x in collect(f.domain)])
end

@doc Markdown.doc"""
    inv(f::GrpGenToGrpGenMor) -> h::GrpGenToGrpGenMor

Assumes that $f$ is an isomorphism. Returns the inverse of $f$.
"""
function inv(f::GrpGenToGrpGenMor)
  return GrpGenToGrpGenMor(f.codomain, f.domain, collect(f.domain)[sortperm(getindex.(f.img))])
end

function Base.show(io::IO, f::GrpGenToGrpGenMor)
  println(io, "Morphism from group\n", f.domain, " to\n", f.codomain)
end

domain(f::GrpGenToGrpGenMor) = f.domain

codomain(f::GrpGenToGrpGenMor) = f.codomain

id_hom(G::GrpGen) = GrpGenToGrpGenMor(G, G, collect(G))

image(GtoH::GrpGenToGrpGenMor) = subgroup(GtoH.codomain, unique(GtoH.img))

function kernel(GtoH::GrpGenToGrpGenMor)
  G = GtoH.domain
  H = GtoH.codomain
  return subgroup(G, getindex.(Ref(G), findall(x-> GtoH(x) == id(H), collect(G))))
end

function issurjective(GtoH::GrpGenToGrpGenMor)
  return order(GtoH.codomain) == length(unique(GtoH.img))
end
#finite groups
isinjective(GtoH::GrpGenToGrpGenMor) = issurjective(GtoH)

isbijective(GtoH::GrpGenToGrpGenMor) = issurjective(GtoH)

(M::GrpGenToGrpAbMor)(g::GrpGenElem) = M.dict[g]

function (M::GrpAbToGrpGenMor)(g::GrpAbElem)
  elem = id(M.codomain)
  for i in 1:length(g.coeff)
    elem = elem * M.I[i]^g.coeff[i]
  end
  return elem
end

################################################################################
#
#  Find Morphisms
#
################################################################################

# TODO: Cache the orders of the generators of the small_groups.
#       Do not recompute it here.
function find_small_group(G::GrpGen; DB = DefaultSmallGroupDB())
  l = order(G)

  D = DB.db

  elements_by_orders = Dict{Int, Array{GrpGenElem, 1}}()

  for i in 1:l
    g = G[i]
    o = order(g)
    if haskey(elements_by_orders, o)
      push!(elements_by_orders[o], g)
    else
      elements_by_orders[o] = [g]
    end
  end

  candidates = Int[]

  local ordershape

  for j in 1:length(D[order(G)])
    ordershape = D[order(G)][j].orderdis

    candidate = true
    for (o, no) in ordershape
      if !haskey(elements_by_orders, o)
        candidate = false
        break
      else
        elts = elements_by_orders[o]
        if length(elts) != no
          candidate = false
          break
        end
      end
     end

     if candidate
        push!(candidates, j)
     end
  end

  @assert length(candidates) > 0

  sort!(candidates, rev = true)

  idG = id(G)

  for j in candidates
    H = D[order(G)][j]

    elbyord = [elements_by_orders[order(o)] for o in H.gens]

    it = Iterators.product(elbyord...)

    words = H.rels

    for poss in it
      is_hom = true
      for w in words
        if eval_word(collect(poss), w) != idG
          is_hom = false
          break
        end
      end

      if is_hom
        if length(closure(collect(poss), *, idG)) == order(G)
          # Found it!
          H = small_group(order(G), j, DB = DB)
          return (order(G), j), H, _spin_up_morphism(gens(H), collect(poss))
        end
      end
    end
  end
  throw(error("Could not identify group"))
end

function eval_word(S, w::Vector{Int})
  g = id(parent(S[1]))
  for i in 1:length(w)
    if w[i] > 0
      g = g * S[w[i]]
    else
      g = g * inv(S[-w[i]])
    end
  end
  return g
end

@doc Markdown.doc"""
    automorphisms(G::GrpGen) -> A::Array{GrpGenToGrpGenMor,1}

Returns all group isomorphisms from $G$ to $G$.
"""
function automorphisms(G::GrpGen)
  Gn, GntoG = find_small_group(G)[2:3]
  auts = _automorphisms(Gn)
  return [inv(GntoG) * aut * GntoG for aut in auts]
end

function _automorphisms(G::GrpGen)
  @assert isfrom_db(G)
  i, j = G.small_group_id
  Gdata = DefaultSmallGroupDB().db[i][j]

  l = order(G)

  elements_by_orders = Dict{Int, Array{GrpGenElem, 1}}()

  # TODO: I think the following is cached somewhere (in the database)
  for i in 1:l
    g = G[i]
    o = order(g)
    if haskey(elements_by_orders, o)
      push!(elements_by_orders[o], g)
    else
      elements_by_orders[o] = [g]
    end
  end

  elbyord = [elements_by_orders[order(G[o])] for o in G.gens]

  it = Iterators.product(elbyord...)

  words::Vector{Vector{Int}} = Gdata.rels

  idG = id(G)

  auts = _aut_group(it, words, idG, order(G))::Vector{Vector{GrpGenElem}}

  # Any element A of auts determines an isomorphism by mapping gens(G)[i] to A[i]

  Ggens = gens(G)

  # TODO: preallocate
  return [_spin_up_morphism(Ggens, a) for a in auts]
end

function _spin_up_morphism(domain::Vector{GrpGenElem}, codomain::Vector{GrpGenElem})
  @assert length(domain) > 0
  @assert length(domain) == length(codomain)
  G = parent(domain[1])
  H = parent(codomain[1])
  pairs = [(domain[i], codomain[i]) for i in 1:length(domain)]
  cl = closure(pairs, (x, y) -> (x[1]*y[1], x[2]*y[2]), (id(G), id(H)))
  img = Vector{GrpGenElem}(undef, length(G))
  for i in 1:length(G)
    img[cl[i][1][]] = cl[i][2]
  end
  phi = GrpGenToGrpGenMor(G, H, img)

  # TODO: Remove this assertion once this is battle tested
  for g in G
    for h in G
      @assert phi(g * h) == phi(g) * phi(h)
    end
  end
  return phi
end

@noinline function _aut_group(it, words, idG, n)
  auts = Vector{GrpGenElem}[]
  for poss in it
    is_hom = true
    for w in words
      if eval_word(poss, w) != idG
        is_hom = false
        break
      end
    end

    if is_hom
      cposs = collect(poss)
      if length(closure(cposs, *, idG)) == n
        push!(auts, cposs)
      end
    end
  end

  return auts
end

function _morphisms_with_gens(G::GrpGen, H::GrpGen, Gens::Array{GrpGenElem,1}, Rels::Array{Array{Int64,1},1})

  l = order(H)

  elements_by_orders = Dict{Int, Array{GrpGenElem, 1}}()

  for i in 1:l
    h = H[i]
    o = order(h)
    for k in multiples(o, order(G))
      if haskey(elements_by_orders, k)
        push!(elements_by_orders[k], h)
      else
        elements_by_orders[k] = [h]
      end
    end
  end

  elbyord = [elements_by_orders[order(o)] for o in Gens]

  it = Iterators.product(elbyord...)

  words::Vector{Vector{Int}} = Rels

  idH = id(H)

  homs = _hom_group(it, words, idH)::Vector{Vector{GrpGenElem}}

  Ggens = Gens

  return [_spin_up_morphism(Ggens, a) for a in homs]
end

@doc Markdown.doc"""
    morphisms(G::GrpGen, H::GrpGen) -> A::Array{GrpGenToGrpGenMor,1}

Returns all group homomorphisms from $G$ to $H$.
"""
function morphisms(G::GrpGen, H::GrpGen)
  Gn, isom = find_small_group(G)[2:3]
  return [inv(isom) * mor for mor in _morphisms(Gn,H)]
end

function _morphisms(G::GrpGen, H::GrpGen)
  @assert isfrom_db(G)
  i, j = G.small_group_id
  Gdata = DefaultSmallGroupDB().db[i][j]

  l = order(H)

  elements_by_orders = Dict{Int, Array{GrpGenElem, 1}}()

  for i in 1:l
    h = H[i]
    o = order(h)
    for k in multiples(o, order(G))
      if haskey(elements_by_orders, k)
        push!(elements_by_orders[k], h)
      else
        elements_by_orders[k] = [h]
      end
    end
  end

  elbyord = [elements_by_orders[order(G[o])] for o in G.gens]

  it = Iterators.product(elbyord...)

  words::Vector{Vector{Int}} = Gdata.rels

  idH = id(H)

  homs = _hom_group(it, words, idH)::Vector{Vector{GrpGenElem}}

  # Any element a of homs determines an homomorphism by mapping gens(G)[i] to A[i]

  Ggens = gens(G)

  # TODO: preallocate
  return [_spin_up_morphism(Ggens, a) for a in homs]
end

@noinline function _hom_group(it, words, idH)
  homs = Vector{GrpGenElem}[]
  for poss in it
    is_hom = true
    for w in words
      if eval_word(poss, w) != idH
        is_hom = false
        break
      end
    end

    if is_hom
      cposs = collect(poss)
      push!(homs, cposs)
    end
  end
  return homs
end

divisors(n::Int64) =  findall(x -> mod(n,x) == 0, 1:n)

multiples(n::Int64, b::Int64) =  [i * n for i in 1:Int64(floor(b/n))]

function isisomorphic(G::GrpGen, H::GrpGen)
  idG, A, AtoG = find_small_group(G)
  idH, B, BtoH = find_small_group(H)
  if idG != idH
    return false, id_hom(G) # I am too lazy
  else
    h = _spin_up_morphism(gens(A), gens(B)) # they must be equal
    return true, _spin_up_morphism(gens(G), [ BtoH(h(preimage(AtoG, g))) for g in gens(G)])
  end
end


text:
function wplot(W, nzs, padding=10)
    h, w = IMSIZE
    n = size(W, 1)
    rows = int(floor(sqrt(n)))
    cols = int(ceil(n / rows))
    halfpad = div(padding, 2)
    dat = zeros(rows * (h + padding), cols * (w + padding))
    for i=1:n
        wt = W[i, :]
        wt = reshape(wt, length(wt))
        wim = map_nonzeros(IMSIZE, wt, nzs)
        wim = wim ./ (maximum(wim) - minimum(wim))
        r = div(i - 1, cols) + 1
        c = rem(i - 1, cols) + 1
        dat[(r-1)*(h+padding)+halfpad+1 : r*(h+padding)-halfpad,
            (c-1)*(w+padding)+halfpad+1 : c*(w+padding)-halfpad] = wim
    end
    view(dat)
    return dat
end

function svm_accuracy(X, y, k=10)
    svmest(idxs) = svmtrain(y[[idxs]], X[:, idxs])
    svmeval(model, idxs) = svmpredict(model, X[:, idxs])
    return cross_validate(svmest, svmeval, 10, Kfold(length(y), 10))

    ## fold_list = collect(Kfold(length(y), k))
    ## for fold_idx in fold_list

    ## end

end

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
using DifferentialEquations
using SimpleDiffEq
using Polynomials
using ModelingToolkit
using Symbolics
using Plots

M = 2
N = 10
f(s,i) = s*i;

@variables a[1:M,1:N];

An = [] # Empty array of Adomian
A₀ = f(a[1,1],a[2,1]) # s₀*i₀
push!(An,A₀)
for n = 1:(N-1)
  A = 0
  for i = 1:M
    for k = 1:n
      A += k*a[i,k+1]*Symbolics.derivative(An[end],a[i,k])
    end
  end
  A /= n
  push!(An,A)
end;

tspan = (0.0,8.0)
trange = 0:0.1:8;

u0 = [20.0 15.0]'; # S, I

p = [0.01, 0.02]; # β, γ

(S,I) = u0'
v = u0
β, γ = p'
Ii = integrate(Polynomial([I],:t))
F = eval(build_function(An[1],a))
SIi = integrate(Polynomial([F(v)],:t))
for i in 1:(N-1)
  vv = [Polynomial(ones(N),:t) Polynomial(ones(N),:t)]'
  vv[1,1] = -β*SIi
  vv[2,1] =  β*SIi - γ*Ii
  global v = hcat(v,vv)
  # Now update
  global I = v[2,i+1]
  global Ii = integrate(I)
  global F = eval(build_function(An[i+1],a))
  global SIi = integrate(F(v))
end

S = sum(v[1,1:end])

I = sum(v[2,1:end])

Sa = [S(t) for t in trange]
Ia = [I(t) for t in trange];

plot(trange,Sa,label="S")
plot!(trange,Ia,label="I")

function sir_ode!(du,u,p,t)
    (S,I) = u
    (β,γ) = p
    @inbounds begin
        du[1] = -β*S*I
        du[2] = β*S*I - γ*I
    end
    nothing
end
prob_ode = ODEProblem(sir_ode!,u0,tspan,p)
sol_ode = solve(prob_ode)
plot(sol_ode,
     label=["S" "I"],
     xlabel="Time",
     ylabel="Number")

Fv = [eval(build_function(An[i],a)) for i in 1:N];

function SIR(u0,p,trange,Fv)
  (S,I) = u0'
  v = u0
  β, γ = p'
  Ii = integrate(Polynomial([I],:t))
  SIi = integrate(Polynomial([Fv[1](v)],:t))
  for i in 1:(N-1)
    vv = [Polynomial(ones(N),:t) Polynomial(ones(N),:t)]'
    vv[1,1] = -β*SIi
    vv[2,1] =  β*SIi - γ*Ii
    v = hcat(v,vv)
    # Now update
    I = v[2,i+1]
    Ii = integrate(I)
    SIi = integrate(Fv[i+1](v))
  end
  S = sum(v[1,1:end])
  I = sum(v[2,1:end])
  Sa = [S(t) for t in trange]
  Ia = [I(t) for t in trange]
  return [trange Sa Ia]
end;

u0 = [990.0 10.0]'
p = [0.0005 0.25]
trange = 0:0.1:40
sol1 = SIR(u0,p,trange,Fv)
plot(sol1[1:end,1],sol1[1:end,2:3],
     label=["S" "I"],
     xlabel="Time",
     ylabel="Number")

u = u0
t = 0.0
sol = [[0.0 u']]
for i in 1:40
  s = SIR(u,p,0.0:0.1:1.0,Fv)
  s[1:end,1] .+= t
  global t += 1.0
  push!(sol,s[2:end,:])
  global u = s[end,2:3]
end
sol = vcat(sol...);

plot(sol[1:end,1],sol[1:end,2:3],
     label=["S" "I"],
     xlabel="Time",
     ylabel="Number")

prob_ode2 = ODEProblem(sir_ode!,u0,(0,40.0),p)
sol_ode2 = solve(prob_ode2)
plot(sol_ode2,
     label=["S" "I"],
     xlabel="Time",
     ylabel="Number")


text:
using RecipesBase
using Setfield
getLensSymbol(lens) = :p
getLensSymbol(lens::Setfield.PropertyLens{F}) where F = F
getLensSymbol(lens::Setfield.ComposedLens) = getLensSymbol(lens.inner)
getLensSymbol(::Setfield.IdentityLens) = :p
getLensSymbol(::Setfield.IndexLens{Tuple{Int64}}) = :p
getLensSymbol(br::AbstractBranchResult) = getLensSymbol(br.lens)

function getLensSymbol(lens1::Lens, lens2::Lens)
	p1 = getLensSymbol(lens1)
	p2 = getLensSymbol(lens2)
	out = p1 == p2 ? (Symbol(String(p1)*"1"), Symbol(String(p2)*"2")) : (p1, p2)
end

function getPlotVars(contres, vars)
	if vars isa Tuple{Symbol, Symbol} || typeof(vars) <: Tuple{Int64, Int64}
		return vars
	else
		return :param, getfirstusertype(contres)
	end
end

# https://github.com/JuliaGraphics/Colors.jl/blob/master/src/names_data.jl
# we don't need to have different colors. Indeed, some bifurcations never occurs together codim1 ∩ codim2 = ∅
const colorbif = Dict(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta, :none => :yellow, :ns => :orange, :pd => :green, :bt => :red, :cusp => :sienna1, :gh => :brown, :zh => :burlywood2, :hh => :green, :user => :darkgoldenrod)

function getColor(sp)
	if sp in keys(colorbif)
		return colorbif[sp]
	else
		return :darkgoldenrod
	end
end

function getAxisLabels(ind1, ind2, br)
	xguide = ""
	yguide = ""
	if ind1 == 1 || ind1 == :param
		xguide = String(getLensSymbol(br))
	elseif ind1 isa Symbol
		xguide = String(ind1)
	end
	if ind2 isa Symbol
		yguide = String(ind2)
	end
	return xguide, yguide
end

# allow to plot a single branch
RecipesBase.@recipe function Plots(contres::AbstractBranchResult;
	plotfold = false,
	putspecialptlegend = true,
	filterspecialpoints = false,
	vars = nothing,
	plotstability = true,
	plotspecialpoints = true,
	branchlabel = "",
	linewidthunstable = 1.0,
	linewidthstable = 2linewidthunstable,
	plotcirclesbif = false,
	applytoY = identity,
	applytoX = identity)
	# Special case labels when vars = (:p,:y,:z) or (:x) or [:x,:y] ...
	ind1, ind2 = getPlotVars(contres, vars)
	xlab, ylab = getAxisLabels(ind1, ind2, contres)
	@series begin
		if hasstability(contres) && plotstability
			linewidth --> map(x -> isodd(x) ? linewidthstable : linewidthunstable, contres.stable)
		end
		xguide --> xlab
		yguide --> ylab
		label --> branchlabel
		map(applytoX, getproperty(contres.branch, ind1)), map(applytoY, getproperty(contres.branch, ind2))
	end

	# display bifurcation points
	bifpt = filter(x -> (x.type != :none) && (plotfold || x.type != :fold) && (x.idx <= length(contres)-1), contres.specialpoint)

	if length(bifpt) >= 1 && plotspecialpoints #&& (ind1 == :param)
		if filterspecialpoints == true
			bifpt = filterBifurcations(bifpt)
		end
		@series begin
			seriestype := :scatter
			seriescolor --> map(x -> getColor(x.type), bifpt)
			markershape --> map(x -> (x.status != :converged) && (plotcirclesbif==false) ? :square : :circle, bifpt)
			markersize --> 3
			markerstrokewidth --> 0
			label --> ""
			xguide --> xlab
			yguide --> ylab
			[applytoX(getproperty(contres[pt.idx], ind1)) for pt in bifpt], [applytoY(getproperty(contres[pt.idx], ind2)) for pt in bifpt]
		end
		# add legend for bifurcation points
		if putspecialptlegend && length(bifpt) >= 1
			bps = unique(x -> x.type, [pt for pt in bifpt if (pt.type != :none && (plotfold || pt.type != :fold))])
			(length(bps) == 0) && return
			for pt in bps
				@series begin
					seriestype := :scatter
					seriescolor --> getColor(pt.type)
					label --> "$(pt.type)"
					markersize --> 3
					markerstrokewidth --> 0
					xguide --> xlab
					yguide --> ylab
					[applytoX(getproperty(contres[pt.idx], ind1))], [applytoY(getproperty(contres[pt.idx], ind2))]
				end
			end
		end
	end
end

# allow to plot branches specified by splatting
RecipesBase.@recipe function Plots(brs::AbstractBranchResult...;
		plotfold = false,
		putspecialptlegend = true,
		filterspecialpoints = false,
		vars = nothing,
		plotstability = true,
		plotspecialpoints = true,
		branchlabel = fill("",length(brs)),
		linewidthunstable = 1.0,
		linewidthstable = 2linewidthunstable,
		applytoY = identity,
		applytoX = identity)
	ind1, ind2 = getPlotVars(brs[1], vars)
	if length(brs) == 0; return; end
	# handle bifurcation points, the issue is to simplify the legend. So we collect all bifurcation points
	bptps = Any[(type = pt.type, param = getproperty(pt, ind1), printsol = getproperty(pt.printsol, ind2)) for pt in brs[1].specialpoint if pt.type != :none]
	for ii=2:length(brs)
		_ind1, _ind2 = getPlotVars(brs[ii], vars)
		for pt in brs[ii].specialpoint
			if pt.type != :none
				push!(bptps, (type = pt.type, param = getproperty(pt, _ind1), printsol = getproperty(pt.printsol, _ind2)))
			end
		end
	end
	bp = unique(x -> x.type, bptps)

	# add legend for bifurcation points
	if putspecialptlegend && length(bp) > 0
		for pt in unique(x -> x.type, bp)
			@series begin
				seriestype := :scatter
				seriescolor --> getColor(pt.type)
				label --> "$(pt.type)"
				markersize --> 3
				markerstrokewidth --> 0
				[applytoX(pt.param)], [applytoY(pt.printsol)]
			end
		end
	end

	for (id, res) in pairs(brs)
		@series begin
			putspecialptlegend --> false
			plotfold --> plotfold
			plotspecialpoints --> plotspecialpoints
			plotstability --> plotstability
			branchlabel --> branchlabel[id]
			linewidthunstable --> linewidthunstable
			linewidthstable --> linewidthstable
			applytoX --> applytoX
			applytoY --> applytoY
			vars --> vars
			if ind1 == 1 || ind1 == :param
				xguide --> String(getLensSymbol(brs[id].lens))
			elseif ind1 isa Symbol
				xguide --> String(ind1)
			end
			if ind2 isa Symbol
				yguide --> String(ind2)
			end
			# collect the values of the bifurcation points to be added in the legend
			ind1, ind2 = getPlotVars(brs[id], vars)
			for pt in res.specialpoint
				# this does not work very well when changing optional argument vars
				pt.type != :none && push!(bp, (type = pt.type, param = getproperty(pt, :param), printsol = getproperty(pt.printsol, ind2)))
			end
			res
		end
	end
end
####################################################################################################
function filterBifurcations(bifpt)
	# this function filters Fold points and Branch points which are located at the same/previous/next point
	length(bifpt) == 0 && return bifpt
	res = [(type = :none, idx = 1, param = 1., printsol = bifpt[1].printsol, status = :guess)]
	ii = 1
	while ii <= length(bifpt) - 1
		if (abs(bifpt[ii].idx - bifpt[ii+1].idx) <= 1) && bifpt[ii].type ∈ [:fold, :bp]
			if (bifpt[ii].type == :fold && bifpt[ii].type == :bp) ||
				(bifpt[ii].type == :bp && bifpt[ii].type == :fold)
				push!(res, (type = :fold, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol, status = bifpt[ii].status) )
			else
				push!(res, (type = bifpt[ii].type, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol, status = bifpt[ii].status) )
				push!(res, (type = bifpt[ii+1].type, idx = bifpt[ii+1].idx, param = bifpt[ii+1].param, printsol = bifpt[ii+1].printsol,status = bifpt[ii].status) )
			end
			ii += 2
		else
			push!(res, (type = bifpt[ii].type, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol, status = bifpt[ii].status) )
			ii += 1
		end
	end
	0 < ii <= length(bifpt) &&	push!(res, (type = bifpt[ii].type, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol, status = bifpt[ii].status) )

	return res[2:end]
end
####################################################################################################
# plot recipes for the bifurcation diagram
RecipesBase.@recipe function f(bd::Vector{BifDiagNode}; code = (), level = (-Inf, Inf))
	for b in bd
		@series begin
			level --> level
			code --> code
			b
		end
	end
end

RecipesBase.@recipe function f(bd::BifDiagNode; code = (), level = (-Inf, Inf))
	if ~hasbranch(bd); return; end
	_bd = getBranch(bd, code)
	@series begin
		level --> level
		code --> ()
		_bd.child
	end
	# !! plot root branch in last so the bifurcation points do not alias, for example a 2d BP would be plot as a 1d BP if the order were reversed
	if level[1] <= _bd.level <= level[2]
		@series begin
			_bd.γ
		end
	end
end

# this might well be type piracy
# TODO try to remove it
RecipesBase.@recipe function f(bd::Nothing)
	nothing
end

####################################################################################################
# plotting function of the periodic orbits
function plotPeriodicPOTrap(x, M, Nx, Ny; ratio = 2, kwargs...)
	@assert ratio > 0 "You need at least one component"
	n = Nx*Ny
	outpo = reshape(x[1:end-1], ratio * n, M)
	po = reshape(x[1:n,1], Nx, Ny)
	rg = 2:6:M
	for ii in rg
		po = hcat(po, reshape(outpo[1:n,ii], Nx, Ny))
	end
	heatmap!(po; color = :viridis, fill=true, xlabel = "space", ylabel = "space", kwargs...)
	for ii in 1:length(rg)
		plot!([ii*Ny, ii*Ny], [1, Nx]; color = :red, width = 3, label = "", kwargs...)
	end
end

function plotPeriodicPOTrap(outpof, n, M; ratio = 2)
	@assert ratio > 0 "You need at least one component"
	outpo = reshape(outpof[1:end-1], ratio * n, M)
	if ratio == 1
		heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis)
	else
		plot(heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis),
			heatmap(outpo[n+2:end,:]', color=:viridis))
	end
end
####################################################################################################
function plotPeriodicShooting!(x, M; kwargs...)
	N = div(length(x), M);	plot!(x; label = "", kwargs...)
	for ii in 1:M
		plot!([ii*N, ii*N], [minimum(x), maximum(x)] ;color = :red, label = "", kwargs...)
	end
end

function plotPeriodicShooting(x, M; kwargs...)
	plot();plotPeriodicShooting!(x, M; kwargs...)
end
####################################################################################################
# plot recipe for codim 2 plot
# TODO Use dispatch for this
RecipesBase.@recipe function Plots(contres::ContResult{Ta, Teigvals, Teigvec, Biftype, Ts, Tparc, Tfunc, Tpar, Tl};
	plotfold = false,
	putspecialptlegend = true,
	filterspecialpoints = false,
	vars = nothing,
	plotstability = true,
	plotspecialpoints = true,
	branchlabel = "",
	linewidthunstable = 1.0,
	linewidthstable = 2linewidthunstable,
	plotcirclesbif = false,
	_basicplot = true,
	applytoY = identity,
	applytoX = identity) where {Ta, Teigvals, Teigvec, Biftype, Ts, Tparc, Tfunc <: Union{FoldProblemMinimallyAugmented, HopfProblemMinimallyAugmented}, Tpar, Tl}
	# Special case labels when vars = (:p,:y,:z) or (:x) or [:x,:y] ...
	ind1, ind2 = getPlotVars(contres, vars)
	xlab, ylab = getAxisLabels(ind1, ind2, contres)
	@series begin
		if hasstability(contres) && false
			linewidth --> map(x -> isodd(x) ? linewidthstable : linewidthunstable, contres.stable)
		end
		xguide --> xlab
		yguide --> ylab
		label --> branchlabel
		map(applytoX, getproperty(contres.branch, ind1)), map(applytoY, getproperty(contres.branch, ind2))
	end

	# display bifurcation points
	bifpt = filter(x -> (x.type != :none) && (plotfold || x.type != :fold), contres.specialpoint)

	if length(bifpt) >= 1 && plotspecialpoints #&& (ind1 == :param)
		if filterspecialpoints == true
			bifpt = filterBifurcations(bifpt)
		end
		@series begin
			seriestype := :scatter
			seriescolor --> map(x -> getColor(x.type), bifpt)
			markershape --> map(x -> (x.status != :converged) && (plotcirclesbif==false) ? :square : :circle, bifpt)
			markersize --> 3
			markerstrokewidth --> 0
			label --> ""
			xguide --> xlab
			yguide --> ylab
			[applytoX(getproperty(contres[pt.idx], ind1)) for pt in bifpt], [applytoY(getproperty(contres[pt.idx], ind2)) for pt in bifpt]
		end
		# add legend for bifurcation points
		if putspecialptlegend && length(bifpt) >= 1
			bps = unique(x -> x.type, [pt for pt in bifpt if (pt.type != :none && (plotfold || pt.type != :fold))])
			(length(bps) == 0) && return
			for pt in bps
				@series begin
					seriestype := :scatter
					seriescolor --> getColor(pt.type)
					label --> "$(pt.type)"
					markersize --> 3
					markerstrokewidth --> 0
					xguide --> xlab
					yguide --> ylab
					[applytoX(getproperty(contres[pt.idx], ind1))], [applytoY(getproperty(contres[pt.idx], ind2))]
				end
			end
		end
	end
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

suitable univariate convex functions defined on ℝ₊₊ and associated oracles
- h_val evaluates (sum of) h
- h_conj_dom_pos is true if domain of conjugate of h is ℝ₊₊, else domain is ℝ
- h_conj evaluates the convex conjugate of (sum of) h
- h_der1/h_der2/h_der3 evaluate 1st/2nd/3rd derivatives of h pointwise
- get_initial_point gives initial point for u,v,λ_i

TODO derive central initial points for get_initial_point and use real type
=#

"""
$(TYPEDEF)

The inverse function ``x \\to x^{-1}``.
"""
struct InvSSF <: SepSpectralFun end

h_val(xs::Vector, ::InvSSF) = sum(inv, xs)
h_conj_dom_pos(::InvSSF) = true
h_conj(xs::Vector, ::InvSSF) = -2 * sum(sqrt, xs)

h_der1(ds::Vector, xs::Vector, ::InvSSF) = (@. ds = -xs^-2)
h_der2(ds::Vector, xs::Vector, ::InvSSF) = (@. ds = 2 * xs^-3)
h_der3(ds::Vector, xs::Vector, ::InvSSF) = (@. ds = -6 * xs^-4)

function get_initial_point(d::Int, ::InvSSF)
    # TODO initial central point
    return (2d, 1, 1)
end

"""
$(TYPEDEF)

The negative logarithm function ``x \\to - \\log(x)``.
"""
struct NegLogSSF <: SepSpectralFun end

h_val(xs::Vector, ::NegLogSSF) = -sum(log, xs)
h_conj_dom_pos(::NegLogSSF) = true
h_conj(xs::Vector, ::NegLogSSF) = -length(xs) - sum(log, xs)

h_der1(ds::Vector, xs::Vector, ::NegLogSSF) = (@. ds = -inv(xs))
h_der2(ds::Vector, xs::Vector, ::NegLogSSF) = (@. ds = xs^-2)
h_der3(ds::Vector, xs::Vector, ::NegLogSSF) = (@. ds = -2 * xs^-3)

function get_initial_point(d::Int, ::NegLogSSF)
    # TODO initial central point
    return (1, 1, 1)
end

"""
$(TYPEDEF)

The negative entropy function ``x \\to x \\log(x)``.
"""
struct NegEntropySSF <: SepSpectralFun end

h_val(xs::Vector, ::NegEntropySSF) = sum(x * log(x) for x in xs)
h_conj_dom_pos(::NegEntropySSF) = false
h_conj(xs::Vector, ::NegEntropySSF) = sum(exp(-x - 1) for x in xs)

h_der1(ds::Vector, xs::Vector, ::NegEntropySSF) = (@. ds = 1 + log(xs))
h_der2(ds::Vector, xs::Vector, ::NegEntropySSF) = (@. ds = inv(xs))
h_der3(ds::Vector, xs::Vector, ::NegEntropySSF) = (@. ds = -xs^-2)

function get_initial_point(d::Int, ::NegEntropySSF)
    # TODO initial central point
    return (1, 1, 1)
end

"""
$(TYPEDEF)

The power function ``x \\to x^p`` parametrized by ``p \\in (1, 2]``. Note for
``p = 2``, it is more efficient to use [`EpiPerSquare`](@ref) on the scaled
triangle.
"""
struct Power12SSF <: SepSpectralFun
    p::Real
    Power12SSF(p::Real) = (@assert 1 < p <= 2; new(p))
end

h_val(xs::Vector, h::Power12SSF) = sum(x^h.p for x in xs)
h_conj_dom_pos(::Power12SSF) = false
function h_conj(xs::Vector, h::Power12SSF)
    p = h.p
    q = p / (p - 1)
    return (p - 1) * sum((x >= 0 ? zero(x) : (abs(x) / p)^q) for x in xs)
end

function h_der1(ds::Vector, xs::Vector, h::Power12SSF)
    p = h.p
    pm1 = p - 1
    @. ds = p * xs^pm1
    return ds
end
function h_der2(ds::Vector, xs::Vector, h::Power12SSF)
    p = h.p
    pm2 = p - 2
    coef = p * (p - 1)
    @. ds = coef * xs^pm2
    return ds
end
function h_der3(ds::Vector, xs::Vector, h::Power12SSF)
    p = h.p
    pm3 = p - 3
    coef = p * (p - 1) * (p - 2)
    @. ds = coef * xs^pm3
    return ds
end

function get_initial_point(d::Int, h::Power12SSF)
    # TODO initial central point
    return (2d, 1, 1)
end


text:
export SomaOrDendrite, DendriteSegment, Neuron, Synapse, NeuralNetwork, SynapseType, excitatory, inhibitory
using DataStructures
abstract type SomaOrDendrite{T} end

mutable struct DefaultParams{T}
    children::Vector{Ref{DefaultParams{T}}}
    epsp_duration::T
    ipsp_duration::T
    plateau_duration::T
    spike_duration::T
    dendritic_threshold::UInt
    synaptic_threshold::UInt
    axon_delay::T
    refractory_period::T
    p_trans::Union{Float64,Vector{Float64}}
    _manually_set::DefaultDict{Symbol,Bool}
    function DefaultParams{T}(; kwargs...) where T
        obj = new{T}(Ref{DefaultParams{T}}[], T(1),T(1),T(100),T(1), 1, 1, T(0), T(1), 1.0, DefaultDict{Symbol,Bool}(false))

        for (key,value) in kwargs
            setproperty!(obj, key, value)
            obj._manually_set[key] = true
        end
        return obj
    end
end

function inherit!(p::DefaultParams, other::DefaultParams; connect=true)
    for (key,value) in other._manually_set
        if !p._manually_set[key]
            setproperty!(p, key, getproperty(other, key))
        end
    end
    for child in p.children
        inherit!(child[], p; connect=false)
    end

    if connect
        push!(other.children, Ref(p))
    end
    nothing
end

##
struct DendriteSegment{T} <: SomaOrDendrite{T}
    children::Vector{DendriteSegment{T}}
    parent::Ref{SomaOrDendrite{T}}
    neuron::Ref{SomaOrDendrite{T}}
    exc_inputs::Vector{Symbol}
    inh_inputs::Vector{Symbol}
    name::Symbol
    parameters::DefaultParams{T}
    function DendriteSegment{T}(children::Vector{DendriteSegment{T}}, exc_inputs, inh_inputs, name; kwargs...) where T

        this = new{T}(children, Ref{SomaOrDendrite{T}}(), Ref{SomaOrDendrite{T}}(), exc_inputs, inh_inputs, name, DefaultParams{T}(;kwargs...))
        for child in children
            child.parent[]=this
            inherit!(child.parameters, this.parameters)
        end
        this
    end
end

function DendriteSegment{T}(children, exc_inputs, inh_inputs, name, kwargs=Dict{Symbol,Any}()) where T
    new_children = DendriteSegment{T}[DendriteSegment{T}(child...) for child in children]
    invoke(DendriteSegment{T}, Tuple{Vector{DendriteSegment{T}},Vector{Symbol},Vector{Symbol},Symbol}, new_children, Vector{Symbol}(exc_inputs), Vector{Symbol}(inh_inputs), Symbol(name); kwargs...)
end
##
struct Neuron{T} <: SomaOrDendrite{T}
    all_segments::Vector{DendriteSegment{T}}
    dendrite::DendriteSegment{T}
    ports::Dict{Symbol,Vector{DendriteSegment{T}}}
    name::Ref{Symbol}
    parameters::DefaultParams{T}

    function Neuron{T}(args...; kwargs...) where T
        dendrite = DendriteSegment{T}(args...)
        params = DefaultParams{T}(;kwargs...)

        # store all segments in BFS-order
        head=0
        all_segments = [dendrite]
        while head < length(all_segments)
            head+=1
            d=all_segments[head]
            append!(all_segments, d.children)
        end

        ports = Dict{Symbol,Vector{DendriteSegment{T}}}()
        for segment in all_segments
            for input in segment.exc_inputs
                pp=get(ports, input, DendriteSegment{T}[])
                push!(pp, segment)
                ports[input]=pp
            end
            for input in segment.inh_inputs
                pp=get(ports, input, DendriteSegment{T}[])
                push!(pp, segment)
                ports[input]=pp
            end
            ports[Symbol("inh_$(segment.name)")] = [segment]
        end

        this = new{T}(all_segments,dendrite,ports, Ref(:neuron), params)

        function recursively_set_neuron!(d::DendriteSegment{T}, n::Neuron{T})
            d.neuron[] = n
            for child in d.children
                recursively_set_neuron!(child, n)
            end
        end

        dendrite.parent[]=this
        recursively_set_neuron!(dendrite, this)
        inherit!(dendrite.parameters, this.parameters)

        this
    end
end

##
struct Axon{T}
    source::Union{Ref{Neuron{T}},Symbol}
    target::Union{Tuple{Ref{Neuron{T}},Symbol},Symbol}
    parameters::DefaultParams{T}
    function Axon{T}(source, target; kwargs...) where T
        params = DefaultParams{T}(;kwargs...)
        if isa(source, Neuron)
            inherit!(params, source.parameters)
        end
        new{T}(source, target, params)
    end
end

struct NeuralNetwork{T}
    neurons::OrderedDict{Symbol,Neuron{T}}
    inputs::Vector{Symbol}
    outputs::Vector{Symbol}
    axons::OrderedDict{Symbol,Vector{Axon{T}}}
    parameters::DefaultParams{T}
    function NeuralNetwork{T}(neurons, inputs, outputs, axons; kwargs...) where T
        params = DefaultParams{T}(;kwargs...)

        for (name,neuron) in pairs(neurons)
            inherit!(neuron.parameters,params)
            for axon in get(axons, name, [])
                inherit!(axon.parameters,neuron.parameters)
            end
        end
        for axon_bundle in values(axons)
            for axon in axon_bundle
                inherit!(axon.parameters,params)
            end
        end

        new{T}(neurons, inputs, outputs, axons, params)
    end
end

function NeuralNetwork(neurons::OrderedDict{Symbol,Neuron{T}}, inputs::Vector{Symbol}, outputs::Vector{Symbol}, axons::OrderedDict{Symbol,Vector{X}}; kwargs...) where {T,X <:Union{Tuple,NTuple,NamedTuple}}
    for (name,neuron) in pairs(neurons)
        neuron.name[] = name
    end

    _axons = OrderedDict{Symbol,Vector{Axon}}()
    for (source_name, terminals) in pairs(axons)
        _terminals = getkey(_axons, source_name, Axon[])
        src = source_name ∈ inputs ? source_name : Ref(neurons[source_name])
        for (target_neuron_name, target_port_name, axon_kwargs...) in terminals
            tgt = target_neuron_name ∈ outputs ? target_neuron_name : (Ref(neurons[target_neuron_name]),target_port_name)
            syn = Axon{T}(src, tgt; axon_kwargs...)
            push!(_terminals, syn)
        end
        _axons[source_name] = _terminals
    end
    NeuralNetwork{T}(neurons, inputs, outputs, _axons; kwargs...)
end


text:
"""
    Molecules.symbol(A::Atom)

Returns a String representing the atomic symbol of the atom.
"""
function symbol(A::Atom)
    return elements[A.Z].symbol
end

"""
    Molecules.nuclear_repulsion(A::Atom, B::Atom)

Returns the repulsion energy between atoms A and B in Atomic units.
"""
function nuclear_repulsion(A::Atom, B::Atom)
    return (A.Z*B.Z)/(√((A.xyz.-B.xyz)⋅(A.xyz.-B.xyz))/bohr_to_angstrom)
end

"""
    Molecules.nuclear_repulsion(atoms::Molecule)

Returns the total nuclear respulsion energy of a group of atoms.
"""
function nuclear_repulsion(atoms::Molecule)
    N = length(atoms)
    E = 0.0
    for i in 1:N
        for j in (i+1):N
            E += nuclear_repulsion(atoms[i], atoms[j])
        end
    end
    return E
end

function ∇nuclear_repulsion(atoms::Molecule, i)
    E = zeros(3)

    A = atoms[i]

    for B in atoms
        B == A ? continue : nothing

        D = √((A.xyz.-B.xyz)⋅(A.xyz.-B.xyz))/bohr_to_angstrom

        E -= (B.Z / D^3) .* (A.xyz - B.xyz)/bohr_to_angstrom
    end

    return A.Z .* E ./ bohr_to_angstrom
end

"""
    Molecules.center_of_mass(atoms::Molecule)

Returns the center of mass of a group if atoms.
"""
function center_of_mass(atoms::Molecule)

    # Compute the mass-weighted XYZ
    cm_xyz = sum(a.xyz * a.mass for a in atoms)

    # Compute total mass
    M = sum(a.mass for a in atoms)

    # Return center of mass
    return cm_xyz ./ M
end


text:
using Test
using MPI

function check_for_query(comm)
    is_message, status = MPI.Iprobe(MPI.Consts.MPI_ANY_SOURCE[], MPI.Consts.MPI_ANY_TAG[], comm)
    if is_message
        recv_id = status.source
        tag_ind = status.tag
        return true, recv_id, tag_ind
    else
        return false, -1, -1
    end
end

MPI.Init()
comm = MPI.COMM_WORLD
myrank = MPI.Comm_rank(comm)
mysize = MPI.Comm_size(comm)

# First rank will send, one-at-a-time, the rank + 4 to each other rank.
#     They will then sum these and test if the sum is correct. Each rank > 0
#     will be waiting at the Ibarrier for messages from 0 in order to test
#     that the Ibarrier is working properly.
# They then communicate back to 0 when they have received all messages to
#     to allow it to reach the barrier.
#
# This is a contrived example but does test to ensure Ibarrier works

if myrank == 0
  for reps in 1:10
    for ii in 1:(mysize-1)
      smsg = MPI.Send(ii + 4, comm; dest=ii, tag=ii)
    end
  end
  for ii in 1:(mysize-1)
    dummy  = [0]
    rmsg   = MPI.Recv!(dummy, comm; source=ii, tag=ii)
  end
end

all_done = false
localsum = 0
msg_num = 0

barrier_req = MPI.Ibarrier(comm)

all_done = MPI.Test(barrier_req)

while !all_done
  global all_done
  global msg_num
  global myrank
  is_request, recv_id, tag_ind = check_for_query(comm)
  if is_request
    dummy    = [0]
    rmsg     = MPI.Recv!(dummy, comm; source=recv_id, tag=tag_ind)
    msg_num += 1
    global localsum += dummy[1]
    if msg_num == 10
      smsg = MPI.Send(tag_ind, comm; dest=0, tag=myrank)
    end
  end # is_request
  all_done = MPI.Test(barrier_req)
end # !all_done

if myrank > 0
  @test localsum == 10 * (myrank + 4)
end

MPI.Finalize()
@test MPI.Finalized()


text:
#function for setting up the linear machine model. Takes training
#  train::Array{StatisticalMoments{Float64}}         Training data set of moments, se StatisticalMoments.jl for definition of moments
#  fm::FeatureMap                Feature map  as defined in ML_models.jl
#  options                       Contains the regularization parameters and other optional arguments
#  kernalname                    Optional kernal name. kernelname = "" means no kernel and kernelname = "gauss" uses a gaussian kernel.
function setup_ml(train::Array{StatisticalMoments{Float64}}, fm::FeatureMap,kernelname::AbstractString,  options, output_type::AbstractString) #,
    # Get feature selection if there is one
    if (isempty(options.aux))
        options.aux = get_feature_L1(fm,output_type);
    end

    #Get all volfractions and particle radiuses
    volfracs= union([td.label[1] for td in train]);
    radiuses = union([td.label[2] for td in train]);
    numlabels = length(volfracs)*length(radiuses);
    #Initialize data matrix.  Each row is a row concatentation of all simulations of a single labels pair.
    featuresize = fm.get_numfeatures(train[1],options);
    X = zeros(numlabels,featuresize);
    ys =  zeros(numlabels,1);
    rowcounts = 1;
    # Initialize Machine Learning model
    ml = ML_model( 0, [0], 0.0, 0.0,[0], [0], [0.0], fm.name,featuresize,numlabels,[0], kernelname,0.0,x->x);
    ml.name = get_model_name(fm,kernelname,output_type,options);
#     if(addftL1toname) #If features were selected by L1, adding f to name
#         ml.name = string( ml.name, "-f")
#     end
    #Build data matrix.  Each row is a row concatentation of all simulations of a single labels pair.
    for mom in train
        #filling in each row of X with a feature vector
        fm.apply(X, rowcounts, mom,options);
        vl = mom.label[1] ;
        rd = mom.label[2] ;
        if(output_type=="radius")
              ys[rowcounts] = rd;
        elseif(output_type=="volumefraction")
              ys[rowcounts] = vl;
        else
              ys[rowcounts] = vl/rd^2;
        end
        rowcounts = rowcounts+1;
    end
    # Data transformations: Centering, scaling, exponential on positive terms
    X = X_transform(X, ml); # println("yvs before undo_transform: ",yvs)
    ys =y_transform(ys , ml);#    println("yvs after undo_transform: ",yvs)
    # Crop data
    s = sample(1:rowcounts-1, rowcounts-1,replace=false);
    X = X[s,:];  # clipping off  the zero rows and re-shuffling
    ys = ys[s];   # clipping off  the zero rows and re-shuffling

    ml.X = X;   # reshuffle the data at this point to avoid data leakage bias?
    ml.ys =ys;
    if(options.plotfigures)
        println("(#data, #features): (",size(ml.X)[1], ", ", size(ml.X)[2],")")
    end
    return ml;
end


text:
# This file is part of the StateMachines.jl Julia package
#
# Paul Bayer
# MIT License
# Github: https://github.com/pbayer/StateMachines.jl
#

"""
    isConsistent

Check if a system (including all of its subsystems) is consistent.
"""
function isConsistent(S::System)::Bool
    # a system must have registered blocks or processes
    # each of its input channels must be connected to a block or a process
    return false
end


text:
##############################################################################
#  function for conversion to image
##############################################################################

"""
Initializes `wkhtmltoimage` in graphics less mode `img_init(0)` or not
`img_init(1)`
"""
function img_init(usegraphics::Int)
  ccall((:wkhtmltoimage_init, libwkhtml),
        Int, (Int,), usegraphics)
end

"""
Closes `wkhtmltoimage`
"""
function img_deinit()
  ccall((:wkhtmltoimage_deinit, libwkhtml),
        Int, ())
end

### Global settings ###

"""
Sets the output parameters by passing pairs `"property" => "values"` (both should be String).

example:
```julia
ps = ImgSettings("out" => "/tmp/example.pdf",
                 "orientation" => "Landscape")
```

##Parameters :

- **crop.left** left/x coordinate of the window to capture in pixels. E.g. "200"
- **crop.top** top/y coordinate of the window to capture in pixels. E.g. "200"
- **crop.width** Width of the window to capture in pixels. E.g. "200"
- **crop.height** Height of the window to capture in pixels. E.g. "200"
- **load.cookieJar** Path of file used to load and store cookies.
- **load.* Page specific settings related to loading content, see Object Specific loading settings.
- **web.* See Web page specific settings.
- **transparent** When outputting a PNG or SVG, make the white background transparent. Must be either "true" or "false"
- **in** The URL or path of the input file, if "-" stdin is used. E.g. "http://google.com"
- **out** The path of the output file, if "-" stdout is used, if empty the content is stored to a internalBuffer.
- **fmt** The output format to use, must be either "", "jpg", "png", "bmp" or "svg".
- **screenWidth** The with of the screen used to render is pixels, e.g "800".
- **smartWidth** Should we expand the screenWidth if the content does not fit? must be either "true" or "false".
- **quality** The compression factor to use when outputting a JPEG image. E.g. "94".

##Web page specific settings

- **web.background** Should we print the background? Must be either "true" or "false".
- **web.loadImages** Should we load images? Must be either "true" or "false".
- **web.enableJavascript** Should we enable javascript? Must be either "true" or "false".
- **web.enableIntelligentShrinking** Should we enable intelligent shrinkng to fit more content on one page? Must be either "true" or "false". Has no effect for wkhtmltoimage.
- **web.minimumFontSize** The minimum font size allowed. E.g. "9"
- **web.printMediaType** Should the content be printed using the print media type instead of the screen media type. Must be either "true" or "false". Has no effect for wkhtmltoimage.
- **web.defaultEncoding** What encoding should we guess content is using if they do not specify it properly? E.g. "utf-8"
- **web.userStyleSheet** Url er path to a user specified style sheet.
- **web.enablePlugins** Should we enable NS plugins, must be either "true" or "false". Enabling this will have limited success.

##Object Specific loading settings

- **load.username** The user name to use when loging into a website, E.g. "bart"
- **load.password** The password to used when logging into a website, E.g. "elbarto"
- **load.jsdelay** The mount of time in milliseconds to wait after a page has done loading until it is actually printed. E.g. "1200". We will wait this amount of time or until, javascript calls window.print().
- **load.zoomFactor** How much should we zoom in on the content? E.g. "2.2".
- **load.customHeaders** TODO
- **load.repertCustomHeaders** Should the custom headers be sent all elements loaded instead of only the main page? Must be either "true" or "false".
- **load.cookies** TODO
- **load.post** TODO
- **load.blockLocalFileAccess** Disallow local and piped files to access other local files. Must be either "true" or "false".
- **load.stopSlowScript** Stop slow running javascript. Must be either "true" or "false".
- **load.debugJavascript** Forward javascript warnings and errors to the warning callback. Must be either "true" or "false".
- **load.loadErrorHandling** How should we handle obejcts that fail to load. Must be one of:
  - **"abort"** Abort the convertion process
  - **"skip"** Do not add the object to the final output
  - **"ignore"** Try to add the object to the final output.
- **load.proxy** String describing what proxy to use when loading the object.
- **load.runScript** TODO

##Header and footer settings

The same settings can be applied for headers and footers, here there are explained in terms of the header.

- **header.fontSize** The font size to use for the header, e.g. "13"
- **header.fontName** The name of the font to use for the header. e.g. "times"
- **header.left** The string to print in the left part of the header, note that some sequences are replaced in this string, see the wkhtmltopdf manual.
- **header.center** The text to print in the center part of the header.
- **header.right** The text to print in the right part of the header.
- **header.line** Whether a line should be printed under the header (either "true" or "false").
- **header.spacing** The amount of space to put between the header and the content, e.g. "1.8". Be aware that if this is too large the header will be printed outside the pdf document. This can be corrected with the margin.top setting.
- **header.htmlUrl** Url for a HTML document to use for the header.

"""
type ImgSettings
  ptr::Ptr{GlobalSettings}
end

# 0 arg case
ImgSettings() = ImgSettings(img_create_global_settings())

# >= 1 case
function ImgSettings(p1::Pair{String,String}, args...)
  ptr = img_create_global_settings()

  img_set_global_setting(ptr, p1.first, p1.second)
  for p in args
    isa(p, Pair{String,String}) || error("argument $p is not Pair{String,String}")
    img_set_global_setting(ptr, p.first, p.second)
  end

  ImgSettings(ptr)
end

function setindex!(settings::ImgSettings, value::String, property::String)
  img_set_global_setting(settings.ptr, property, value)
end

function getindex(settings::ImgSettings, property::String)
  img_get_global_setting(settings.ptr, property)
end

function img_create_global_settings()
  ccall((:wkhtmltoimage_create_global_settings, libwkhtml),
        Ptr{GlobalSettings}, ())
end

function img_get_global_setting(gs::Ptr{GlobalSettings}, property::String)
  cval = zeros(UInt8, 100)
  ccall((:wkhtmltoimage_get_global_setting, libwkhtml),
        Int,
        (Ptr{GlobalSettings}, Cstring, Cstring, Int),
        gs, convert(Cstring, pointer(property)),
            convert(Cstring, pointer(cval)), sizeof(cval)-1 )
  unsafe_string(pointer(cval))
end

function img_set_global_setting(gs::Ptr{GlobalSettings}, property::String, value::String)
  ccall((:wkhtmltoimage_set_global_setting, libwkhtml),
        Int,
        (Ptr{GlobalSettings}, Cstring, Cstring),
        gs,
        convert(Cstring, pointer(property)),
        convert(Cstring, pointer(value)) )
end

### conversion ###

type ImgConverter
  ptr::Ptr{Converter}
end

function ImgConverter(settings::ImgSettings)
  obj = ImgConverter(img_create_converter(settings.ptr))

  finalizer(obj, o -> img_destroy_converter(o.ptr))

  obj
end

function run(conv::ImgConverter)
  img_wkconvert(conv.ptr)
end

function img_create_converter(gs::Ptr{GlobalSettings})
  ccall((:wkhtmltoimage_create_converter, libwkhtml),
        Ptr{Converter},
        (Ptr{GlobalSettings}, Cstring),
        gs, C_NULL )
end

function img_destroy_converter(conv::Ptr{Converter})
  ccall((:wkhtmltoimage_destroy_converter, libwkhtml),
        Void,
        (Ptr{Converter},),
        conv )
end

function img_wkconvert(conv::Ptr{Converter})
  ccall((:wkhtmltoimage_convert, libwkhtml),
        Int, (Ptr{Converter},), conv )
end

function img_add_object(conv::Ptr{Converter}, os::Ptr{ObjectSettings},
                    data::Union{String, Ptr{Void}})
  ccall((:wkhtmltoimage_add_object, libwkhtml),
        Void,
        (Ptr{Converter}, Ptr{ObjectSettings}, Cstring),
        conv, os,
        isa(data, String) ? convert(Cstring, pointer(data)) : data)
end


text:
using BenchmarkTools
include("ZZDiffusionBridge.jl")

#comment out runall() in example_sin.jl
include("./examples/example_sin.jl")

# see notes section: tune velocities
# the prefactor (β + 1.0) is just made to make a
# bit more comparable the experiments

function tune_velocity( L::Int64, β = 0.0)
    θ = fill(0.0, 2<<L - 1)
    k = 1
    for i in 0:L
        for j in 0:2^i-1
            θ[k] =  (3*β + 0.5)*2.0^(-i*β)
            k += 1
        end
    end
    return(θ)
end

#just little differences from the original zz_sampler (src/zz_sampler).
function zz_sampler(S::System, X::AbstractModel, T::Float64, L::Int64, u::Float64, v::Float64, clock::Float64, τ0::Float64 , n0::Int64)
    t = 0.0
    Ξ = Skeleton[]
    while t < clock
        τ0 , n0 = findmin(S.τ)
        S.ξ .+=  S.θ*τ0
        t += τ0
        if acc_rej(n0, S, X, sampling_scheme(X), u, v, t)
            S.θ[n0] = -S.θ[n0]
            push!(Ξ, (Skeleton(copy(S.ξ), t)))
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), true, u, v, t)
        else
            update_events!(n0, τ0, S, X, sampling_scheme(X), dependence_strucute(X), false, u, v, t)
        end
    end
    Ξ, τ0, n0
end

# the the average of the path over one realization.
# since the path is piecewise linear. We take the midpoint.
function piec_lin_integration(y, L, T)
    N = 2<<L
    somma = 0.5*(y[1] + y[end]) + sum(y[2:(end-1)])
    return somma*T/N
end

#given an array of skeleton, compute the sample avergae
# `sampling _frequency`. Discretization given the continuous path
# L truncation, u, v initial and final points
function compute_average(X::Array{Skeleton,1}, sampling_frequncy::Float64, u, v, L, T)
    dt = 0.5:sampling_frequncy:T-0.5 #if you start at 0.0 and finish at T the function Interpolate gives errors
    M = 0
    count = 0
    for i in dt
        ξ = FindCoordinates(X, i).ξ
        ξfe = fs_expansion(ξ, u, v, L, T)
        path = piec_lin_integration(ξfe, L, T)
        M = M + path
        count += 1
    end
    return M/count
end

# experiment: timing each ZigZag of length `clock`. repeat it
# `rep` times. return the an array with the single repretition and the time
# β determines the decay of the velocities
# L truncation level

function vel_experiment(β, L, clock, rep)
    α = 0.5
    T = 200.0
    u = -3π
    v = 3π
    τ0 = 0.0
    n0 = 0
    ξ = fill(0.0, 2<<L - 1)
    θ = tune_velocity(L, β)
    S = System(L, T, ξ, θ)
    X = SinSDE(α)
    Y = []
    TIME = []
    tot_time = 0
    for i in 1:rep
        (Ξ, τ0, n0), time = @timed  zz_sampler(S, X, T, L, u, v, clock, τ0, n0)
        #S::System, X::AbstractModel, T::Float64, L::Int64, u, v, clock, τ0 , n0)
        tot_time += time
        push!(Y, copy(Ξ) )
        push!(TIME, tot_time )
    end
    return Y, TIME
end

function runall_velocity()
    rep = 20
    clock = 100
    Y = fill(0.0,(6,rep))
    Time = fill(0.0,(6,rep))
    row = 1
    L = 10
    p = plot()
    for β in 0.0:0.1:0.5
        X, time = vel_experiment(β, L, 100.0, rep)
        ave = cumsum(compute_average.(X, 1.0, -3π,  3π, L, T))
        Y[row,:] =  [ave[i]/i for i in 1:rep]
        Time[row, :] = time
        plot!(p, Time[row,:], abs.(Y[row,:]), label =("beta =  $β"))
        row += 1
    end
    display(p)
    return Y, Time
end

X, T = runall_velocity()

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
"""
This struct contains the technical data of a renewable generator.
    type: Generator technology type.
    active_power_limits: Named tuple of minimum and maximum generation limits.
    ramp_limits: Ramp limits (MW/hr).
    operation_cost: Operation cost based on PowerSystems.jl definitions.
    bus: Name of the bus on which the generator is located.
    zone: Name of zone in which the generator is located.
"""
struct RenewableTech
    type::String
    active_power_limits::NamedTuple{(:min, :max),Tuple{Float64,Float64}}
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    operation_cost::Union{Nothing, PSY.TwoPartCost}
    bus::Int64
    zone::String
end

get_type(tech::RenewableTech) = tech.type
get_active_power_limits(tech::RenewableTech) = tech.active_power_limits
get_ramp_limits(tech::RenewableTech) = tech.ramp_limits
get_operation_cost(tech::RenewableTech) = tech.operation_cost
get_bus(tech::RenewableTech) = tech.bus
get_zone(tech::RenewableTech) = tech.zone

"""
This struct contains all the data for a renewable generator.
    name: Generator's name.
    tech: Struct containing technical data.
    decision_year: The year in which the investment decision is made.
    construction_year: The year in which construction is completed.
    retirement_year: The year in which the generator is retired.
    end_life_year: Age-based retirement year.
    products: Vector containing the data of all products provided by the generator.
    finance_data: Struct containing financial data.
"""
mutable struct RenewableGenEMIS{T <: BuildPhase} <: GeneratorEMIS{T}
    name::String
    tech::RenewableTech
    decision_year::Int64
    construction_year::Int64
    retirement_year::Int64
    end_life_year::Int64
    products::Vector{Product}
    finance_data::Finance
end


text:
#comparison.jl -- comparison operators for valids.

function Base.:(==)(x::Valid{N,ES}, y::Valid{N,ES}) where {N,ES}
    if (next(x.upper) == x.lower) & (next(y.upper) == y.lower)

        x_null = (x.upper == zero(Vnum{N,ES})) | (x.lower == zero(Vnum{N,ES}))
        y_null = (y.upper == zero(Vnum{N,ES})) | (y.lower == zero(Vnum{N,ES}))

        return x_null == y_null
    end
    return (x.upper == y.upper) && (x.lower == y.lower)
end

function Base.max(x::Sigmoid{N,ES,:exact}, y::Sigmoid{N,ES,:lower}) where {N,ES}
    x == Sigmoid{N,ES,:exact}(Inf) && return x
    (@s x) <= (@s y) ?  y : x
end
function Base.max(x::Sigmoid{N,ES,:exact}, y::Sigmoid{N,ES,:upper}) where {N,ES}
    x == Sigmoid{N,ES,:exact}(Inf) && return x
    (@s x) >= (@s y) ?  x : y
end
Base.max(x::Sigmoid{N,ES,:lower}, y::Sigmoid{N,ES,:exact}) where {N,ES} = max(y,x)
Base.max(x::Sigmoid{N,ES,:upper}, y::Sigmoid{N,ES,:exact}) where {N,ES} = max(y,x)

function Base.min(x::Sigmoid{N,ES,:exact}, y::Sigmoid{N,ES,:lower}) where {N,ES}
    x == Sigmoid{N,ES,:exact}(Inf) && return x
    (@s x) <= (@s y) ?  x : y
end
function Base.min(x::Sigmoid{N,ES,:exact}, y::Sigmoid{N,ES,:upper}) where {N,ES}
    x == Sigmoid{N,ES,:exact}(Inf) && return x
    (@u s) >= (@s y) ?  y : x
end
Base.min(x::Sigmoid{N,ES,:lower}, y::Sigmoid{N,ES,:exact}) where {N,ES} = min(y,x)
Base.min(x::Sigmoid{N,ES,:upper}, y::Sigmoid{N,ES,:exact}) where {N,ES} = min(y,x)

function Base.abs(x::Valid{N,ES}) where {N,ES}
    if roundsinf(x)
        return containszero(x) ? zero(Vnum{N,ES}) → Vnum{N,ES}(Inf) : min(x.lower, -x.upper) → Vnum{N,ES}(Inf)
    elseif containszero(x)
        return zero(Vnum{N,ES}) → max(-x.lower, x.upper)
    else
        return nonpositive(x) ? -x : x
    end
end

#in order to get lu factorization working we need to implement a very special
#form of the > function.

Base.isless(x::Valid{N,ES}, y::Valid{N,ES}) where {N,ES} = x.upper < y.upper
Base.:<(x::Valid{N,ES}, y::Valid{N,ES}) where {N,ES} = Base.isless(x,y)

#not-nowhere-equal
function nowhere_equal(a::Valid{N,ES}, b::Valid{N,ES}) where {N,ES}
    if roundsinf(a)
        roundsinf(b) && return false
        (a.lower > b.upper) && (a.upper < b.lower)
    elseif roundsinf(b)
        (a.lower > b.upper) && (a.upper < b.lower)
    else
        (a.lower > b.upper) || (b.lower > a.upper)
    end
end

≸(a::Valid{N,ES}, b::Valid{N,ES}) where {N,ES} = !nowhere_equal(a,b)
≹(a::Valid{N,ES}, b::Valid{N,ES}) where {N,ES} = ≸(a, b)
#≮(a::Valid{N,ES}, b::Valid{N,ES}) = !(a < b)
#≯(a::Valid{N,ES}, b::Valid{N,ES}) = !(a > b)

export ≸,≹,≯,≮

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
#!/usr/bin/julia

f = open("2019/03/input.txt", "r")
inputs = [line for line in readlines(f)]
close(f)

#println(inputs)

function manhattan(c1, c2)
    return abs(c1[1] - c2[1]) + abs(c1[2] - c2[2])
end

function manhattan_from_zero(c)
    return abs(c[1]) + abs(c[2])
end

function collision(a, b, x, y)
    collide = false
    # ab horizontal and xy vertical
    if a[2] == b[2] && x[1] == y[1]
        #println("ab horizontal and xy vertical")
        min_x = min(a[1], b[1])
        max_x = max(a[1], b[1])
        min_y = min(x[2], y[2])
        max_y = max(x[2], y[2])
        if min_x < x[1] < max_x && min_y < a[2] < max_y
            collide = true
        end
    # ab vertical and xy horizontal
    elseif a[1] == b[1] && x[2] == y[2]
        #println("ab vertical and xy horizontal")
        min_x = min(x[1], y[1])
        max_x = max(x[1], y[1])
        min_y = min(a[2], b[2])
        max_y = max(a[2], b[2])
        if min_x < a[1] < max_x && min_y < x[2] < max_y
            collide = true
        end
    end
    return collide
end

function get_col_coords(a, b, x, y)
    # if ab horizontal and xy vertical
    if a[2] == b[2] && x[1] == y[1]
        return [x[1], a[2]]
    # if ab vertical and xy horizontal
    elseif a[1] == b[1] && x[2] == y[2]
        return[a[1], x[2]]
    end
end

function get_min_dist(coords)
    min_dist = typemax(Int)
    for c in coords
        m = manhattan_from_zero(c)
        if m < min_dist
            min_dist = m
        end
    end
    return min_dist
end

function find_col_coords(wire_points)
    collision_coords = []
    for z in collect(2:length(wire_points))
        w1 = wire_points[z-1]
        w2 = wire_points[z]
        for i in collect(2:length(w1))
            w1_1 = w1[i-1]
            w1_2 = w1[i]
            for j in collect(2:length(w2))
                w2_1 = w2[j-1]
                w2_2 = w2[j]
                #println(w1_1, ",", w1_2, " vs ", w2_1, ",", w2_2)
                if collision(w1_1, w1_2, w2_1, w2_2)
                    colcords = get_col_coords(w1_1, w1_2, w2_1, w2_2)
                    #println("Collide! At the disco ", colcords, ": ", w1_1, ",", w1_2, " ", w2_1, ",", w2_2)
                    push!(collision_coords, colcords)
                end
            end
        end
    end
    return collision_coords
end

function get_wire_points(inputs)
    wire_points = [[[0, 0]] for x in inputs]
    for i in eachindex(inputs)
        x, y = 0, 0
        for move in split(inputs[i], ",")
            dir = first(move)
            dist = parse(Int, move[2:end])
            if dir == 'D'
                y -= dist
            elseif dir == 'U'
                y += dist
            elseif dir == 'R'
                x += dist
            elseif dir == 'L'
                x -= dist
            end
            push!(wire_points[i], [x, y])
        end
    end
    return wire_points
end

function get_delays(wire_points, collision_coords)
    delays = []
    for wire in wire_points
        wire_delays = [typemax(Int) for x in collision_coords]
        wire_steps = 0
        position = [0,0]
        for i in eachindex(wire)[2:end]
            dest = wire[i]
            for z in collect(1:manhattan(position, dest))
                if position in collision_coords
                    # check for known number of steps
                    ind = findfirst(isequal(position), collision_coords)
                    if wire_steps < wire_delays[ind]
                        wire_delays[ind] = wire_steps
                    end
                end

                if position == dest
                    # destination reached
                    println(position, " reached, steps=", wire_steps+1)
                elseif position[1] == dest[1]
                    # moving in y dir
                    if position[2] < dest[2]
                        position[2] += 1
                    elseif position[2] > dest[2]
                        position[2] -= 1
                    else
                        println("pos = dest check broken")
                    end
                elseif position[2] == dest[2]
                    # moving in x dir
                    if position[1] < dest[1]
                        position[1] += 1
                    elseif position[1] > dest[1]
                        position[1] -= 1
                    else
                        println("pos = dest check broken v2")
                    end
                else
                    println("uh on spaghetti o no no no")
                end
                wire_steps += 1
            end
        end
        #println(wire_steps)
        push!(delays, wire_delays)
    end
    return delays
end

function get_min_delay(delays)
    combined = delays[1] + delays[2]
    return min(combined...)
end

wire_points = get_wire_points(inputs)

#println(wire_points)

collision_coords = find_col_coords(wire_points)
#println(collision_coords)

println("Min manhattan distance: ", get_min_dist(collision_coords))

d = get_delays(wire_points, collision_coords)
#println(d)

println("Min delay distance: ", get_min_delay(d))

max_stars_repo_name: chaserobertson/advent

max_issues_repo_name: chaserobertson/advent

max_forks_repo_name: chaserobertson/advent


text:
#############################################################################
# entropy.jl
# entropy (ie, sum_i( -x_i log (x_i) ) of an expression x
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

export entropy
export sign, curvature, monotonicity, evaluate

### Entropy: sum_i -x_i log (x_i)

# TODO: make this work for a *list* of inputs, rather than just for scalar/vector/matrix inputs

# Entropy atom: -xlogx entrywise
mutable struct EntropyAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function EntropyAtom(x::AbstractExpr)
        if sign(x) == ComplexSign()
            error("The argument should be real but it's instead complex")
        else
            children = (x,)
            # TODO check positivity or enforce it
            return new(:entropy, hash(children), children, size(x))
        end
    end
end

function sign(x::EntropyAtom)
    return NoSign()
end

function monotonicity(x::EntropyAtom)
    return (NoMonotonicity(),)
end

function curvature(x::EntropyAtom)
    return ConcaveVexity()
end

function evaluate(x::EntropyAtom)
    c = evaluate(x.children[1])
    return -c .* log.(c)
end

function conic_form!(e::EntropyAtom, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, e)
        # -x log x >= t  <=>  x exp(t/x) <= 1  <==>  (t,x,1) in exp cone
        t = Variable(e.size)
        x = e.children[1]
        objective = conic_form!(t, unique_conic_forms)
        for i = 1:size(x, 1)
            for j = 1:size(x, 2)
                conic_form!(ExpConstraint(t[i, j], x[i, j], 1), unique_conic_forms)
            end
        end
        cache_conic_form!(unique_conic_forms, e, objective)
    end
    return get_conic_form(unique_conic_forms, e)
end

entropy(x::AbstractExpr) = sum(EntropyAtom(x))


text:
include("types.jl")
include("events.jl")
include("requests.jl")
include("responses.jl")
include("dap_message_defs.jl")


text:
module IntoTheJungleCodeModRemovePlatformVerticalCollisionController

using ..Ahorn, Maple

@mapdef Entity "IntoTheJungleCodeMod/RemovePlatformVerticalCollisionController" RemovePlatformVerticalCollisionController(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Remove Platform Vertical Collision Controller\n(Into The Jungle code mod)" => Ahorn.EntityPlacement(
        RemovePlatformVerticalCollisionController
    )
)

function Ahorn.selection(entity::RemovePlatformVerticalCollisionController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::RemovePlatformVerticalCollisionController, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.northernLights, -12, -12)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper