text:
# Indexing at a point
@inline function (itp::GriddedInterpolation{T,N})(x::Vararg{Number,N}) where {T,N}
    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
    wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
    coefficients(itp)[wis...]
end
@propagate_inbounds function (itp::GriddedInterpolation{T,N})(x::Vararg{Number,M}) where {T,M,N}
    inds, trailing = split_trailing(itp, x)
    @boundscheck (check1(trailing) || Base.throw_boundserror(itp, x))
    @assert length(inds) == N
    itp(inds...)
end
@inline function (itp::GriddedInterpolation)(x::Vararg{UnexpandedIndexTypes})
    xis = to_indices(itp, x)
    xis == x && error("evaluation not supported for GriddedInterpolation at positions $x")
    itp(xis...)
end

@inline function gradient(itp::GriddedInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
    wis = weightedindexes((value_weights, gradient_weights), itpinfo(itp)..., x)
    SVector(map(inds->coefficients(itp)[inds...], wis))
end

itpinfo(itp::GriddedInterpolation) = (tcollect(itpflag, itp), itp.knots)

roundbounds(x::Integer, knotvec::AbstractVector) = gridded_roundbounds(x, knotvec)
roundbounds(x::Number, knotvec::AbstractVector) = gridded_roundbounds(x, knotvec)
function gridded_roundbounds(x, knotvec::AbstractVector)
    i = find_knot_index(knotvec, x)
    iclamp = max(i, first(axes1(knotvec)))
    inext = min(iclamp+1, last(axes1(knotvec)))
    ifelse(i < iclamp, i+1, ifelse(x - knotvec[iclamp] < knotvec[inext] - x, i, inext))
end

floorbounds(x::Integer, knotvec::AbstractVector) = gridded_floorbounds(x, knotvec)
floorbounds(x, knotvec::AbstractVector) = gridded_floorbounds(x, knotvec)
function gridded_floorbounds(x, knotvec::AbstractVector)
    i = find_knot_index(knotvec, x)
    max(i, first(axes1(knotvec)))
end

ceilbounds(x::Integer, knotvec::AbstractVector) = gridded_ceilbounds(x, knotvec)
ceilbounds(x, knotvec::AbstractVector) = gridded_ceilbounds(x, knotvec)
function gridded_ceilbounds(x, knotvec::AbstractVector)
    i = find_knot_index(knotvec, x)
    iclamp = max(i, first(axes1(knotvec)))
    min(iclamp+1, last(axes1(knotvec)))
end

@inline find_knot_index(knotv, x) = searchsortedfirst(knotv, x, Base.Order.ForwardOrdering()) - 1
@inline find_knot_index(knotv, x::AbstractVector) = searchsortedfirst_vec(knotv, x) .- 1

@inline function weightedindex_parts(fs::F, mode::Gridded, knotvec::AbstractVector, x) where F
    i = find_knot_index(knotvec, x)
    weightedindex_parts2(fs, mode, knotvec, x, i)
end

@inline function weightedindex_parts2(fs::F, mode::Gridded, knotvec::AbstractVector, x, i) where F
    ax1 = axes1(knotvec)
    iclamp = clamp(i, first(ax1), last(ax1)-1)
    weightedindex(fs, degree(mode), knotvec, x, iclamp)
end

function weightedindex(fs::F, deg::Constant, knotvec, x, iclamp) where F
    pos, δx = positions(deg, knotvec,  x)
    (position=pos, coefs=fmap(fs, deg, δx))
end
function weightedindex(fs::F, deg::Degree, knotvec, x, iclamp) where F
    @inbounds l, u = knotvec[iclamp], knotvec[iclamp+1]
    δx = ratio(x - l, u - l)
    (position=iclamp, coefs=rescale_gridded(fs, fmap(fs, deg, δx), u-l))
end

rescale_gridded(fs::F, coefs, Δx) where F =
    (rescale_gridded(fs[1], coefs[1], Δx), rescale_gridded(Base.tail(fs), Base.tail(coefs), Δx)...)
rescale_gridded(::Tuple{}, ::Tuple{}, Δx) = ()
rescale_gridded(::typeof(value_weights), coefs, Δx) = coefs
rescale_gridded(::typeof(gradient_weights), coefs, Δx) = coefs./Δx
rescale_gridded(::typeof(hessian_weights), coefs, Δx) = coefs./Δx.^2

@inline function (itp::GriddedInterpolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
    itps = tcollect(itpflag, itp)
    if x[1] isa AbstractVector
        wis = dimension_wis_vec(value_weights, itps, itp.knots, x)
    else
        wis = dimension_wis(value_weights, itps, itp.knots, x)
    end
    coefs = coefficients(itp)
    ret = [coefs[i...] for i in Iterators.product(wis...)]
    reshape(ret, shape(wis...))
end

function dimension_wis(f::F, itps, knots, xs) where F
    itpflag, knotvec, x = itps[1], knots[1], xs[1]
    function makewi(y)
        pos, coefs = weightedindex_parts((f,), itpflag, knotvec, y)
        maybe_weightedindex(pos, coefs[1])
    end
    (makewi.(x), dimension_wis(f, Base.tail(itps), Base.tail(knots), Base.tail(xs))...)
end

function dimension_wis_vec(f::F, itps, knots, xs) where F
    itpflag, knotvec, x = itps[1], knots[1], xs[1]
    ivec = find_knot_index(knotvec, x)
    function makewi(y, i)
        pos, coefs = weightedindex_parts2((f,), itpflag, knotvec, y, i)
        maybe_weightedindex(pos, coefs[1])
    end
    (makewi.(x, ivec), dimension_wis(f, Base.tail(itps), Base.tail(knots), Base.tail(xs))...)
end

function dimension_wis(f::F, itps::Tuple{NoInterp,Vararg{Any}}, knots, xs) where F
    (Int.(xs[1]), dimension_wis(f, Base.tail(itps), Base.tail(knots), Base.tail(xs))...)
end
dimension_wis(f, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()

function getindex_return_type(::Type{GriddedInterpolation{T,N,TCoefs,IT,K}}, argtypes) where {T,N,TCoefs,IT<:DimSpec{Gridded},K}
    Tret = TCoefs
    for a in argtypes
        Tret = Base.promote_op(*, Tret, a)
    end
    Tret
end

Base.@propagate_inbounds function searchsortedfirst_exp_left(v, xx, lo, hi)
    for i in 0:4
        ind = lo + i
        ind > hi && return ind
        xx <= v[ind] && return ind
    end
    n = 3
    tn2 = 2^n
    tn2m1 = 2^(n-1)
    ind = lo + tn2
    while ind <= hi
        xx <= v[ind] && return searchsortedfirst(v, xx, lo + tn2 - tn2m1, ind, Base.Order.Forward)
        tn2 *= 2
        tn2m1 *= 2
        ind = lo + tn2
    end
    return searchsortedfirst(v, xx, lo + tn2 - tn2m1, hi, Base.Order.Forward)
end

function searchsortedfirst_vec(v::AbstractVector, x::AbstractVector)
    issorted(x) || return searchsortedfirst.(Ref(v), x)
    out = zeros(Int, length(x))
    lo = 1
    hi = length(v)
    @inbounds for i in 1:length(x)
        xx = x[i]
        y = searchsortedfirst_exp_left(v, xx, lo, hi)
        out[i] = y
        lo = min(y, hi)
    end
    return out
end


text:
using Distributions
using StatUtils
using Test
using StableRNGs

function test_bootstrapper()
    RNG = StableRNG(1234)
    x = rand(RNG, 1000)
    B = StatUtils.run_bootstrap(1000,sum,x;RNG=RNG)
    @test B.μ ≈ 508.1411134730461
    @test B.σ ≈ 9.045400206380133
    println("Bootstrapper test passed")

end

function test_fdr()
	pvalues = [0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.00]
	pidx = StatUtils.fdr(pvalues)
	@test pidx == [1,2,3,4]
	println("FDR test passed")
end

function test_groupby()
	groups = StatUtils.groupby([1,2,3], [1,2,3])
	@test Dict{Int64, Array{Int64,1}}(1=>[1], 2=>[2], 3=>[3]) == groups
	println("groupby test passed")
end

@testset "Bootstrapper" begin
    test_bootstrapper()
end
@testset "False Discovery Rate" begin
    test_fdr()
end
@testset "Group by" begin
    test_groupby()
end

@testset "Bootstrap median" begin
    RNG = StableRNG(1234)
    x1 = rand(RNG, Float64, (1,100))
    μ,σ = StatUtils.bootstrap_median(x1, 10_000, RNG)
    @test μ[1] ≈ 0.4977072659335834
    @test σ[1] ≈ 0.04304170028498684
end

@testset "Spearman pv" begin
    #no correlation
    RNG = StableRNG(1234)
    x1 = rand(RNG, 100)
    x2 = rand(RNG, 100)
    cc,pv = StatUtils.spearmanr(x1,x2, :both, RNG)
    @test cc ≈ -0.023438343834383438
    @test pv ≈ 0.999

    #small non-significant correlation
    RNG = StableRNG(1234)
    x1 = rand(RNG, 100)
    x2 = 0.1*x1 +  0.9*rand(RNG, 100)
    cc,pv = StatUtils.spearmanr(x1,x2, :right, RNG)
    @test cc ≈ 0.09610561056105611
    @test pv ≈ 0.189
    #strong significant correlation
    RNG = StableRNG(1234)
    x1 = rand(RNG, 100)
    x2 = 0.5*x1 +  0.5*rand(RNG, 100)
    cc,pv = StatUtils.spearmanr(x1,x2, :right, RNG)
    @test cc ≈ 0.6623702370237023
    @test pv ≈ 0.0
end

@testset "Percentile" begin
    @test_throws ArgumentError Percentile(-0.1)
    @test_throws ArgumentError Percentile(101.0)
    @test typemax(Percentile) == 100.0
    @test typemin(Percentile) == 0.0
end

@testset "Generalized Poisson" begin
    G = GeneralizedPoisson(1.5370664608871476,-0.2658194383776511)
    @test pdf(G,0) ≈ 0.21501092012028916
    @test log(pdf(G,0)) ≈ logpdf(G,0)
    @test log(pdf(G,1)) ≈ logpdf(G,1)
    @test mean(G) ≈ 1.2142857142857142
    @test var(G) ≈ 0.7578397212543557
    @test sum(pdf.(G,0:7)) ≈ 1.0000005196278676
end

@testset "Bootstrap regression" begin
    #generate data
    RNG = StableRNG(1234)
    x = range(0.0, stop=1.0, length=50)
    y = 0.5 .+ 0.5*x .+ 0.05*randn(RNG, length(x))
    μ, σ, xx = StatUtils.bootstrap_regression(x,y;RNG=RNG)
    sse = sum(abs2, μ-y)
    sq = sum(abs2, y .- 0.5)
    @test sse/sq < 0.05
end

@testset "Robust regression" begin
    RNG = StableRNG(1234)
    x = range(0.0, stop=1.0, length=20)
    y = 0.5 .+ 0.1*x .+ 0.3*randn(RNG, length(x))
    β0 = rand(RNG, Float64, 2)
    β = StatUtils.robust_regression(x,y, β0)
    @test β[1] ≈ 0.23571370032430217
    @test β[2] ≈ 0.38013900282542895
end

@testset "Angular histogram" begin
    θ = [0.0, π/4, π/2, -π, -π/4, -π/2, -3*π/4]
    counts, binidx, bins = StatUtils.angular_histogram(θ, π/4)
    @test bins ≈ [-π, -3*π/4, -π/2, -π/4, 0.0, π/4, π/2, 3π/4]
    @test counts == [1, 1, 1, 1, 1, 1, 1, 0]
    @test binidx == [5, 6, 7, 1, 4, 3, 2]
end

@testset "Grouping" begin
    x = [1,3,3,1,5,7,7,1,5]
    y = StatUtils.compress(x)
    @test y == [1,2,2,1,3,4,4,1,3]
    X = fill(0.0, 3,4,10,6)
    grouping = fill(0, 10)
    grouping[1:5] .= 1
    grouping[6:10] .= 2
    X[:,:,findall(grouping.==1),:] .= reshape([1.0,2.0,3.0,4.0,5.0], 1,1,5,1)
    X[:,:,findall(grouping.==2),:] .= reshape([6.0, 7.0, 8.0, 9.0, 10.0], 1, 1, 5,1)
    sidx = [7, 10, 2, 5, 8, 6, 1, 3, 9, 4]
    X = X[:,:,sidx,:]
    μ = mean(X, grouping[sidx])
    @test size(μ) == (3,4,2,6)
    @test μ[:,:,1,:] ≈ fill(3.0, 3,4,6)
    @test μ[:,:,2,:] ≈ fill(8.0, 3,4,6)
end

ext: jl


text:
using TensorOperations

function propagate(LHS, X, Y)
    @tensor R[6,7,8] := LHS[1,2,3]*X[1,4,6]*Y[2,5,7,4]*X[3,5,8]
end


text:
DE = de_rand_1_bin({
  :SearchSpace => symmetric_search_space(1, (0.0, 10.0)),
  :Population => reshape([1.0:10.0], 10, 1),
  "f" => 0.4,
  "cr" => 0.5,
  "NumParents" => 3})

facts("Differential evolution optimizer") do

context("random_sampler") do
  for(i in 1:NumTestRepetitions)
    numSamples = rand(1:8)
    sampled = BlackBoxOptim.random_sampler(DE, numSamples)

    @fact length(sampled) => numSamples

    popindices = 1:popsize(DE)

    # All sampled indices are indices into the population
    @fact all([in(index, popindices) for index in sampled]) => true
  end
end

context("radius_limited_sampler") do
  DE = de_rand_1_bin({
    :SearchSpace => symmetric_search_space(1, (0.0, 10.0)),
    :Population => rand(100,1),
    "f" => 0.4, "cr" => 0.5, "NumParents" => 3})

  psize = popsize(DE)
  popindices = 1:popsize(DE)
  sampler_radius = DE.options["SamplerRadius"]

  for(i in 1:NumTestRepetitions)
    numSamples = rand(1:sampler_radius)
    sampled = BlackBoxOptim.radius_limited_sampler(DE, numSamples)

    @fact length(sampled) => numSamples

    # All sampled indices are indices into the population
    @fact all([in(index, popindices) for index in sampled]) => true

    mini, maxi = minimum(sampled), maximum(sampled)
    if (maxi - mini) > max(numSamples+2, sampler_radius)
      @fact mini + psize <= maxi + sampler_radius => true
    end
  end
end

context("rand_bound_from_target!") do
  context("does nothing if within bounds") do
    @fact BlackBoxOptim.rand_bound_from_target!([0.0], [0.0], [(0.0, 1.0)]) => [0.0]

    @fact BlackBoxOptim.rand_bound_from_target!([0.0, 11.4], [0.1, 12.3],
      RangePerDimSearchSpace([(0.0, 1.0), (10.0, 15.0)])) => [0.0, 11.4]
  end

  context("bounds if lower than min bound") do
    @fact BlackBoxOptim.rand_bound_from_target!([-0.1], [0.0], [(0.0, 1.0)]) => [0.0]

    res = BlackBoxOptim.rand_bound_from_target!([-0.1], [0.5], [(0.0, 1.0)])
    @fact 0.5 >= res[1] >= 0.0 => true

    res = BlackBoxOptim.rand_bound_from_target!([-11.1, 0.5], [-10.8, 0.5], [(-11.0, 1.0), (0.0, 1.0)])
    @fact -10.8 >= res[1] >= -11.0 => true
    @fact res[2] => 0.5

    res = BlackBoxOptim.rand_bound_from_target!([50.4, -103.1], [49.6, -101.4], [(30.0, 60.0), (-102.0, -1.0)])
    @fact res[1] => 50.4
    @fact -101.4 >= res[2] >= -102.0 => true
  end

  context("bounds if higher than max bound") do
    @fact BlackBoxOptim.rand_bound_from_target!([1.1], [1.0], [(0.0, 1.0)]) => [1.0]

    res = BlackBoxOptim.rand_bound_from_target!([97.0], [95.0], [(-10.0, 96.0)])
    @fact 95.0 <= res[1] <= 96.0 => true
  end
end

context("de_crossover_binomial") do
  context("always copies from donor if length is 1") do
    @fact BlackBoxOptim.de_crossover_binomial(DE, [0.0], 1, [1.0]) => [1.0]

    @fact BlackBoxOptim.de_crossover_binomial(DE, [-10.0], 1, [42.42]) => [42.42]
  end

  context("always copies at least one element from donor") do
    for(i in 1:NumTestRepetitions)
      len = rand(1:100)
      target, donor = rand(1, len), rand(1, len)
      res = BlackBoxOptim.de_crossover_binomial(DE, target, 1, donor)
      @fact any([ in(x, donor) for x = res ]) => true
    end
  end

  context("unlikely to copy everything if vectors are large") do
    for(i in 1:NumTestRepetitions)
      len = 50
      target, donor = rand(1, len), rand(1, len)
      res = BlackBoxOptim.de_crossover_binomial(DE, target, 1, donor)
      @fact any([ in(x, target) for x = res ]) => true
    end
  end
end

context("de_mutation_rand_1") do
  @fact ndims(BlackBoxOptim.de_mutation_rand_1(DE, 1, [4,9,8])) => 2

  @fact BlackBoxOptim.de_mutation_rand_1(DE, 1, [1, 2, 3])[1] => (3.0 + (0.4 * (1.0 - 2.0)))
  @fact BlackBoxOptim.de_mutation_rand_1(DE, 2, [2, 3, 1])[1] => (1.0 + (0.4 * (2.0 - 3.0)))
  @fact BlackBoxOptim.de_mutation_rand_1(DE, 3, [3, 2, 1])[1] => (1.0 + (0.4 * (3.0 - 2.0)))
  @fact BlackBoxOptim.de_mutation_rand_1(DE, 4, [1, 3, 5])[1] => (5.0 + (0.4 * (1.0 - 3.0)))
  @fact BlackBoxOptim.de_mutation_rand_1(DE, 5, [4, 9, 8])[1] => (8.0 + (0.4 * (4.0 - 9.0)))

  de2 = de_rand_1_bin({:SearchSpace => symmetric_search_space(2, (0.0, 10.0)),
    :Population => reshape([1.0:8.0], 4, 2),
    "f" => 0.6, "cr" => 0.5, "NumParents" => 3}
  )

  res = BlackBoxOptim.de_mutation_rand_1(de2, 10, [1,2,3])
  @fact ndims(res) => 2
  @fact res[1] => (3.0 + (0.6 * (1.0 - 2.0)))
  @fact res[2] => (7.0 + (0.6 * (5.0 - 6.0)))

  res2 = BlackBoxOptim.de_mutation_rand_1(de2, 2, [1,2,4])
  @fact res2[1] => (4.0 + (0.6 * (1.0 - 2.0)))
  @fact res2[2] => (8.0 + (0.6 * (5.0 - 6.0)))
end

context("ask") do
  for(i in 1:NumTestRepetitions)
    res = BlackBoxOptim.ask(DE)

    @fact length(res) => 2
    trial, trialIndex = res[1]
    target, targetIndex = res[2]

    @fact ndims(trial) => 2
    @fact 1 <= trialIndex <= length(DE.population) => true
    @fact isinspace(trial, DE.search_space) => true

    @fact ndims(target) => 2
    @fact 1 <= targetIndex <= length(DE.population) => true
    @fact isinspace(target, DE.search_space) => true

    @fact trialIndex == targetIndex => true
  end
end

end


text:
export AffineSubspace,
    AffineExtrinsic,
    AffineIntrinsic,
    Coordinates,
    Intrinsic,
    Extrinsic,
    intrinsic,
    extrinsic,
    dim,
    codim,
    ambient_dim,
    rand_affine_subspace,
    coord_change,
    geodesic_distance,
    geodesic

# References
const _LKK19 = """Lim, Lek-Heng, Ken Sze-Wai Wong, and Ke Ye. "Numerical algorithms on the affine Grassmannian." SIAM Journal on Matrix Analysis and Applications 40.2 (2019): 371-393"""

"""
    Coordinates

A type used for encoding the used coordinates and for performing coordinate changes.

Currently supported coordinates are:
- [`Intrinsic`](@ref)
- [`Extrinsic`](@ref)
"""
struct Coordinates{T} end
Base.broadcastable(C::Coordinates) = Ref(C)

"""
    Intrinsic <: Coordinates

Indicates the use of the intrinsic description of an affine subspace. See also [`AffineIntrinsic`](@ref).
"""
const Intrinsic = Coordinates{:Intrinsic}()

"""
    Extrinsic <: Coordinates

Indicates the use of the extrinsic description of an affine subspace. See also [`AffineExtrinsic`](@ref).
"""
const Extrinsic = Coordinates{:Extrinsic}()
# These are not part of the public API
const IntrinsicStiefel = Coordinates{:IntrinsicStiefel}()
const ExtrinsicStiefel = Coordinates{:ExtrinsicStiefel}()

"""
    AffineExtrinsic

Extrinsic description of an ``m``-dimensional affine subspace ``L`` in ``n``-dimensional space.
That is ``L = \\{ x | A x = b \\}``.
"""
struct AffineExtrinsic{T}
    A::Matrix{T}
    b::Vector{T}
end
(A::AffineExtrinsic)(x::AbstractVector) = A.A * x - A.b

"""
    dim(A::AffineExtrinsic)

Dimension of the affine subspace `A`.
"""
dim(A::AffineExtrinsic) = size(A.A, 2) - size(A.A, 1)

"""
    codim(A::AffineExtrinsic)

Codimension of the affine subspace `A`.
"""
codim(A::AffineExtrinsic) = size(A.A, 1)

Base.:(==)(A::AffineExtrinsic, B::AffineExtrinsic) = A.A == B.A && A.b == B.b
function Base.copy!(A::AffineExtrinsic, B::AffineExtrinsic)
    copy!(A.A, B.A)
    copy!(A.b, B.b)
    A
end
Base.copy(A::AffineExtrinsic) = AffineExtrinsic(copy(A.A), copy(A.b))

function Base.show(io::IO, A::AffineExtrinsic{T}) where {T}
    println(io, "AffineExtrinsic{$T}:")
    println(io, "A:")
    show(io, A.A)
    println(io, "\nb:")
    show(io, A.b)
end

Base.broadcastable(A::AffineExtrinsic) = Ref(A)

"""
    AffineIntrinsic

Intrinsic description of an ``m``-dimensional affine subspace ``L`` in ``n``-dimensional space.
That is ``L = \\{ u | A u + b₀ \\}``. Here, ``A`` and ``b₀`` are in orthogonal coordinates.
That is, the columns of ``A`` are orthonormal and ``A' b₀ = 0``.
"""
struct AffineIntrinsic{T}
    # orthogonal coordinates
    A::Matrix{T}
    b₀::Vector{T}
    # stiefel coordinates
    Y::Matrix{T}
end
AffineIntrinsic(A::Matrix, b₀::Vector) = AffineIntrinsic(A, b₀, stiefel_coordinates(A, b₀))

(A::AffineIntrinsic)(u::AbstractVector, ::Coordinates{:Intrinsic}) = A.A * u + A.b₀
(A::AffineIntrinsic)(u::AbstractVector, ::Coordinates{:IntrinsicStiefel}) = A.Y * u

function stiefel_coordinates(A::AbstractMatrix, b::AbstractVector)
    γ = sqrt(1 + sum(abs2, b))
    n, k = size(A)
    Y = zeros(eltype(A), n + 1, k + 1)
    Y[1:n, 1:k] .= A
    Y[1:n, k+1] .= b ./ γ
    Y[n+1, k+1] = 1 / γ
    Y
end

"""
    dim(A::AffineIntrinsic)

Dimension of the affine subspace `A`.
"""
dim(I::AffineIntrinsic) = size(I.A, 2)

"""
    codim(A::AffineIntrinsic)

Codimension of the affine subspace `A`.
"""
codim(I::AffineIntrinsic) = size(I.A, 1) - size(I.A, 2)

function Base.:(==)(A::AffineIntrinsic, B::AffineIntrinsic)
    A.A == B.A && A.b₀ == B.b₀ && A.Y == B.Y
end

function Base.show(io::IO, A::AffineIntrinsic{T}) where {T}
    println(io, "AffineIntrinsic{$T}:")
    println(io, "A:")
    show(io, A.A)
    println(io, "\nb₀:")
    show(io, A.b₀)
end

function Base.copy!(A::AffineIntrinsic, B::AffineIntrinsic)
    copy!(A.A, B.A)
    copy!(A.b₀, B.b₀)
    copy!(A.Y, B.Y)
    A
end
Base.copy(A::AffineIntrinsic) = AffineIntrinsic(copy(A.A), copy(A.b₀), copy(A.Y))

Base.broadcastable(A::AffineIntrinsic) = Ref(A)

function AffineIntrinsic(E::AffineExtrinsic)
    svd = LA.svd(E.A; full = true)
    b₀ = svd \ E.b
    m, n = size(E.A)
    A = Matrix((@view svd.Vt[m+1:end, :])')
    AffineIntrinsic(A, b₀)
end

function AffineExtrinsic(I::AffineIntrinsic)
    svd = LA.svd(I.A; full = true)
    m, n = size(I.A)
    A = Matrix((@view svd.U[:, n+1:end])')
    b = A * I.b₀
    AffineExtrinsic(A, b)
end

"""
    AffineSubspace(A, b)

An ``m``-dimensional affine subspace ``L`` in ``n``-dimensional space given by
the extrinsic description ``L = \\{ x | A x = b \\}``.

```julia-repl
julia> A = AffineSubspace([1 0 3; 2 1 3], [5, -2])
1-dim. affine subspace {x|Ax=b} with eltype Float64:
A:
2×3 Array{Float64,2}:
 1.0  0.0  3.0
 2.0  1.0  3.0
b:
2-element Array{Float64,1}:
  5.0
 -2.0

julia> dim(A)
1

julia> codim(A)
2

julia> ambient_dim(A)
3
```

An `AffineSubspace` holds always its [`extrinsic`](@ref) description, see also [`AffineIntrinsic`](@ref),
as well as its [`intrinsic`](@ref) description, see [`AffineExtrinsic`](@ref).

```julia-repl
julia> intrinsic(A)
AffineIntrinsic{Float64}:
A:
3×1 Array{Float64,2}:
 -0.6882472016116853
  0.6882472016116853
  0.22941573387056186
b₀:
3-element Array{Float64,1}:
 -3.0526315789473677
 -3.947368421052632
  2.684210526315789
```

An `AffineSubspace` can be evaluated with either using [`Intrinsic`](@ref)
or [`Extrinsic`](@ref) coordinates.

```julia-repl
julia> u = [0.5]
1-element Array{Float64,1}:
 0.5

julia> x = A(u, Intrinsic)
3-element Array{Float64,1}:
 -3.3967551797532103
 -3.6032448202467893
  2.79891839325107

julia> A(x, Extrinsic)
  2-element Array{Float64,1}:
   0.0
   0.0
```

To change the used coordinates you can use [`coord_change`](@ref).
```julia-repl
julia> coord_change(A, Extrinsic, Intrinsic, x)
1-element Array{Float64,1}:
 0.49999999999999994
```
"""
struct AffineSubspace{T}
    extrinsic::AffineExtrinsic{T}
    intrinsic::AffineIntrinsic{T}
end

AffineSubspace(I::AffineIntrinsic) = AffineSubspace(AffineExtrinsic(I), I)
AffineSubspace(E::AffineExtrinsic) = AffineSubspace(E, AffineIntrinsic(E))

function AffineSubspace(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
    size(A, 1) == length(b) || throw(ArgumentError("Size of A and b not compatible."))
    0 < size(A, 1) < size(A, 2) ||
        throw(ArgumentError("Affine subspace has to be given in extrinsic coordinates, i.e., by A x = b."))

    AffineSubspace(AffineExtrinsic(Matrix(float.(A)), Vector(float.(b))))
end

Base.broadcastable(A::AffineSubspace) = Ref(A)

"""
    dim(A::AffineSubspace)

Dimension of the affine subspace `A`.
"""
dim(A::AffineSubspace) = dim(A.intrinsic)

"""
    codim(A::AffineSubspace)

Codimension of the affine subspace `A`.
"""
codim(A::AffineSubspace) = codim(A.intrinsic)

"""
    ambient_dim(A::AffineSubspace)

Dimension of ambient space of the affine subspace `A`.
"""
ambient_dim(A::AffineSubspace) = dim(A) + codim(A)

function Base.show(io::IO, A::AffineSubspace{T}) where {T}
    println(io, "$(dim(A))-dim. affine subspace {x|Ax=b} with eltype $T:")
    println(io, "A:")
    show(io, A.extrinsic.A)
    println(io, "\nb:")
    show(io, A.extrinsic.b)
end

"""
    intrinsic(A::AffineSubspace)

Obtain the intrinsic description of `A`, see also [`AffineIntrinsic`](@ref).
"""
intrinsic(A::AffineSubspace) = A.intrinsic

"""
    extrinsic(A::AffineSubspace)

Obtain the extrinsic description of `A`, see also [`AffineExtrinsic`](@ref).
"""
extrinsic(A::AffineSubspace) = A.extrinsic

function Base.copy!(A::AffineSubspace, B::AffineSubspace)
    copy!(A.intrinsic, B.intrinsic)
    copy!(A.extrinsic, B.extrinsic)
    A
end
Base.copy(A::AffineSubspace) = AffineSubspace(copy(A.extrinsic), copy(A.intrinsic))

function Base.:(==)(A::AffineSubspace, B::AffineSubspace)
    intrinsic(A) == intrinsic(B) && extrinsic(A) == extrinsic(B)
end

function (A::AffineSubspace)(x::AbstractVector, ::Coordinates{:Intrinsic})
    intrinsic(A)(x, Intrinsic)
end
function (A::AffineSubspace)(x::AbstractVector, ::Coordinates{:IntrinsicStiefel})
    intrinsic(A)(x, IntrinsicStiefel)
end
function (A::AffineSubspace)(x::AbstractVector, ::Coordinates{:Extrinsic})
    extrinsic(A)(x)
end

"""
    rand_affine_subspace(n::Integer; dim | codim, real = false)

Generate a random [`AffineSubspace`](@ref) with given dimension `dim` or codimension `codim`
(one of them has to be provided) in ambient space of dimension `n`.
If `real` is `true`, then the extrinsic description is real.
The subspace is generated by drawing each entry of the extrinsic description indepdently
from a normal distribuation using `randn`.

## Example
```julia-repl
julia> rand_affine_subspace(3; dim = 1)
1-dim. affine subspace {x|Ax=b} with eltype Complex{Float64}:
A:
2×3 Array{Complex{Float64},2}:
  -1.73825+1.27987im   -0.0871343+0.840408im  -0.551957+0.106397im
 -0.597132-0.343965im   -0.122543-0.172715im   -1.04949+0.370917im
b:
2-element Array{Complex{Float64},1}:
  0.47083334430689394 + 0.8099804422599071im
 -0.12018696822943896 + 0.11723026326952792im

julia> rand_affine_subspace(4; codim = 1)
3-dim. affine subspace {x|Ax=b} with eltype Complex{Float64}:
A:
1×4 Array{Complex{Float64},2}:
 0.345705+0.0893881im  -0.430867-0.663249im  0.979969-0.569378im  -0.29722-0.192493im
b:
1-element Array{Complex{Float64},1}:
 0.7749708228192062 + 0.9762873764567546im
```
"""
function rand_affine_subspace(
    n::Integer;
    dim::Union{Nothing,Integer} = nothing,
    codim::Union{Nothing,Integer} = nothing,
    real::Bool = false,
)
    !isnothing(dim) ||
        !isnothing(codim) ||
        throw(ArgumentError("Neither `dim` or `codim` specified."))

    if !isnothing(dim)
        0 < dim < n || throw(ArgumentError("`dim` has to be between 0 and `n`."))
        k = dim
    else
        0 < codim < n || throw(ArgumentError("`codim` has to be between 0 and `n`."))
        k = n - codim
    end
    T = real ? Float64 : ComplexF64
    A = randn(T, n - k, n)
    b = randn(T, n - k)

    AffineSubspace(A, b)
end

# Coordinate changes
"""
    coord_change(A::AffineSubspace, C₁::Coordinates, C₂::Coordinates, p)

Given an affine subspace `A` and a point `p` in coordinates `C₁` compute the point
`x` describing p in coordinates `C₂`.

## Example
```julia-repl
julia> A = AffineSubspace([1 0 3; 2 1 3], [5, -2]);

julia> u = [1.25];

julia> x = coord_change(A, Intrinsic, Extrinsic, u)
3-element Array{Float64,1}:
 -3.9129405809619744
 -3.087059419038025
  2.9709801936539915

julia> A(x, Extrinsic)
2-element Array{Float64,1}:
 0.0
 0.0

julia> x - A(u, Intrinsic)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0
```
"""
coord_change(A::AffineSubspace, ::C, ::C, x) where {C<:Coordinates} = x
coord_change(A::AffineSubspace, ::Coordinates{:Intrinsic}, ::Coordinates{:Extrinsic}, u) =
    A(u, Intrinsic)
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:Intrinsic},
    ::Coordinates{:IntrinsicStiefel},
    u,
)
    [
        u
        1 / A.intrinsic.Y[end, end]
    ]
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:Intrinsic},
    ::Coordinates{:ExtrinsicStiefel},
    u,
)
    v = coord_change(A, Intrinsic(), IntrinsicStiefel(), u)
    coord_change(A, IntrinsicStiefel(), ExtrinsicStiefel(), v)
end

coord_change(A::AffineSubspace, ::Coordinates{:Extrinsic}, ::Coordinates{:Intrinsic}, x) =
    A.intrinsic.A' * (x - A.intrinsic.b₀)
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:Extrinsic},
    ::Coordinates{:IntrinsicStiefel},
    x,
)
    u = coord_change(A, Extrinsic(), Intrinsic(), x)
    coord_change(A, Intrinsic(), IntrinsicStiefel(), u)
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:Extrinsic},
    ::Coordinates{:ExtrinsicStiefel},
    x,
)
    [x; 1]
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:IntrinsicStiefel},
    ::Coordinates{:Intrinsic},
    u,
)
    γ = A.intrinsic.Y[end, end]
    [u[i] / (γ * u[end]) for i = 1:(length(u)-1)]
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:IntrinsicStiefel},
    ::Coordinates{:Extrinsic},
    u,
)
    x = coord_change(A, IntrinsicStiefel(), ExtrinsicStiefel(), u)
    coord_change(A, ExtrinsicStiefel(), Extrinsic(), x)
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:IntrinsicStiefel},
    ::Coordinates{:ExtrinsicStiefel},
    u,
)
    A(u, IntrinsicStiefel)
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:ExtrinsicStiefel},
    ::Coordinates{:Intrinsic},
    x,
)
    u = coord_change(A, ExtrinsicStiefel(), IntrinsicStiefel(), x)
    coord_change(A, IntrinsicStiefel(), Intrinsic(), u)
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:ExtrinsicStiefel},
    ::Coordinates{:Extrinsic},
    x,
)
    x[1:end-1] ./ x[end]
end
function coord_change(
    A::AffineSubspace,
    ::Coordinates{:ExtrinsicStiefel},
    ::Coordinates{:IntrinsicStiefel},
    x,
)
    A.intrinsic.Y' * x
end

"""
    geodesic_distance(A::AffineSubspace, B::AffineSubspace)

Compute the geodesic distance between `A` and `B` in the affine Grassmanian `Graff(k, n)`
where `k = dim(A)` and `n` is the amebient dimension.
This follows the derivation in [^LKK19].

[^LKK19]: $_LKK19.
"""
geodesic_distance(A::AffineSubspace, B::AffineSubspace) =
    geodesic_distance(A.intrinsic, B.intrinsic)
function geodesic_distance(A::AffineIntrinsic, B::AffineIntrinsic)
    sqrt(sum(σᵢ -> acos(min(σᵢ, 1.0))^2, LA.svdvals(A.Y' * B.Y)))
end

# geodesic
"""
    geodesic_svd(A::AffineSubspace, B::AffineSubspace)

Computes the factors ``Q``, ``Θ`` and ``U`` from Corollary 4.3 in [^LKK19].
These values are necessary to construct the geodesic between `A` and `B`.

[^LKK19]: $_LKK19
"""
geodesic_svd(A::AffineSubspace, B::AffineSubspace) = geodesic_svd(A.intrinsic, B.intrinsic)
function geodesic_svd(A::AffineIntrinsic, B::AffineIntrinsic)
    U, Σ, V = LA.svd(A.Y' * B.Y)
    # Θ = acos.(min.(Σ, 1.0))
    # We have acos(1-eps()) = 2.1073424255447017e-8
    # So this is numerically super unstable if the singular value is off by only eps()
    # We try to correct this fact by treating σ >= 1 - 2eps() as 1.0
    Θ = map(σ -> σ + 2 * eps() > 1.0 ? 0.0 : acos(σ), Σ)

    n, k = size(A.Y)
    # M = (LA.I - A.Y * A.Y') * B.Y * inv(A.Y' * B.Y)
    # Have to compute an SVD of M s.t. M = Q tanΘ U'
    # Equivalently M * U = Q tan(Θ)
    # We can achieve this by using a *pivoted* QR
    # since then R will be a diagonal matrix s.t. the absolute value of R_ii is θ_{k-i}

    # inv(A.Y' * B.Y) * U = V * LA.diagm(inv.(Σ))
    MU = (LA.I - A.Y * A.Y') * B.Y * V * LA.diagm(inv.(Σ))
    Q, R = LA.qr(MU, Val(true))
    # correct signs and ordering of Q
    Q′ = Q[:, k:-1:1]
    for j = 1:k
        # Look if we need to flip signs
        real(R[k-j+1, k-j+1]) < 0 || continue
        for i = 1:n
            Q′[i, j] = -Q′[i, j]
        end
    end

    Q′, Θ, U
end

"""
    geodesic(A::AffineSubspace, B::AffineSubspace)

Returns the geodesic ``γ(t)`` connecting `A` and `B` in the Grassmanian ``Gr(k+1,n+1)``
where ``k`` is the dimension of ``A`` and ``n`` is the ambient dimension.
See also Corollary 4.3 in [^LKK19].

[^LKK19]: $_LKK19
"""
geodesic(A::AffineSubspace, B::AffineSubspace) = geodesic(A.intrinsic, B.intrinsic)
function geodesic(A::AffineIntrinsic, B::AffineIntrinsic)
    Q, Θ, U = geodesic_svd(A, B)
    t -> A.Y * U * LA.diagm(cos.(t .* Θ,)) * U' + Q * LA.diagm(sin.(t .* Θ,)) * U'
end
#


text:
#=
ZygoteRules.@adjoint function ODESolution(u,args...)
  function ODESolutionAdjoint(ȳ)
    (ȳ,ntuple(_->nothing, length(args))...)
  end
  ODESolution(u,args...), ODESolutionAdjoint
end
=#

ZygoteRules.@adjoint function ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11
                }(u,args...) where {T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}

                function ODESolutionAdjoint(ȳ)
                  (ȳ,ntuple(_->nothing, length(args))...)
                end

                ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}(u,args...),ODESolutionAdjoint
end

ZygoteRules.@adjoint function ZygoteRules.literal_getproperty(sol::AbstractTimeseriesSolution, ::Val{:u})
  function solu_adjoint(Δ)
        zerou = zero(sol.prob.u0)
        _Δ = @. ifelse(Δ == nothing,(zerou,),Δ)
        (DiffEqBase.build_solution(sol.prob,sol.alg,sol.t,_Δ),)
  end
  sol.u,solu_adjoint
end

ZygoteRules.@adjoint function ZygoteRules.literal_getproperty(sol::AbstractNoTimeSolution, ::Val{:u})
  function solu_adjoint(Δ)
        zerou = zero(sol.prob.u0)
        _Δ = @. ifelse(Δ == nothing,(zerou,),Δ)
        (DiffEqBase.build_solution(sol.prob,sol.alg,_Δ,sol.resid),)
  end
  sol.u,solu_adjoint
end

ZygoteRules.@adjoint function (f::ODEFunction)(u,p,t)
  if f.vjp === nothing
    ZygoteRules._pullback(f.f,u,p,t)
  else
    f.vjp(u,p,t)
  end
end

ZygoteRules.@adjoint! function (f::ODEFunction)(du,u,p,t)
  if f.vjp === nothing
    ZygoteRules._pullback(f.f,du,u,p,t)
  else
    f.vjp(du,u,p,t)
  end
end

ZygoteRules.@adjoint function SciMLBase.tmap(f, args::Union{AbstractArray,Tuple}...)
  ∇tmap(__context__, f, args...)
end

ZygoteRules.@adjoint function SciMLBase.responsible_map(f, args::Union{AbstractArray,Tuple}...)
  ∇responsible_map(__context__, f, args...)
end

ZygoteRules.@adjoint function DiffEqBase.EnsembleSolution(sim,time,converged)
  out = EnsembleSolution(sim,time,converged)
  function EnsembleSolution_adjoint(p̄::AbstractArray{T,N}) where {T,N}
    arrarr = [[p̄[ntuple(x->Colon(),Val(N-2))...,j,i] for j in 1:size(p̄)[end-1]] for i in 1:size(p̄)[end]]
    (EnsembleSolution(arrarr, 0.0, true),nothing,nothing)
  end
  function EnsembleSolution_adjoint(p̄::EnsembleSolution)
    (p̄,nothing,nothing)
  end
  out,EnsembleSolution_adjoint
end
ZygoteRules.@adjoint ZygoteRules.literal_getproperty(sim::EnsembleSolution, ::Val{:u}) = sim.u, p̄ -> (EnsembleSolution(p̄, 0.0, true),)

#=
ChainRulesCore.frule(f::ODEFunction,u,p,t)
  if f.jvp === nothing
    ChainRulesCore.frule(f.f,u,p,t)
  else
    function ode_jvp(f,du,dp,dt)
      f.jvp_u(du,u,p,t) + f.jvp_p(dp,u,p,t) + f.tgrad(u,p,t)*dt
    end
    f.f(u,p,t),ode_jvp
end
=#

function ChainRulesCore.rrule(f::ODEFunction,u,p,t)
  if f.vjp === nothing
    ChainRulesCore.rrule(f.f,u,p,t)
  else
    f.vjp(u,p,t)
  end
end

ZygoteRules.@adjoint numargs(f) = (numargs(f),df->(nothing,))
ChainRulesCore.rrule(::typeof(numargs),f) = (numargs(f),df->(nothing,))

# Until https://github.com/FluxML/Zygote.jl/issues/664 is fixed
ZygoteRules.@adjoint function Base.pairs(x::T) where T
         y = Base.pairs(x)
         back(dx::NamedTuple) = (dx.data,)
         function back(dx::Dict)
           T <: AbstractDict && return (dx,)
           z = zero(x)
           for (k,v) in dx
             z[k] = v
           end
           (z,)
         end
         y, back
       end


text:
module OEIS

using Requests, Iterators

import Base: show
export oeis

# Parse OEIS internal format
# see https://oeis.org/eishelp1.html for a description
# we are interested in:
#   %I id number (A followed by 6-digit number)
#   %S, %T, %U absolute values of terms (%T, %U optional)
#   %V, %W, %X signed values (only if sequence contains negative numbers)
#   %N name

type IntegerSequence{T<:Integer}
    id::Symbol
    name::ASCIIString
    values::Vector{T}
end

function show(io::IO,seq::IntegerSequence)
    println(io, string(typeof(seq)), ": ", string(seq.id))
    println(io, "  ", seq.name)
    if length(seq.values) <= 20
        print(io, "  ", join(seq.values,','))
    else
        print(io, "  ", join(seq.values[1:20],','),",...")
    end
end

# there has to be a better way than this...
bestint(s::String) = macroexpand(parse(s))

function parse_oeis_page(s::String)
    if contains(s,"Too many results. Please narrow search.")
        error("Too many results. Please narrow search.")
    end
    seqs = IntegerSequence[]
    lines = filter(l -> length(l)>=10 && l[1]=='%', split(s,'\n') )
    for r in groupby(l -> l[4:10], lines)
        push!(seqs, parse_oeis_record(r))
    end
    seqs
end

function parse_oeis_record(r)
    local id, name, vec
    for line in r
        C = line[2]
        if C == 'I'
            id = line[4:10]
        elseif C == 'N'
            name = line[12:end]
        elseif C in ('S','V')
            vec = Integer[bestint(s) for s in split(line[12:end],',',0,false)]
        elseif C in ('T','U','W','X')
            append!(vec,Integer[bestint(s) for s in split(line[12:end],',',0,false)])
        end
    end
    if all(i -> isa(i,Int), vec)
        vec = Int[i for i in vec]
    end
    IntegerSequence(symbol(id),bytestring(name),vec)
end

function oeis(query::String)
    r = get("http://oeis.org/search"; query={"fmt" => "text", "q" => query})
    parse_oeis_page(r.data)
end

function oeis(id::Symbol)
    r = get("http://oeis.org/search"; query={"fmt" => "text", "q" => string("id:",id)})
    parse_oeis_page(r.data)[1]
end

function oeis(seq::Array)
    r = get("http://oeis.org/search"; query={"fmt" => "text", "q" => join(seq,',')})
    parse_oeis_page(r.data)
end

end # module


text:
# ---
# title: 122. Best Time to Buy and Sell Stock II
# id: problem122
# author: AquaIndigo
# date: 2020-11-08
# difficulty: Easy
# categories: Array, Greedy
# link: <https://leetcode.com/problems/best-time-to-buy-and-sell-stock-ii/description/>
# hidden: true
# ---
#
# Say you have an array `prices` for which the _i_ th element is the price of a
# given stock on day _i_.
#
# Design an algorithm to find the maximum profit. You may complete as many
# transactions as you like (i.e., buy one and sell one share of the stock
# multiple times).
#
# **Note:** You may not engage in multiple transactions at the same time (i.e.,
# you must sell the stock before you buy again).
#
# **Example 1:**
#
#
#
#     Input: [7,1,5,3,6,4]
#     Output: 7
#     Explanation: Buy on day 2 (price = 1) and sell on day 3 (price = 5), profit = 5-1 = 4.
#                  Then buy on day 4 (price = 3) and sell on day 5 (price = 6), profit = 6-3 = 3.
#
#
# **Example 2:**
#
#
#
#     Input: [1,2,3,4,5]
#     Output: 4
#     Explanation: Buy on day 1 (price = 1) and sell on day 5 (price = 5), profit = 5-1 = 4.
#                  Note that you cannot buy on day 1, buy on day 2 and sell them later, as you are
#                  engaging multiple transactions at the same time. You must sell before buying again.
#
#
# **Example 3:**
#
#
#
#     Input: [7,6,4,3,1]
#     Output: 0
#     Explanation: In this case, no transaction is done, i.e. max profit = 0.
#
#
#
# **Constraints:**
#
#   * `1 <= prices.length <= 3 * 10 ^ 4`
#   * `0 <= prices[i] <= 10 ^ 4`
#
#
## @lc code=start
using LeetCode

function max_profit(prices::Vector{Int})::Int
    res = 0
    for i in 2:length(prices)
        res += (prices[i] > prices[i - 1]) ? (prices[i] - prices[i - 1]) : 0
    end
    return res
end
## @lc code=end


text:
using LinearAlgebra
using DynamicPolynomials
using SwitchOnSafety
using Combinatorics
using SparseArrays
using JuMP, MosekTools
using SpecialFunctions

include("../src/RandomTrajectories.jl")
include("../src/AlgebraicLift.jl")
include("../src/ScenarioOpti.jl")
include("../src/ProbabilisticCertificates.jl")
include("../src/WhiteBoxAnalysis.jl")

dim = 6; numMode = 6; dimIn = Integer(floor(dim/2)+1)

numScen_budget = 40000

#beta = 0.01
#batchsize = 1000
tol=1e-3

A = [2*rand(Float64, (dim, dim)).-1 for i=1:numMode]#generate_switched_linear_systems(numMode,dim)
B = 2*rand(Float64, (dim, dimIn)).-1

#jsrboundopen= white_box_jsr(A)
#println("JSR open loop: $jsrboundopen")
gaWb = white_box_stabilization_quad(A,B)
println("White-box stabilization: $gaWb")

(state0_budget,state_budget) = generate_trajectories(1,A,numScen_budget)

jsrMat = []

K0 = zeros(dimIn,dim)

for N in 20000:1000:numScen_budget
    state0 = state0_budget[:,1:N]
    state = state_budget[:,1:N]
    K,jsr_bound = probabilistc_stability_certificate(state0,state;B=B,numMode=numMode,d=1,batchsize=1000,K0=K0,beta=0.01,tol=1e-3)
    push!(jsrMat,jsr_bound)
    println(jsrMat)
    println("White-box stabilization: $gaWb")
end

#n,M,m = 3,3,2
#Any[1.6905964015508155, 1.4949948318912094, 1.4578897728330178, 1.4449970570350583, 1.3661954310327227, 1.3557547083722759, 1.3519850357471235, 1.366442095956408, 1.3613413759324937, 1.3429185951062852, 1.3492377476476396, 1.353041725614804, 1.350881519910859, 1.351024882306086, 1.332142623005733, 1.3332099256037162, 1.3304895958561995, 1.3293168479495931, 1.3257398747945923, 1.3242527214300455]
#White-box stabilization: (1.3046225525634858, [-0.08704353863380242 -0.014885289542307709 -0.049119276805454554; -0.507914925287981 0.346359451257569 1.0595139493317263], [0.4730751424314884 0.05238239264562874 -0.11595043774240237; 0.05238239264562874 0.47706738397162957 -0.06627547776094828; -0.11595043774240237 -0.06627547776094828 0.5394594622004449])
#1000:1000:40000
#Any[1.5729841161676168, 1.0811554227453197, 1.1204878768144788, 0.9379157526294873, 1.0825047220475408, 0.9230141022655761, 0.9067207134210432, 0.9086742873941583, 0.9032261428457313, 0.9074603880590151, 0.9025114517834851, 0.9188062149754335, 0.879041002412628, 0.8870911981713733, 0.8926992150562685, 0.9244034116395815, 0.9213403526016543, 0.9186221154258161, 0.8857493686041322, 0.8845971086454935, 0.8796659173161556, 0.8787210471488426, 0.8778585397515366, 0.8770655056628548, 0.8763348706011853, 0.8756594770044023, 0.8816237449635299, 0.8810259381333286, 0.8803501376108269, 0.8799417718551741, 0.8794587848269098, 0.8790053814525175, 0.8785788973412633, 0.8752015696756477, 0.8748238216588682, 0.860271590366133, 0.8598788107531674, 0.8595062987770098, 0.8591372241094922, 0.8767243697570494]
#White-box stabilization: (0.8140390017107871, [0.4613815512089252 -0.5074407448126327 0.6781102279468121; -0.4054564358304964 0.9172867615097666 -0.5910686302058721], [0.1575553288978803 0.02705742403492229 -0.09788071638611223; 0.02705742403492229 0.3447930410671564 -0.04099837698782956; -0.09788071638611222 -0.04099837698782956 0.636873987762721])

#n,M,m = 4,4,3
#Any[Inf, Inf, Inf, Inf, 8.180216224894956, 6.799812206868395, 3.9527055267977347, 4.59174170956012, 3.0846800593837282, 2.6548561907838444, 2.511122991060089, 2.29880083863334, 2.161273525017097, 2.155982590641091, 2.0759785750943713, 2.6898719706804557, 2.3931499521503685, 2.048817168640575, 1.9826781075826072, 2.167094635668247]
#White-box stabilization: (1.2912508949641772, [-0.29467773970869054 0.45120789516804116 -0.044747113301878305 -1.3412127586682072; 0.29469719061888044 -0.12455867244325437 0.47017019594572057 -0.812536496343961; 0.46212127494925737 0.30630374660913356 0.03271426927287058 -1.5218113173193089], [0.3425385468310902 0.008731320963757413 0.13821082035445126 0.18629022929560726; 0.008731320963757404 0.27313481061590716 0.081610108315615 0.027369650943003167; 0.13821082035445126 0.08161010831561501 0.2951065970488265 0.08589689280299384; 0.1862902292956073 0.02736965094300317 0.08589689280299384 0.4106130096540772])
#1000:1000:40000
#Any[Inf, Inf, 7.875913638698062, 7.17815675416856, 3.4416932736968016, 19.394289213014666, 2.583723965590616, 4.935647364390182, 1.984015772805858, 1.8674821390119283, 1.851433761685506, 2.2556567432257753, 2.1083617293205275, 2.0562587423170604, 1.720919007398438, 1.7461732618949821, 1.782734889401107, 1.7091518413970257, 1.7326683668613923, 1.6920558239323553, 1.7738463864124947, 1.8309626179003315, 1.7467378224996795, 1.6560635313126297, 1.650624547634913, 1.82414577393583, 1.6025444213063642, 1.6289801783733817, 1.745663331430371, 1.609830928758939, 1.597559164787974, 1.5904782457385636, 1.5838032484067968, 1.577498470309936, 1.5715323643953338, 1.6089367460744595, 1.553410839521955, 1.657015603207103, 1.5468122074359238, 1.546056871905595]
#White-box stabilization: (1.2337473357054138, [-0.007571043628938606 0.09068410059566147 -0.008151181227635546 -0.23763939318748917; 0.5431810653267845 0.25227603574433566 -0.1337727867212335 0.044797330378517915; 0.24885446041091436 0.4054576299301648 -0.3586486818168962 0.1656272016965966], [0.24836613958342035 0.09944045818289737 -0.16527107495145638 -0.12057885330835569; 0.09944045818289737 0.3576493794229119 -0.04223338581847068 -0.05938473714569832; -0.16527107495145635 -0.04223338581847068 0.428282720662584 0.07401563630157972; -0.12057885330835569 -0.05938473714569833 0.07401563630157973 0.5431548947140289])

#n,M,m = 4,6,3
#Any[Inf, Inf, Inf, 6.737095856110309, 3.675632937907526, 4.711663285904102, 3.6036786944899486, 3.245532901636809, 3.179705851172089, 2.928079175769009, 2.786632175891354, 2.5987061969012726, 2.690047270463758, 2.56824929071184, 2.3637115663198616, 3.0243999806253585, 3.0800995588355096, 2.7689719140449247, 2.640541451762192, 2.5793649540455488]
#White-box stabilization: (1.411688319459564, [0.18773192345444614 0.330514311888444 0.08303501025300539 0.6688093068290499; 0.096991858004839 -0.1376269451782156 0.16258760593673052 -0.28712484054435716; -0.05179889939849871 -0.21974552189934024 0.019847110276567503 -0.574729653178076], [0.3048746641659247 0.08704701508549992 -0.022825373701927137 0.03258573105931016; 0.08704701508549993 0.3155721631766691 -0.033540595048654746 0.01081474042260808; -0.022825373701927144 -0.033540595048654746 0.4273011813464707 -0.19692622371349458; 0.03258573105931017 0.01081474042260808 -0.19692622371349458 0.4285114388792512])
#1000:1000:40000
#Any[Inf, Inf, Inf, 25.084757321839195, 8.770136588631802, 9.976749384135678, 6.877930810422303, 6.939754176855923, 5.970890225057548, 5.158585131735316, 7.026520008542251, 3.9068211007996236, 2.7116440033548512, 2.6401441896367848, 2.6546224672270173, 2.479164798426021, 2.4708216864342907, 2.460497586092857, 2.429913835468321, 2.3758813963902914, 2.3461590747282215, 2.5836376981383578, 2.5429050265257245, 2.509091638025923, 2.4782728932067286, 2.4500540776504156, 2.4383169713775903, 2.4144773087877156, 2.3329332189333867, 2.251096846759749, 2.238871368744152, 2.2247113742266023, 2.17326016039637, 2.210486799215429, 2.198078155410338, 2.186357054727009, 2.2937024450686985, 2.268810905462152, 2.2561233423195572, 2.244082132489508]
#White-box stabilization: (1.8347012279895227, [1.2413485353963243 -0.13936635314223478 0.24391546443723508 -0.6272971813182449; 0.11882472083532031 0.006835856413333752 0.033442637097538024 -0.017462245766754272; 0.4572432136771731 -0.14943479419562547 0.11084913373271818 0.1330692144943799], [0.5288623013381373 -0.06955297816624266 0.07739492604379934 -0.14945226904731068; -0.06955297816624267 0.5133593931137286 0.024751588463856253 -0.06508656604630439; 0.07739492604379934 0.024751588463856253 0.4582612910506892 0.02558559352254509; -0.14945226904731068 -0.06508656604630439 0.02558559352254509 0.5703552027639498])

#n,M,m = 6,6,4
#30000:1000:40000
#ny[8.072715651100154, 7.768234690251989, 8.180570339667353, 7.886210946675976, 8.14582975086884, 7.871660888641335, 7.1033528127184, 7.390710250162713, 7.431287049154341, 5.567010043658964, 5.467431313638779]
#White-box stabilization: (1.8494425061800728, [0.012221018509647168 0.5312774075791682 0.038976706966226404 0.24115751812320538 0.29963266504137354 -0.24193177579705996; -0.28880253871014694 -0.23324503363629637 -0.030309376528313354 0.007505709936567113 -0.12073558590288352 0.15846373101244968; -0.20699525532892935 -0.4544479559088058 0.11982743997670932 -0.039258464791142526 -0.3181309400642283 0.3150143828761531; -0.0139831817182914 0.5874777519813043 -0.10178037874387551 0.12249922444250158 0.49381698904747273 -0.023525512008606104], [0.2629781490641416 0.0003384808471250383 -0.014510051589145008 -0.03191164850633556 0.013452938609376693 0.006781828774616212; 0.00033848084712503846 0.29904793232982807 0.0283758893731878 -0.00703246955168064 -0.05022634263568734 0.05898316898301458; -0.014510051589145012 0.0283758893731878 0.32192795235731525 -0.02850577360183669 -0.02556378959490166 0.02713605333508494; -0.03191164850633556 -0.007032469551680639 -0.028505773601836697 0.3370007449715669 -0.02596961925645344 0.06670614725088554; 0.013452938609376693 -0.05022634263568735 -0.02556378959490166 -0.02596961925645344 0.3207606635040778 -0.07275525717780275; 0.006781828774616212 0.05898316898301458 0.027136053335084936 0.06670614725088554 -0.07275525717780275 0.2882382320544729])
#20000:1000:40000
#Any[107.0843979270694, 59.342126663772156, 22.53861363197418, 19.35338764565481, 25.84527685861678, 21.968851042329128, 19.689300883307723, 49.614306286747734, 37.50510470585297, 10.732596162611243, 12.937874087207211, 8.87777827552819, 7.665576388216101, 7.2410821009579, 7.041489494504242, 6.996118755748454, 6.822953474887835, 6.505336929771839, 7.494983714450516, 6.310508993177824, 6.189122196716375]
#White-box stabilization: (1.9120856741024894, [0.2397554519049773 -0.10131199227907627 0.23889381969404355 -0.46808685505580094 -0.2830979655692059 -0.20420265706522728; 0.34396199478718237 0.04113041958533565 0.46428253981645234 -0.0810038935859582 -0.7453144451913242 -0.29377911940563817; 0.3294039638069217 -0.35268064974948693 0.423490743615612 -0.27429288114813877 -0.4997957061387201 -0.2536533942553235; 0.3459759820407083 -0.10603493805154796 0.38431492409875784 -0.16744635480413356 -1.045721513883651 0.0009893083224673024], [0.46126631090810905 0.0138872832302601 0.21125527828557872 0.12615429354805274 0.03787693562846611 -0.1251841989029774; 0.013887283230260098 0.41321690740308814 0.04124657326637017 -0.04856644977200873 0.05513722724728301 -0.11742322672944541; 0.21125527828557875 0.041246573266370164 0.6074761020842684 0.12138653021576508 -0.027014667936934723 -0.024310740534286227; 0.12615429354805274 -0.04856644977200873 0.12138653021576508 0.4420956139160551 0.06134324952849663 -0.07777124466288073; 0.03787693562846611 0.05513722724728302 -0.027014667936934723 0.061343249528496624 0.37137744074642076 -0.09637505917692687; -0.1251841989029774 -0.11742322672944543 -0.024310740534286227 -0.07777124466288071 -0.09637505917692687 0.44820687724218183])

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
function _obj_function(p,parameter_dictionary::Dict{String,Any}, v_variable_data::Array{Float64,1},
    experimental_data_array::Array{Float64,2})

    parameter_dictionary["KD"] = p[1]
    parameter_dictionary["n"] = p[2]

    # compute the model v-variable -
    model_v_variable = compute_activity_variable(experimental_data_array[:,1], parameter_dictionary)

    # compute the residual between the data and model estimated v-variable -
    error = sum((v_variable_data .- model_v_variable).^2)

    # return -
    return error
end

function estimate_parameters(parameter_dictionary::Dict{String,Any}, v_variable_data::Array{Float64,1},
    experimental_data_array::Array{Float64,2})

    try

        # get unkown parameter guess from config -
        KD = parameter_dictionary["KD"]
        n = parameter_dictionary["n"]

        # setup the calculation -
        xinitial = [KD n]
        OF(p) = _obj_function(p,parameter_dictionary, v_variable_data, experimental_data_array)

        # call the optimizer -
        opt_result = optimize(OF,xinitial,BFGS())

        # return -
        return Optim.minimizer(opt_result)
    catch error
        throw(error)
    end
end


text:
using XGBoost

include("utils.jl")

# we load in the agaricus dataset
# In this example, we are aiming to predict whether a mushroom can be eated

# we use auxiliary function to read LIBSVM format into julia Matrix
train_X, train_Y = readlibsvm("../data/agaricus.txt.train", (6513, 126))
test_X, test_Y = readlibsvm("../data/agaricus.txt.test", (1611, 126))

#-------------Basic Training using XGBoost-----------------
# note: xgboost naturally handles sparse input
# use sparse matrix when your feature is sparse(e.g. when you using one-hot encoding vector)
# model paramters can be set as paramter for ```xgboost``` function, or use an Array{(ASCIIString, Any), 1}/Dict()
num_round = 2

print("training xgboost with dense matrix\n")
# you can directly pass julia's matrix or sparse matrix as data,
#   by calling xgboost(data, num_round, label=label, training-parameters)
bst = xgboost(train_X, num_round, label = train_Y, eta=1, max_depth=2, objective="binary:logistic")

print("training xgboost with sparse matrix\n")
sptrain = sparse(train_X)
# alternatively, you can pass parameters in as a map
param = Dict("max_depth"=>2, "eta"=>1, "objective"=>"binary:logistic")
bst = xgboost(sptrain, num_round, label = train_Y, param=param)

# you can also put in xgboost's DMatrix object
# DMatrix stores label, data and other meta datas needed for advanced features
print("training xgboost with DMatrix")
dtrain = DMatrix(train_X, label = train_Y)
bst = xgboost(dtrain, num_round, eta = 1, objective = "binary:logistic")

# you can also specify data as file path to a LibSVM format input
bst = xgboost("../data/agaricus.txt.train", num_round, max_depth = 2, eta = 1, objective = "binary:logistic")

#--------------------basic prediction using XGBoost--------------
# you can do prediction using the following line
# you can put in Matrix, SparseMatrix or DMatrix
preds = predict(bst, test_X)
print("test-error=", sum((preds .> 0.5) .!= test_Y) / float(size(preds)[1]), "\n")

#-------------------save and load models-------------------------
# save model to binary local file
save(bst, "xgb.model")
# load binary model to julia
bst2 = Booster(model_file = "xgb.model")
preds2 = predict(bst2, test_X)
print("sum(abs(pred2-pred))=", sum(abs(preds2 .- preds)), "\n")

#----------------Advanced features --------------
# to use advanced features, we need to put data in xgb.DMatrix
dtrain = DMatrix(train_X, label = train_Y)
dtest = DMatrix(test_X, label = test_Y)

#---------------Using watchlist----------------
# watchlist is a list of DMatrix, each of them tagged with name
# DMatrix in watchlist should have label (for evaluation)
watchlist  = [(dtest,"eval"), (dtrain,"train")]
# we can change evaluation metrics, or use multiple evaluation metrics
bst = xgboost(dtrain, num_round, param=param, watchlist=watchlist, metrics=["logloss", "error"])

# we can also save DMatrix into binary file, then we can load it faster next time
save(dtest, "dtest.buffer")
save(dtrain, "dtrain.buffer")

# load model and data in
dtrain = DMatrix("dtrain.buffer")
dtest = DMatrix("dtest.buffer")
bst = Booster(model_file = "xgb.model")

# information can be extracted from DMatrix using get_info
label = get_info(dtest, "label")
pred = predict(bst, dtest)
print("test-error=", sum((pred .> 0.5) .!= label) / float(size(pred)[1]), "\n")

# You can dump the tree you learned using dump_model into a text file
dump_model(bst, "dump.raw.txt")
# If you have feature map file, you can dump the model in a more readable way
dump_model(bst, "dump.nice.txt", fmap="../data/featmap.txt", with_stats=true)

# You can also get information about feature importances in the model
dump(importance(bst, fmap="../data/featmap.txt"))


text:
function Input(texttodisplay)
    # function that allows for user input in the CLI
    println(texttodisplay)
    readline()
end

function runpage()
    options = "
    Choose model version to run:
    - 1 PAGE-ICE (10 timestep)
    - 2 PAGE-ICE with Growth Effects (10 timestep)
    - 3 PAGE-ICE (annual)
    - 4 PAGE-ICE with independent Interannual Temperature Variability (annual)
    - 5 PAGE-ICE with autoregressive Interannual Temperature Variability (annual)
    - 6 PAGE-ICE with Growth Effects (annual)
    - 7 PAGE-ICE with Growth Effects and autoregressive Interannual Temperature Variability (annual)

    Input number:
    "

    version = Input(options)

    if isequal(version, "1") || isequal(version, "PAGE-ICE (10 timestep)")
        include("../analysis/allscc/runmodel.jl")
    elseif isequal(version, "2") || isequal(version, "PAGE-ICE with Growth Effects (10 timestep)")
        include("../analysis/allscc/runmodel_growth.jl")
    elseif isequal(version, "3") || isequal(version, "PAGE-ICE (annual)")
        include("../analysis/allscc/runmodel_annual.jl")
    elseif isequal(version, "4") || isequal(version, "PAGE-ICE with independent Interannual Temperature Variability (annual)")
        include("../analysis/allscc/runmodel_variability.jl")
    elseif isequal(version, "5") || isequal(version, "PAGE-ICE with Interannual Temperature Variability with Autoregression (annual)")
        include("../analysis/allscc/runmodel_ARvariability.jl")
    elseif isequal(version, "6") || isequal(version, "PAGE-ICE with Growth Effects (annual)")
        include("../analysis/allscc/runmodel_annualGrowth.jl")
    elseif isequal(version, "7") || isequal(version, "PAGE-ICE with Growth Effects and autoregressive Interannual Temperature Variability (annual)")
        include("../analysis/allscc/runmodel_annualGrowth_ARvariability.jl")
    else
        println("WARNING: No valid model input provided. Please provide a valid model choice")
    end
end

runpage()


text:
using Distributed
using JLD2

include("../BioClim/src/required.jl")

## Get & prepare data
## Get & prepare data
@time @everywhere begin
    # Load data from CSV files
    #df = CSV.read("/Users/kiristern/Documents/GitHub/SDM/data/pred_prey/Puma_concolor.csv", header=true)
    # Prepare data (select columns, arrange values)
    #df = prepare_gbif_data(df)
    df = pred_prey
    # Separate species
    taxa_occ = [df[df.species .== u,:] for u in unique(df.species)]
    #
    # # Define coordinates range
    # lon_range = (-145.0, -50.0)
    # lat_range = (20.0, 75.0)
end

# ## Get the worldclim data
# @time wc_vars = pmap(x -> worldclim(x), 1:19);
# temp = wc_vars[1]

## Load predictions
@load "data/predictions.jld2" predictions
pred = predictions[1]
p1 = predictions[2]
p2 = predictions[3]
p3 = predictions[4]
p4 = predictions[5]
p5 = predictions[6]
p6 = predictions[7]
p7 = predictions[8]
# Combine layers
layers = [pred, p1, p2,p3,p4,p5,p6,p7]

# Create function
function expand_layers(layers::Array{SDMLayer{Float64},1})
    # Get minimum coordinates
    min_lon = min(map(x -> x.left, layers)...)
    max_lon = max(map(x -> x.right, layers)...)
    min_lat = min(map(x -> x.bottom, layers)...)
    max_lat = max(map(x -> x.top, layers)...)

    # Get grid size (rounding should prevent problems with last decimal)
    grid_size_lons = 1/round(1/stride(layers[1],1))
    grid_size_lats = 1/round(1/stride(layers[1],2))

    # Get coordinate range of newlayer -> original layers must have same stride
    lons_newlayers = min_lon+grid_size_lons:2*grid_size_lons:max_lon-grid_size_lons
    lats_newlayers = min_lat+grid_size_lats:2*grid_size_lats:max_lat-grid_size_lats

    # Create expanded layers
    newlayers = []
    for i in 1:length(layers)
        # Get position of original layers in new layer
        ## find where the point = 0 (= min value; all others are the absolute value)
        m_lon = findmin(abs.(layers[i].left+grid_size_lons .- lons_newlayers))[2]
        M_lon = findmin(abs.(layers[i].right-grid_size_lons .- lons_newlayers))[2]
        m_lat = findmin(abs.(layers[i].bottom+grid_size_lats .- lats_newlayers))[2]
        M_lat = findmin(abs.(layers[i].top-grid_size_lats .- lats_newlayers))[2]
        # Create new grid
        newgrid = fill(NaN, length(lats_newlayers), length(lons_newlayers))
        # Fill in original values
        newgrid[(m_lat:M_lat), (m_lon:M_lon)] .= layers[i].grid
        # Convert to SDMLayer
        newlayer = SDMLayer(newgrid, min_lon, max_lon, min_lat, max_lat)
        # Export result
        push!(newlayers, newlayer)
    end
    return newlayers
end

@time newlayers = expand_layers(layers)

# Plot result
plotn1 = plotSDM(newlayers[1])
plotn2 = plotSDM(newlayers[2])
plotn3 = plotSDM(newlayers[3])

# Combine heatmaps
sdm_plot = plotSDM(newlayers[1])
heatmap!(
    sdm_plot,
    longitudes(newlayers[5]), latitudes(newlayers[5]), # layer range
    newlayers[5].grid, # evenness values
    c=cgrad([:lightpurple, :purple, :darkpurple]), # ~color palette
    clim=(0.0, maximum(filter(!isnan, newlayers[5].grid)))) # colorbar limits
sdm_plot
# plot!(
#     sdm_plot,
#     xlims=(pred.left, pred.right),
#     ylims=(pred.bottom, pred.top),
#     aspectratio=92.60/60.75
# )
# Compare to p2
plotp1

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
function run(s)
    # Your code here
    return 0
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
using LanguageDynamics
using Documenter

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

makedocs(;
    modules=[LanguageDynamics],
    authors="Henri Kauhanen <mail@hkauhanen.fi>",
    repo="https://github.com/hkauhanen/LanguageDynamics.jl/blob/{commit}{path}#{line}",
    sitename="LanguageDynamics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://hkauhanen.github.io/LanguageDynamics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Getting started" => "getting-started.md",
        "Philosophy" => "philosophy.md",
        "Examples" => "examples.md",
        "Reference" => ["Populations" => "populations.md",
                        "Interactors" => "interactors.md",
                        "Replicators" => "replicators.md",
                        "Abstract types" => "abstract-types.md",
                        "Auxiliaries" => "auxiliaries.md",
                        "All functions and types" => "all.md",
                        "The type hierarchies" => "type-hierarchies.md"]
    ],
)

deploydocs(;
    repo="github.com/hkauhanen/LanguageDynamics.jl",
    devbranch="main",
)


text:
using Documenter
using CounterTools

makedocs(
    modules     = [CounterTools],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename    = "CounterTools",
    doctest     = true,
    pages       = Any[
        "Getting Started"   => "index.md",
        "Manual" => Any[
            "monitors.md",
            "records.md",
            "counters.md",
        ],
        "Examples" => Any[
            "core/example.md",
            "uncore/imc.md",
        ],
    ],
)

deploydocs(
    repo        = "github.com/hildebrandmw/CounterTools.jl.git",
    target      = "build",
)


text:
# ---
# title: 47. Permutations II
# id: problem47
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/permutations-ii/description/>
# hidden: true
# ---
#
# Given a collection of numbers, `nums`, that might contain duplicates, return
# _all possible unique permutations **in any order**._
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,1,2]
#     Output:
#     [[1,1,2],
#      [1,2,1],
#      [2,1,1]]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3]
#     Output: [[1,2,3],[1,3,2],[2,1,3],[2,3,1],[3,1,2],[3,2,1]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 8`
#   * `-10 <= nums[i] <= 10`
#
#
## @lc code=start
using LeetCode

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


text:
#=
Original code by Thomas Moore
@denbigh
included in https://github.com/ypaul21/Clapeyron.jl/pull/56
=#
"""
    DETPFlash(;numphases = 2;max_steps = 1e4*(numphases-1),population_size =50,time_limit = Inf)

Method to solve non-reactive multicomponent flash problem by finding global minimum of Gibbs Free Energy via Differential Evolution.

User must assume a number of phases, `numphases`. If true number of phases is smaller than numphases, model should predict either (a) identical composition in two or more phases, or (b) one phase with negligible total number of moles. If true number of phases is larger than numphases, a thermodynamically unstable solution will be predicted.

The optimizer will stop at `max_steps` evaluations or at `time_limit` seconds
"""
Base.@kwdef struct DETPFlash <: TPFlashMethod
    numphases::Int = 2
    max_steps::Int = 1e4*(numphases-1)
    population_size::Int = 50
    time_limit::Float64 = Inf
end

#z is the original feed composition, x is a matrix with molar fractions, n is a matrix with molar amounts
function partition!(dividers,n,x,nvals)
    numphases, numspecies = size(x)
    @inbounds for j = 1:numspecies
        nvals[1,j] = dividers[1,j]*n[j]
        for i = 2:numphases - 1
            Σ = sum(@view(nvals[1:(i-1), j]))
            nvals[i,j] = dividers[i,j]*(n[j] - Σ)
        end
        Σphases = sum(@view(nvals[1:(numphases-1), j]))
        nvals[numphases, j] = n[j] - Σphases
    end
#Calculate mole fractions xij
for i = 1:numphases
    ni =  @view(nvals[i, :])
    invn = 1/sum(ni)
    xi = @view(x[i, :])
    xi  .= ni
    xi .*= invn
end
end

function tp_flash_impl(model::EoSModel, p, T, n, method::DETPFlash)

    numspecies = length(model)
    TYPE = typeof(p+T+first(n))
    numphases = method.numphases
    x = zeros(TYPE,numphases, numspecies)
    nvals = zeros(TYPE,numphases, numspecies)

    GibbsFreeEnergy(dividers) = Obj_de_tp_flash(model,p,T,n,dividers,numphases,x,nvals)

    options = Metaheuristics.Options(time_limit = method.time_limit,iterations = method.max_steps,seed = UInt(373))
    algorithm = Metaheuristics.DE(N=method.population_size,options=options)

    #Minimize Gibbs Free Energy
#     result = bboptimize(GibbsFreeEnergy; SearchRange = (0.0, 1.0),
#         NumDimensions = numspecies*(numphases-1), MaxSteps=MaxSteps, PopulationSize = PopulationSize,
#         TraceMode = TraceMode)

    bounds = vcat(zeros(TYPE,(1,numspecies*(numphases-1))),ones(TYPE,1,numspecies*(numphases-1)))
    result = Metaheuristics.optimize(GibbsFreeEnergy,bounds,algorithm)

    dividers = reshape(Metaheuristics.minimizer(result),
            (numphases - 1, numspecies))

    #Initialize arrays xij and nvalsij,
    #where i in 1..numphases, j in 1..numspecies
    #xij is mole fraction of j in phase i.
    #nvals is mole numbers of j in phase i.
    partition!(dividers,n,x,nvals)

    return (x, nvals, Metaheuristics.minimum(result))
end
"""
    Obj_de_tp_flash(model,p,T,z,dividers,numphases)

Function to calculate Gibbs Free Energy for given partition of moles between phases.
This is a little tricky.

We must find a way of uniquely converting a vector of numbers,
each in (0, 1), to a partition. We must be careful that
the mapping is 1-to-1, not many-to-1, as if many inputs
map to the same physical state in a redundant way, there
will be multiple global optima, and the global optimization
will perform poorly.

Our approach is to specify (numphases-1) numbers in (0,1) for
each species. We then scale these numbers systematically in order to partition
the species between the phases. Each set of (numphases - 1) numbers
will result in a unique partition of the species into the numphases
phases.
"""
function Obj_de_tp_flash(model,p,T,n,dividers,numphases,x,nvals)
    _0 = zero(p+T+first(n))
    numspecies = length(n)
    TYPE  = typeof(_0)
    bignum = TYPE(1e300)
    dividers = reshape(dividers,
        (numphases - 1, numspecies))
    #Initialize arrays xij and nvalsij,
    #where i in 1..numphases, j in 1..numspecies
    #xij is mole fraction of j in phase i.
    #nvals is mole numbers of j in phase i.
    #x = zeros(TYPE,numphases, numspecies)
    #nvals = zeros(TYPE,numphases, numspecies)
    #Calculate partition of species into phases
    partition!(dividers,n,x,nvals)
    #Calculate Overall Gibbs Free Energy (J)
    #If any errors are encountered, return a big number, ensuring point is discarded
    #by DE Algorithm
    G = _0
    for i ∈ 1:numphases
            G += gibbs_free_energy(model, p, T, @view(nvals[i, :]))
            #calling with PTn calls the internal volume solver
            #if it returns an error, is a bug in our part.
    end
    return ifelse(isnan(G),bignum,G/R̄/T)
end

numphases(method::DETPFlash) = method.numphases

export DETPFlash


text:
module GrpcerrorsClients
using gRPCClient

include("grpcerrors.jl")
using .grpcerrors

import Base: show

# begin service: grpcerrors.GRPCErrors

export GRPCErrorsBlockingClient, GRPCErrorsClient

struct GRPCErrorsBlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::GRPCErrorsBlockingStub

    function GRPCErrorsBlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = GRPCErrorsBlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct GRPCErrorsClient
    controller::gRPCController
    channel::gRPCChannel
    stub::GRPCErrorsStub

    function GRPCErrorsClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = GRPCErrorsStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::GRPCErrorsBlockingClient) = print(io, "GRPCErrorsBlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::GRPCErrorsClient) = print(io, "GRPCErrorsClient(", client.channel.baseurl, ")")

import .grpcerrors: SimpleRPC
"""
    SimpleRPC

- input: grpcerrors.Data
- output: grpcerrors.Data
"""
SimpleRPC(client::GRPCErrorsBlockingClient, inp::grpcerrors.Data) = SimpleRPC(client.stub, client.controller, inp)
SimpleRPC(client::GRPCErrorsClient, inp::grpcerrors.Data, done::Function) = SimpleRPC(client.stub, client.controller, inp, done)

import .grpcerrors: StreamResponse
"""
    StreamResponse

- input: grpcerrors.Data
- output: Channel{grpcerrors.Data}
"""
StreamResponse(client::GRPCErrorsBlockingClient, inp::grpcerrors.Data) = StreamResponse(client.stub, client.controller, inp)
StreamResponse(client::GRPCErrorsClient, inp::grpcerrors.Data, done::Function) = StreamResponse(client.stub, client.controller, inp, done)

import .grpcerrors: StreamRequest
"""
    StreamRequest

- input: Channel{grpcerrors.Data}
- output: grpcerrors.Data
"""
StreamRequest(client::GRPCErrorsBlockingClient, inp::Channel{grpcerrors.Data}) = StreamRequest(client.stub, client.controller, inp)
StreamRequest(client::GRPCErrorsClient, inp::Channel{grpcerrors.Data}, done::Function) = StreamRequest(client.stub, client.controller, inp, done)

import .grpcerrors: StreamRequestResponse
"""
    StreamRequestResponse

- input: Channel{grpcerrors.Data}
- output: Channel{grpcerrors.Data}
"""
StreamRequestResponse(client::GRPCErrorsBlockingClient, inp::Channel{grpcerrors.Data}) = StreamRequestResponse(client.stub, client.controller, inp)
StreamRequestResponse(client::GRPCErrorsClient, inp::Channel{grpcerrors.Data}, done::Function) = StreamRequestResponse(client.stub, client.controller, inp, done)

# end service: grpcerrors.GRPCErrors

end # module GrpcerrorsClients


text:
tempSum = 0
numOfCells = 0

function countAvgTemp(r::Robot)
    dir = Ost
    move_while_can!(r, dir)
    dir = inverse(dir)
    while !isborder(r,Nord)
      move!(r, Nord)
      move_while_can!(r, dir)
      dir = inverse(dir)
    end
    if numOfCells > 0
        return tempSum/numOfCells
    else
        return 0
    end
end

#Двигает робота в сторону side, пока нет перегородки
function move_while_can!(r::Robot, side)
    checkCell(r)
    while !isborder(r, side)
        move!(r, side)
        checkCell(r)
    end
end

function checkCell(r::Robot)
   global tempSum, numOfCells
   if ismarker(r)
      tempSum += temperature(r)
      numOfCells += 1
   end
end

function inverse(side)
    return HorizonSide(mod(Int(side) + 2, 4))
end


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--gpu"
            help = "use gpu"
            action = :store_true
        "--iter"
            help = "number of iterations"
            arg_type = Int
            default = 1000
        "--step"
            help = "interval of checking monitored values"
            arg_type = Int
            default = 100
        "--seed"
            help = "seed"
            arg_type = Int
            default = 777
        "--Float32"
            help = "use Float32 instead of Float64"
            action = :store_true
        "--set_zero_subnormals"
            help = "set subnormals to zero"
            action = :store_true
        "--init_from_master"
            help = "use centralized random initialization (costly)"
            action = :store_true
        "--eval_obj"
            help = "evaluate objective function (costly). Maximum difference from the previous iteration is printed otherwise."
            action = :store_true
    end

    return s
end

function parse_commandline_nmf()
    s = parse_commandline()
    @add_arg_table s begin
        "--rows"
            help = "number of rows"
            arg_type = Int
            default = 10_000
        "--cols"
            help = "number of cols"
            arg_type = Int
            default = 10_000
        "--r"
            help = "intermediate size"
            arg_type = Int
            default = 20
    end
    return parse_args(s)
end

function parse_commandline_mds()
    parse_commandline_nmf()
end

function parse_commandline_cox_base()
    s = parse_commandline()
    @add_arg_table s begin
        "--lambda"
            help = "regularization parameter"
            arg_type = Float64
            default=0.0
    end
    return s
end

function parse_commandline_cox()
    s = parse_commandline_cox_base()
    @add_arg_table s begin
        "--rows"
            help = "number of rows"
            arg_type = Int
            default = 10_000
        "--cols"
            help = "number of cols"
            arg_type = Int
            default = 10_000
        "--censor_rate"
            help = "rate of censored subjects"
            arg_type = Float64
            default=0.5
    end

    return parse_args(s)
end

function parse_commandline_cox_ukbk()
    s = parse_commandline_cox_base()
    @add_arg_table s begin
        "--prefix"
            help = "prefix for the output file"
            arg_type = String
            default = "test"
    end
    return parse_args(s)
end

function parse_commandline_pet()
    s = parse_commandline()
    @add_arg_table s begin
        "--data"
            help = "path to datafile"
            arg_type = String
        "--reg"
            help = "regularization parameter"
            arg_type = Float64
            default=0.0
    end
    return s
end

function parse_commandline_pet_l2()
    s = parse_commandline_pet()
    return parse_args(s)
end
function parse_commandline_pet_l1()
    s = parse_commandline_pet()
    @add_arg_table s begin
        "--sigma"
            help = "dual stepsize"
            arg_type = Float64
            default=1/3
        "--tau"
            help = "primal stepsize"
            arg_type = Float64
            default=1/3
    end
    return parse_args(s)
end


text:
module Atom

include("common.jl")
export UP, DOWN
export Operator
export Float

include("hilbert.jl")
export Hilbert
# export Basis.Hilbert, Basis.Fock
# export Basis.BinaryFock, Basis.creation

include("green.jl")
export Green

include("hubbard.jl")
export Hubbard

end # module


text:
using JuMP
using Ipopt

function Solve_Dist_SOC(buses, lines, generators)
info("Starting DistSOC Model")

bus_set = collect(keys(buses))
line_set = collect(keys(lines))
gen_set = collect(keys(generators))
gen_bus = []
for g in keys(generators)
    push!(gen_bus, generators[g].bus_idx)
end
root_bus = 0
for b in keys(buses)
    if buses[b].is_root
        root_bus = b
    end
end

lines_to = Dict()
for l in keys(lines)
    lines_to[lines[l].to_node] = lines[l]
end

v_root = 1

m = Model(solver=IpoptSolver(print_level=3))

@variable(m, v[bus_set] >= 0) #variable for voltage square
@variable(m, fp[bus_set]) #variable for active power flow
@variable(m, fq[bus_set]) # variable for reactive power flow
@variable(m, gp[bus_set]) #variable for active power generation
@variable(m, gq[bus_set]) #variable for reactive power generation
@variable(m, l[bus_set] >= 0) # variable for quadratic current

@objective(m, Min, sum((gp[b])*buses[b].generator.cost for b in gen_bus))

for b in gen_bus
# All buses with generation
    @constraint(m, gp[b] <= buses[b].generator.g_P_max) #Upper limit constraint for active power generation
    @constraint(m, gp[b] >= 0) #Lower limit constraint for active power generation
    @constraint(m, gq[b] <= buses[b].generator.g_Q_max) #Upper limit constraint for reactive power generation
    @constraint(m, gq[b] >= -buses[b].generator.g_Q_max) #Lower limit constraint for reactive power generation
end

for b in setdiff(bus_set, gen_bus)
# All buses without generation
    @constraint(m, gp[b] == 0)
    @constraint(m, gq[b] == 0)
end

for b in bus_set
# All buses
    @constraint(m, buses[b].d_P - gp[b] + sum((fp[k] + l[k]*lines_to[k].r) for k in buses[b].children) == fp[b])
    @constraint(m, buses[b].d_Q - gq[b] + sum((fq[k] + l[k]*lines_to[k].x) for k in buses[b].children) == fq[b])

    @constraint(m, v[b] <= buses[b].v_max)
    @constraint(m, v[b] >= buses[b].v_min)
end

for b in setdiff(bus_set, [root_bus])
# All buses without root
    b_ancestor = buses[b].ancestor[1]
    R = lines_to[b].r
    X = lines_to[b].x
    @constraint(m, v[b] == v[b_ancestor] - 2*(R*fp[b] + X*fq[b]) - l[b]*(R+X))
    @constraint(m, fp[b]^2 + fq[b]^2 <= lines_to[b].s_max^2)
    @constraint(m, (fp[b] + l[b]*R)^2 + (fq[b] + l[b]*X)^2 <= lines_to[b].s_max^2)
    @NLconstraint(m, (fp[b]^2 + fq[b]^2)/l[b] <= v[b])
end

#Voltage and flow constraints for root node
@constraint(m, v[root_bus] == v_root)
@constraint(m, l[root_bus] == 0)
@constraint(m, fp[root_bus] == 0)
@constraint(m, fq[root_bus] == 0)

solve(m)

# Prepare Results
bus_results = DataFrame(bus = Any[], d_P = Any[], d_Q = Any[], g_P = Any[], g_Q = Any[], v_squared = Any[])
line_results = DataFrame(line = Any[], from = Any[], to = Any[], f_P = Any[], f_Q = Any[], a_squared = Any[])
for b in bus_set
    res = [b, buses[b].d_P, buses[b].d_Q, getvalue(gp[b]), getvalue(gq[b]), getvalue(v[b])]
    push!(bus_results, res)

    # Calc square current on line
    a_res = 0
    fp_res = getvalue(fp[b])
    fq_res = getvalue(fq[b])
    if b != root_bus
        a_res = getvalue(l[b])
        lres = [lines_to[b].index, buses[b].ancestor[1], b, fp_res, fq_res, a_res]
        push!(line_results, lres)
    end
end

sort!(bus_results, cols=:bus)
sort!(line_results, cols=:to)

return getobjectivevalue(m), bus_results, line_results

end


text:
"""
$(TYPEDEF)

Central difference scheme with given step-size.

### Fields

- `Δt` -- step-size
"""
struct CentralDifference{N} <: AbstractSolver
    Δt::N
end

CentralDifference(; Δt::N) where N = CentralDifference(Δt)

function _init(alg::CentralDifference)
    Δt = alg.Δt

    # compute integration constants
    a₀ = 1/Δt^2
    a₁ = 1/(2*Δt)
    a₂ = 2*a₀
    a₃ = 1/a₂

    return a₀, a₁, a₂, a₃
end

function _build_solution(alg::CentralDifference{N}, U, NSTEPS) where {N}
    t = range(zero(N), step=alg.Δt, length=(NSTEPS+1))
    return Solution(alg, U, nothing, nothing, t)
end

function _solve(alg::CentralDifference{N},
                ivp::InitialValueProblem{ST, XT},
                NSTEPS::Int; kwargs...) where {N, VT, ST, XT<:Tuple{VT, VT}}

    sys = system(ivp)
    (U₀, U₀′) = initial_state(ivp)
    Δt = alg.Δt
    IMAX = NSTEPS + 1
    M, C, K, R = _unwrap(sys, IMAX)

    U₀′′ = M \ (R[1] - C * U₀′ - K * U₀)

    # integration constants
    a₀, a₁, a₂, a₃ = _init(alg)
    U⁻ = U₀ - Δt * U₀′ + a₃ * U₀′′
    M̂ = a₀*M + a₁*C
    M̂⁻¹ = factorize(M̂)
    M̂₋ = a₀*M - a₁*C

    U = Vector{VT}(undef, IMAX+1)
    U[1] = U⁻
    U[2] = U₀

    @inbounds for i in 2:IMAX
        R̂ᵢ = R[i] - (K - a₂ * M) * U[i] - M̂₋ * U[i-1]
        U[i+1] = M̂⁻¹ \ R̂ᵢ
    end

    return _build_solution(alg, view(U, 2:IMAX+1), NSTEPS)
end

function _solve(alg::CentralDifference{N},
                ivp::InitialValueProblem{ST, XT},
                NSTEPS::Int;
                t0=zero(N), kwargs...) where {N, VT, ST<:SecondOrderContinuousSystem, XT<:Tuple{VT, VT}}

    sys = system(ivp)
    (U₀, U₀′) = initial_state(ivp)
    Δt = alg.Δt
    IMAX = NSTEPS + 1
    @unpack M, C, fi, fe = sys

    # initial acceleration, obtained from the equilibrium
    # condition Mu'' + Cu' + fint(u) = fext(t), at time t = t0
    U₀′′ = M \ (fe(t0) - C * U₀′ - fi(U₀))

    # integration constants
    a₀, a₁, a₂, a₃ = _init(alg)
    U⁻ = U₀ - Δt * U₀′ + a₃ * U₀′′
    M̂ = a₀*M + a₁*C
    M̂⁻¹ = factorize(M̂)
    M̂₋ = a₀*M - a₁*C

    U = Vector{VT}(undef, IMAX+1)
    U[1] = U⁻
    U[2] = U₀

    t = t0
    @inbounds for i in 2:IMAX
        R̂ᵢ = fe(t) - fi(U[i]) + a₂ * M * U[i] - M̂₋ * U[i-1]
        U[i+1] = M̂⁻¹ \ R̂ᵢ
        t += Δt
    end

    return _build_solution(alg, view(U, 2:IMAX+1), NSTEPS)
end


text:
module RelPerm

module Corey
	function S_wn(S_w, S_wir, S_nir)
		return (S_w - S_wir) / (1 - S_wir - S_nir)
	end
	function kr_w(S_w, S_wir, S_nir, N_w, k0_w)
		if S_w < S_wir
			return 0.0
		elseif S_w > 1 - S_nir
			return k0_w
		else
			return k0_w * (S_wn(S_w, S_wir, S_nir)) ^ N_w
		end
	end
	function kr_n(S_n, S_wir, S_nir, N_n, k0_n)
		if S_n < S_nir
			return 0.0
		elseif S_n > 1 - S_wir
			return k0_n
		else
			return k0_n * (1 - S_wn(1 - S_n, S_wir, S_nir)) ^ N_n
		end
	end
end

end


text:
# # Spacecraft
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Spacecraft.ipynb)
#
#md # !!! note "Overview"
#md #     System type: hybrid system with nonlinear dynamics\
#md #     State dimension: 4 + 1\
#md #     Application domain: Flight Dynamics, Switched Controllers
#
# ## Model description
#
# Spacecraft rendezvous is a perfect use case for formal veriﬁcation of hybrid
# systems with nonlinear dynamics since mission failure can cost lives and is
# extremely expensive. This benchmark is taken from [1].

# The nonlinear dynamic equations describe the two-dimensional, planar motion of
# the space-craft on an orbital plane towards a space station:
#
# ```math
#   \left\{ \begin{array}{lcl}
#   \dot{x} & = & v_x \\
#   \dot{y} & = & v_y \\
#   \dot{v_x} & = & n^2x + 2nv_y + \frac{\mu}{r^2} - \frac{\mu}{r^3} (r +x) + \frac{u_x}{m_c} \\
#   \dot{v_y} & = & n^2y - 2nv_x - \frac{\mu}{r^3_c}y + \frac{u_y}{m_c}
#   \end{array} \right.
# ```
#
# The model consists of position (relative to the target) ``x``, ``y`` [m], time
# ``t`` [min], as well as horizontal and vertical velocity ``v_x``, ``v_y`` [m / min].
# The parameters are ``µ = 3.986 × 10^{14} × 60^2`` [``m^3 / min^2``],
# ``r = 42164 × 10^3 [m], m_c = 500 [kg], n = \sqrt{\frac{\mu}{r^3}}`` and
# ``r_c = \sqrt{(r + x)^2 + y^2}``.
#
# The hybrid nature of this benchmark originates from a switched controller. In
# particular, the modes are approaching (``x ∈ [−1000, −100] [m]``), rendezvous
# attempt (``x ≥ −100 [m]``), and aborting. A transition to mode aborting occurs
# nondeterministically at ``t ∈ [120, 150] [min]``.
#
# The linear feedback controllers for the diﬀerent modes are deﬁned as
# ``\binom{u_x}{u_y} = K_1\underline{x}`` for mode approaching, and
# ``\binom{u_x}{u_y} = K_2\underline{x}`` for mode rendezvous attempt, where
# ``\underline{x} = (x, y, v_x, v_y)^T`` is the vector of system states. The
# feedback matrices ``K_i`` were determined with an LQR-approach applied to the
# linearized system dynamics, which resulted in the following numerical values:
# ```math
# K_1 =
# \begin{pmatrix}
# −28.8287 & 0.1005 & −1449.9754 & 0.0046 \\
# −0.087 & −33.2562 & 0.00462 & −1451.5013
# \end{pmatrix}
# ```
# ```math
# K_2 =
# \begin{pmatrix}
# −288.0288 & 0.1312 & −9614.9898 & 0 \\
# −0.1312 & −288 & 0 & −9614.9883
# \end{pmatrix}
# ```
# In the mode aborting, the system is uncontrolled ``\binom{u_x}{u_y} = \binom{0}{0}``.

using ReachabilityAnalysis

const μ = 3.986e14 * 60^2
const r = 42164.0e3
const r² = r^2
const mc = 500.0
const n² = μ / r^3
const n = sqrt(n²)

const two_n = 2*n
const μ_r² = μ/r²;

#-
const K₁ = [-28.8287 0.1005 -1449.9754 0.0046 0.0;
            -0.087 -33.2562 0.00462 -1451.5013 0.0]
const K₂ = [-288.0288 0.1312 -9614.9898 0.0 0.0;
            -0.1312 -288.0 0.0 -9614.9883   0.0]

const K₁mc = K₁/mc
const K₂mc = K₂/mc

function mymul!(v, A, x)
    @inbounds for ind = 1:length(v)
        v[ind] = zero(x[1])
        for jind = 1:length(x)
            v[ind] += A[ind, jind] * x[jind]
        end
    end
    return nothing
end

# ### Dynamics in the 'approaching' mode

@taylorize function spacecraft_approaching!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    uxy = Vector{typeof(x)}(undef, 2)
    mymul!(uxy, K₁mc, u)

    ## x' = vx
    du[1] = vx

    ## y' = vy
    du[2] = vy

    ## vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    du[3] = (n²*x + two_n*vy) + ((μ_r² - μ_rc³*rx) + uxy[1])

    ## vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    du[4] = (n²*y - two_n*vx) - (μ_rc³*y - uxy[2])

    ## t' = 1
    du[5] = one(x)

    return du
end

# ### Dynamics in the 'rendezvous attempt' mode

@taylorize function spacecraft_attempt!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    uxy = Vector{typeof(x)}(undef, 2)
    mymul!(uxy, K₂mc, u)

    ## x' = vx
    du[1] = vx

    ## y' = vy
    du[2] = vy

    ## vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    du[3] = (n²*x + two_n*vy) + ((μ_r² - μ_rc³*rx) + uxy[1])

    ## vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    du[4] = (n²*y - two_n*vx) - (μ_rc³*y - uxy[2])

    ## t' = 1
    du[5] = one(x)

    return du
end

# ### Dynamics in the 'aborting' mode

@taylorize function spacecraft_aborting!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    ## x' = vx
    du[1] = vx

    ## y' = vy
    du[2] = vy

    ## vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x)
    du[3] = (n²*x + two_n*vy) + (μ_r² - μ_rc³*rx)

    ## vy' = n²y - 2n*vx - μ/(rc^3)y
    du[4] = (n²*y - two_n*vx) - μ_rc³*y

    ## t' = 1
    du[5] = one(x)

    return du
end

# ### Hybrid System setup

# To model the system as a hybrid automaton, it is useful to work with symbolic
# state variables.

using Symbolics

const var = @variables x y vx vy t

function spacecraft(; abort_time=(120.0, 150.0))

    n = 4 + 1  # number of variables
    t_abort_lower, t_abort_upper = abort_time

    automaton = LightAutomaton(3)

    ## mode 1 "approaching"
    invariant = HalfSpace(x <= -100, var)
    approaching = @system(x' = spacecraft_approaching!(x), dim:5, x ∈ invariant)

    ## mode 2 "rendezvous attempt"
    invariant = HalfSpace(x >= -100, var)
    attempt = @system(x' = spacecraft_attempt!(x), dim:5, x ∈ invariant)

    ## mode 3 "aborting"
    invariant = Universe(n)
    aborting = @system(x' = spacecraft_aborting!(x), dim: 5, x ∈ invariant)

    ## transition "approach" -> "attempt"
    add_transition!(automaton, 1, 2, 1)
    guard = HalfSpace(x >= -100, var)
    t1 = @map(x -> x, dim: n, x ∈ guard)

    ## transition "approach" -> "abort"
    add_transition!(automaton, 1, 3, 2)
    guard_time = HPolyhedron([t >= t_abort_lower, t <= t_abort_upper], var)
    t2 = @map(x -> x, dim: n, x ∈ guard_time)

    ## transition "attempt" -> "abort"
    add_transition!(automaton, 2, 3, 3)
    t3 = @map(x -> x, dim: n, x ∈ guard_time)

    H = HybridSystem(automaton=automaton,
                     modes=[approaching, attempt, aborting],
                     resetmaps=[t1, t2, t3])

     ## initial condition in mode 1
     X0 = Hyperrectangle([-900., -400., 0., 0., 0.], [25., 25., 0., 0., 0.])
     init = [(1, X0)]

     return InitialValueProblem(H, init)
end

# ## Reachability settings
#
# The spacecraft starts from the initial set ``x ∈ [−925, −875] [m]``,
# ``y ∈ [−425, −375] [m], vx = 0 [m/min]`` and ``vy = 0 [m/min]``. For the
# considered time horizon of ``t ∈ [0, 200] [min]``, the following speciﬁcations
# have to be satisﬁed:
#
# - Line-of-sight: In mode rendezvous attempt, the spacecraft has to stay
#      inside line-of-sight cone ``L = {\binom{x}{y} | (x ≥ −100) ∧ (y ≥ x \tan(30°)) ∧ (−y ≥ x \tan(30°))}``.
#
# - Collision avoidance: In mode aborting, the spacecraft has to avoid a
#      collision with the target, which is modeled as a box ``B`` with ``0.2m``
#      edge length and the center placed at the origin.
#
# - Velocity constraint: In mode rendezvous attempt, the absolute velocity
#      has to stay below ``3.3 [m/min]``: ``\sqrt{v^2_x + v^2_y} ≤ 3.3 [m/min]``.
#
# !!! note "Remark on velocity constraint"
#     In the original benchmark [1], the constraint
#     on the velocity was set to 0.05 m/s, but it can be shown (by a counterexample)
#     that this constraint cannot be satisﬁed. We therefore use the relaxed
#     constraint ``0.055 [m/s] = 3.3 [m/min]``.

const tan30 = tand(30)

LineOfSightCone = HPolyhedron([x >= -100, y >= x*tan30, -y >= x*tan30], var)

Target = Hyperrectangle(zeros(2), [0.2, 0.2]);

# ### Auxiliary functions

function line_of_sight(sol)
    all_idx = findall(x -> x == 2, location.(sol))  # attempt
    for idx in all_idx
        verif = all(set(R) ⊆ LineOfSightCone for R in sol[idx][2:end])
        !verif && return false
    end
    return true
end

function absolute_velocity(R)
    vx, vy = 3, 4
    vx2 = set(overapproximate(project(R, vars=(vx,)), Interval))
    vy2 = set(overapproximate(project(R, vars=(vy,)), Interval))
    return sqrt(vx2.dat + vy2.dat)
end

function velocity_constraint(sol)
    all_idx = findall(x -> x == 2, location.(sol)) # attempt mode
    for idx in all_idx
        ## maximum velocity measured in m/min
        verif = all(absolute_velocity(R) < 0.055 * 60. for R in sol[idx])
        !verif && return false
    end
    return true
end

function target_avoidance(sol)
    all_idx = findall(x -> x == 3, location.(sol)) # aborting mode
    for idx in all_idx
        verif = all(is_intersection_empty(set(Projection(R, [x, y])), Target) for R in sol[idx])
        !verif && return false
    end
    return true
end

# ## Results

# We used ``n_Q = 1`` and ``n_T = 5``. The transition to the aborting mode is
# handled by clustering the flowpipe into 29 boxes.

using Plots

function solve_spacecraft(prob; k=25, s=missing)

    ## template directions
    boxdirs = BoxDirections(5)

    ## transition from mode 1 to mode 2
    sol12 = solve(prob,
                tspan=(0.0, 200.0),
                alg=TMJets20(abstol=1e-5, maxsteps=10_000, orderT=5, orderQ=1, disjointness=BoxEnclosure()),
                max_jumps=1,
                intersect_source_invariant=false,
                intersection_method=TemplateHullIntersection(boxdirs),
                clustering_method=LazyClustering(1),
                disjointness_method=BoxEnclosure())

    sol12jump = overapproximate(sol12[2](120 .. 150), Zonotope)
    t0 = tstart(sol12jump[1])
    sol12jump_c = cluster(sol12jump, 1:length(sol12jump), BoxClustering(k, s))

    ## transition from mode 2 to mode 3
    H = system(prob)
    sol3 = solve(IVP(mode(H, 3), [set(X) for X in sol12jump_c]),
                 tspan=(t0, 200.0),
                 alg=TMJets20(abstol=1e-10, orderT=7, orderQ=1, disjointness=BoxEnclosure()))
    d = Dict{Symbol, Any}(:loc_id => 3)

    F12 = [fp for fp in sol12.F]
    F23 = [Flowpipe(fp.Xk, d) for fp in sol3.F]
    return HybridFlowpipe(vcat(F12, F23))
end

prob = spacecraft()
sol = solve_spacecraft(prob; k=25, s=missing)
solz = overapproximate(sol, Zonotope);

#----------------

plot(solz, vars=(1, 2), xlab="x", ylab="y")

#-

# Note that in the previous plot, it is not specified which reach-set corresponds to
# which mode. After finding the location (mode of the hybrid system) associated
# to each reach-set of the solution, we can plot each one using different colors.
# The plotting code is more sophisticated than before, but is an example of finer
# control of the options used to plot the flowpipes.

idx_approaching = findall(x -> x == 1, location.(solz))
idx_attempt = findall(x -> x == 2, location.(solz))
idx_aborting = findall(x -> x == 3, location.(solz))

using Plots, Plots.PlotMeasures, LaTeXStrings

fig = plot(legend=:bottomright, tickfont=font(10, "Times"), guidefontsize=15,
           xlab=L"x", ylab=L"y", lw=0.0,
           xtick=[-750, -500, -250, 0, 250.], ytick=[-400, -300, -200, -100, 0.],
           xlims=(-1000.0, 400.0), ylims=(-450.0, 0.0), size=(600, 600),
           bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)

plot!(fig, solz[idx_approaching[1]], vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1., lab="Approaching")
for k in 2:length(idx_approaching)
    plot!(fig, solz[idx_approaching[k]], vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1.)
end

plot!(fig, solz[idx_attempt[1]], vars=(1, 2), lw=0.0, color=:red, alpha=1., lab="Attempt")
for k in 2:length(idx_attempt)
    plot!(fig, solz[idx_attempt[k]], vars=(1, 2), lw=0.0, color=:red, alpha=1.)
end

plot!(fig, solz[idx_aborting[1]], vars=(1, 2), lw=0.0, color=:cyan, alpha=1., lab="Aborting")
for k in 2:length(idx_aborting)
    plot!(fig, solz[idx_aborting[k]], vars=(1, 2), lw=0.0, color=:cyan, alpha=1.)
end

fig

# ## References

# [1] N. Chan and S. Mitra. *Verifying safety of an autonomous spacecraft
# rendezvous mission.* In ARCH17. 4th International Workshop on Applied
# Veriﬁcation of Continuous and Hybrid Systems, collocated with Cyber-Physical
# Systems Week (CPSWeek) on April 17, 2017 in Pittsburgh, PA, USA, pages 20–32, 2017.
#


text:
function isHappy(s)::Bool
    a::Int = s isa Number ? s : (try
                                 parse(Int64, s)
                                 catch e
                                 throw(error("Wrong type $e"))
                                 end)
    innerHappy(a, Set([]))
end

function eachBit(n::Int)
    result = []
    while n != 0
        let left = n % 10
            push!(result, left)
            n = Int((n-left)/10)
        end
    end
    result
end

function innerHappy(n,se)
    n in se && return false
    next = sum(eachBit(n) .^ 2)
    next == 1 ? (return true) : (se = union(se, n))
    innerHappy(next,se)
end

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
data(f::Taylor1Model) = (f.n, f.x0, f.I)

function +(f::Taylor1Model, g::Taylor1Model)
    @assert data(f) == data(g)

    return Taylor1Model(f, f.p + g.p, f.Δ + g.Δ)
end

function -(f::Taylor1Model, g::Taylor1Model)
    @assert data(f) == data(g)

    return Taylor1Model(f, f.p - g.p, f.Δ - g.Δ)
end

+(α::Real, f::Taylor1Model) = Taylor1Model(f, f.p + α, f.Δ)
+(f::Taylor1Model, α::Real)= α + f

+(f::Taylor1Model, α::Real)= α + f

-(f::Taylor1Model) = (-1 * f)
-(α::Real, f::Taylor1Model) = α + (-f)
-(f::Taylor1Model, α::Real) = f + (-α)

*(α::Real, f::Taylor1Model) = Taylor1Model(f, α*f.p, α*f.Δ)
*(f::Taylor1Model, α::Real) = α * f

function *(f::Taylor1Model, g::Taylor1Model)
    @assert data(f) == data(g)

    n, x0, I = f.n, f.x0, f.I

    c = zeros(eltype(g.p), 2n+1)

    a = f.p
    b = g.p

    for i in 0:n
        for j in 0:n
            c[i+j+1] += a[i] * b[j]
        end
    end

    d = zeros(eltype(g.p), 2n+1)

    d[n+2:end] = c[n+2:end]

    B = bound(Taylor1(d), x0, I)
    Bf = bound(a, x0, I)
    Bg = bound(b, x0, I)

    Δ = B + (f.Δ * Bg) + (g.Δ * Bf) + (f.Δ * g.Δ)

    return Taylor1Model(f, Taylor1(c[1:n+1]), Δ)

end