text:
mutable struct PhotometryStructure
           pokes::DataFrames.AbstractDataFrame
           streaks::DataFrames.AbstractDataFrame
           traces::DataFrames.AbstractDataFrame
       end

"""
`verify_names`
check that all the Dataframes in an Array of PhotometryStructures
have the same names for a given field
"""
function verify_names(data::Array{Flipping.PhotometryStructure,1},field::Symbol)
   verify = Array{Bool,1}(0)
   for i = 2:size(data,1)
       first = names(getfield(data[i], field))
       second = names(getfield(data[i-1], field))
       push!(verify,first != second)
   end
   return findall(verify)
end

"""
`convertin_DB`
convert an Array of PhotometryStructure to a Indexed table
"""
function convertin_DB(data::Array{Flipping.PhotometryStructure},field::Symbol)
    provisory = []
    for i = 1:size(data,1)
        ongoing = JuliaDB.table(getfield(data[i],field))
        if isempty(provisory)
            provisory = ongoing
        else
            provisory = merge(provisory,ongoing)
        end
    end
    return provisory
end
#= old version convertin_DB
provisory = DataFrame()
for i = 1:size(data,1)
    t = deepcopy(getfield(data[i],field))
    if isempty(provisory)
        provisory = t
    else
        append!(provisory,t)
    end
end
converted = JuliaDB.table(provisory)
return converted
=#

function Base.unique(data::Array{Flipping.PhotometryStructure,1},column::Symbol, field::Symbol)
    result = []
    for i=1:size(data,1)
        session = getfield(data[i], field)
        values = unique(session[column])
        for value in values
            push!(result,value)
        end
    end
    result = unique(result)
end

function Base.maximum(data::Array{Flipping.PhotometryStructure,1},column::Symbol, field::Symbol)
    result = []
    for i=1:size(data,1)
        session = getfield(data[i], field)
        value = maximum(session[column])
        push!(result,value)
    end
    result = maximum(result)
end

function Base.minimum(data::Array{Flipping.PhotometryStructure,1},column::Symbol, field::Symbol)
    result = []
    for i=1:size(data,1)
        session = getfield(data[i], field)
        value = minimum(session[column])
        push!(result,value)
    end
    result = minimum(result)
end

function Base.copy(data::Flipping.PhotometryStructure)
    PhotometryStructure(copy(data.pokes,data.streaks,data.traces))
end


text:
mults = [1:1000]
@time  res = +(filter(x -> (x % 3 == 0) | (x % 5 == 0), mults)...)
res = 0
@time for m in mults
    if (m % 3 == 0) | (m % 5 == 0)
        res += m
    end
end
println("The answer to Euler Problem 1 is $res")

ext: jl

max_stars_repo_path: Problem 1.jl

max_stars_repo_name: andrew-christianson/Polyglot-Euler

max_issues_repo_path: Problem 1.jl

max_issues_repo_name: andrew-christianson/Polyglot-Euler

max_forks_repo_path: Problem 1.jl

max_forks_repo_name: andrew-christianson/Polyglot-Euler


text:
function detect_anagrams(subject::AbstractString, candidates::AbstractArray)

end


text:
using FractalDevelopment
using Test

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

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
"""
```
type Observable
```

### Fields

- `key::Symbol`
- `input_series::Vector{Symbol}`: vector of mnemonics, each in the form
  `:MNEMONIC__SOURCE` (e.g. `:GDP__FRED`). This vector is parsed to determine
  source (e.g. per-capita consumption gets population and consumption).
- `fwd_transform::Function`: Extracts appropriate `input_series` from a
  DataFrame of levels, and transforms data to model units (for example, computes
  per-capita growth rates from levels).
- `rev_transform::Function`: Transforms a series from model units into
  observable units. May take kwargs.
- `name::String`: e.g. \"Real GDP growth\"
- `longname::String`: e.g. \"Real GDP growth per capita\"
"""
type Observable
    key
    input_series::Vector{Symbol} # (vector of Mnemonics, e.g. GDPC1@FRED)
                                 # parse this to determine source
                                 # (eg consumption per cap gets population and consumption)

    fwd_transform::Function      # corresponds to data_transforms
                                 # field. Operates on dataframe of levels, that has every data series
                                 # requested by any observable

    rev_transform::Function      # corresponds to getyp, or getyp4q. May take kwargs (untransformed,
                                 # 1q transforms, 4q transforms, levels/growth rates, etc)

    name::String
    longname::String
end

"""
```
type PseudoObservable
```

### Fields

- `key::Symbol`
- `name::String`: e.g. \"Flexible Output Growth\"
- `longname::String`: e.g. \"Output that would prevail in a flexible-price
  economy\"
- `rev_transform::Function`: Transforms a series from model units into
  observable units. May take kwargs.
"""
type PseudoObservable
    key::Symbol
    name::String
    longname::String
    rev_transform::Function
end

function PseudoObservable(k::Symbol)
    PseudoObservable(k, "", "", identity)
end

function check_mnemonics(levels::DataFrame, mnemonics::Symbol)
    for mnemonic in mnemonics
        @assert in(mnemonic, names(levels)) "Dataframe is missing $(mnemonic)"
    end
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+-1*a*cos(x)^2)^-1*sin(x)^6, x) == :(3//8*x*a^-1+-3//8*a^-1*cos(x)*sin(x)+-1//4*a^-1*sin(x)^3*cos(x))
@test integrate((a+-1*a*cos(x)^2)^-1*sin(x)^5, x) == :(-1*a^-1*cos(x)+1//3*a^-1*cos(x)^3)
@test integrate((a+-1*a*cos(x)^2)^-1*sin(x)^4, x) == :((1/2)*x*a^-1+-1//2*a^-1*cos(x)*sin(x))
@test integrate((a+-1*a*cos(x)^2)^-1*sin(x)^3, x) == :(-1*a^-1*cos(x))
@test integrate((a+-1*a*cos(x)^2)^-1*sin(x)^2, x) == :(x*a^-1)
@test integrate((a+-1*a*cos(x)^2)^-1*sin(x), x) == :(-1*a^-1*arctanh(cos(x)))
@test integrate((a+-1*a*cos(x)^2)^-1*csc(x), x) == :(-1//2*a^-1*arctanh(cos(x))+-1//2*a^-1*cot(x)*csc(x))
@test integrate((a+-1*a*cos(x)^2)^-1*csc(x)^2, x) == :(-1*a^-1*cot(x)+-1//3*a^-1*cot(x)^3)
@test integrate((a+-1*a*cos(x)^2)^-1*csc(x)^3, x) == :(-3//8*a^-1*arctanh(cos(x))+-3//8*a^-1*cot(x)*csc(x)+-1//4*a^-1*csc(x)^3*cot(x))
@test integrate((a+b*cos(x)^2)^-1*sin(x)^7, x) == :(1//5*b^-1*cos(x)^5+b^-3*(a^2+3*b^2+3*a*b)*cos(x)+-1//3*b^-2*cos(x)^3*(a+3b)+-1*a^-1//2*b^-7//2*(a+b)^3*arctan(a^-1//2*b^(1/2)*cos(x)))
@test integrate((a+b*cos(x)^2)^-1*sin(x)^5, x) == :(-1//3*b^-1*cos(x)^3+b^-2*(a+2b)*cos(x)+-1*a^-1//2*b^-5//2*(a+b)^2*arctan(a^-1//2*b^(1/2)*cos(x)))
@test integrate((a+b*cos(x)^2)^-1*sin(x)^3, x) == :(b^-1*cos(x)+-1*a^-1//2*b^-3//2*(a+b)*arctan(a^-1//2*b^(1/2)*cos(x)))
@test integrate((a+b*cos(x)^2)^-1*sin(x), x) == :(-1*a^-1//2*b^-1//2*arctan(a^-1//2*b^(1/2)*cos(x)))
@test integrate((a+b*cos(x)^2)^-1*csc(x), x) == :(-1*(a+b)^-1*arctanh(cos(x))+-1*a^-1//2*b^(1/2)*(a+b)^-1*arctan(a^-1//2*b^(1/2)*cos(x)))
@test integrate((a+b*cos(x)^2)^-1*csc(x)^3, x) == :(-1*(2a+2b)^-1*cot(x)*csc(x)+-1//2*(a+b)^-2*(a+3b)*arctanh(cos(x))+-1*a^-1//2*b^3//2*(a+b)^-2*arctan(a^-1//2*b^(1/2)*cos(x)))
@test integrate((a+b*cos(x)^2)^-1*csc(x)^5, x) == :(-1*(4a+4b)^-1*csc(x)^3*cot(x)+-1//8*(a+b)^-3*(3*a^2+15*b^2+10*a*b)*arctanh(cos(x))+-1*a^-1//2*b^5//2*(a+b)^-3*arctan(a^-1//2*b^(1/2)*cos(x))+-1//8*(a+b)^-2*(3a+7b)*cot(x)*csc(x))
@test integrate((a+b*cos(x)^2)^-1*sin(x)^6, x) == :(-1//8*x*b^-3*(8*a^2+15*b^2+20*a*b)+1//4*b^-1*sin(x)^3*cos(x)+-1*a^-1//2*b^-3*(a+b)^5//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x))+1//8*b^-2*(4a+7b)*cos(x)*sin(x))
@test integrate((a+b*cos(x)^2)^-1*sin(x)^4, x) == :((1/2)*b^-1*cos(x)*sin(x)+-1//2*x*b^-2*(2a+3b)+-1*a^-1//2*b^-2*(a+b)^3//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*sin(x)^2, x) == :(-1*x*b^-1+-1*a^-1//2*b^-1*(a+b)^(1/2)*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1, x) == :(-1*a^-1//2*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*csc(x)^2, x) == :(-1*(a+b)^-1*cot(x)+-1*b*a^-1//2*(a+b)^-3//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*csc(x)^4, x) == :(-1*(3a+3b)^-1*cot(x)^3+-1*(a+b)^-2*(a+2b)*cot(x)+-1*a^-1//2*b^2*(a+b)^-5//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*csc(x)^6, x) == :(-1*(5a+5b)^-1*cot(x)^5+-1*(a+b)^-3*(a^2+3*b^2+3*a*b)*cot(x)+-1//3*(a+b)^-2*cot(x)^3*(2a+3b)+-1*a^-1//2*b^3*(a+b)^-7//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((4+-3*cos(x)^3)^-1*sin(x), x) == :(-1//72*6^2//3*log(2*2^1//3+3^2//3*cos(x)^2+2^2//3*3^1//3*cos(x))+1//36*6^2//3*log(2^2//3+-1*3^1//3*cos(x))+-1//12*2^2//3*3^1//6*arctan(1//3*3^(1/2)*(1+6^1//3*cos(x))))
@test integrate((1+-1*cos(x)^2)^-1, x) == :(-1*cot(x))
@test integrate((1+-1*cos(x)^2)^-2, x) == :(-1*cot(x)+-1//3*cot(x)^3)
@test integrate((1+-1*cos(x)^2)^-3, x) == :(-1*cot(x)+-2//3*cot(x)^3+-1//5*cot(x)^5)
@test integrate((a+b*cos(x)^2)^-1*cos(x)^7, x) == :(1//5*b^-1*sin(x)^5+b^-3*(a^2+b^2+-1*a*b)*sin(x)+1//3*b^-2*sin(x)^3*(a+-2b)+-1*a^3*b^-7//2*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*sin(x)))
@test integrate((a+b*cos(x)^2)^-1*cos(x)^5, x) == :(-1//3*b^-1*sin(x)^3+-1*b^-2*(a+-1b)*sin(x)+a^2*b^-5//2*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*sin(x)))
@test integrate((a+b*cos(x)^2)^-1*cos(x)^3, x) == :(b^-1*sin(x)+-1*a*b^-3//2*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*sin(x)))
@test integrate((a+b*cos(x)^2)^-1*cos(x), x) == :(b^-1//2*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*sin(x)))
@test integrate((a+b*cos(x)^2)^-1*sec(x), x) == :(a^-1*arctanh(sin(x))+-1*a^-1*b^(1/2)*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*sin(x)))
@test integrate((a+b*cos(x)^2)^-1*sec(x)^3, x) == :((1/2)*a^-1*sec(x)*tan(x)+(1/2)*a^-2*(a+-2b)*arctanh(sin(x))+a^-2*b^3//2*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*sin(x)))
@test integrate((a+b*cos(x)^2)^-1*sec(x)^5, x) == :(1//4*a^-1*sec(x)^3*tan(x)+1//8*a^-3*(3*a^2+8*b^2+-4*a*b)*arctanh(sin(x))+-1*a^-3*b^5//2*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*sin(x))+1//8*a^-2*(-4b+3a)*sec(x)*tan(x))
@test integrate((a+b*cos(x)^2)^-1*cos(x)^6, x) == :(1//4*b^-1*cos(x)^3*sin(x)+1//8*x*b^-3*(3*b^2+8*a^2+-4*a*b)+a^5//2*b^-3*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x))+-1//8*b^-2*(-3b+4a)*cos(x)*sin(x))
@test integrate((a+b*cos(x)^2)^-1*cos(x)^4, x) == :((1/2)*b^-1*cos(x)*sin(x)+-1//2*x*b^-2*(-1b+2a)+-1*a^3//2*b^-2*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*cos(x)^2, x) == :(x*b^-1+a^(1/2)*b^-1*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1, x) == :(-1*a^-1//2*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*sec(x)^2, x) == :(a^-1*tan(x)+b*a^-3//2*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*sec(x)^4, x) == :(1//3*a^-1*tan(x)^3+a^-2*(a+-1b)*tan(x)+-1*a^-5//2*b^2*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-1*sec(x)^6, x) == :(1//5*a^-1*tan(x)^5+a^-3*(a^2+b^2+-1*a*b)*tan(x)+1//3*a^-2*tan(x)^3*(-1b+2a)+a^-7//2*b^3*(a+b)^-1//2*arctan(a^-1//2*(a+b)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^2)^-2, x) == :(-1//2*a^-3//2*(a+b)^-3//2*(b+2a)*arctan(a^-1//2*(a+b)^(1/2)*cot(x))+-1//2*b*a^-1*(a+b)^-1*(a+b*cos(x)^2)^-1*cos(x)*sin(x))
@test integrate((a+b*cos(x)^2)^-3, x) == :(-1//8*a^-5//2*(a+b)^-5//2*(3*b^2+8*a^2+8*a*b)*arctan(a^-1//2*(a+b)^(1/2)*cot(x))+-1//4*b*a^-1*(a+b)^-1*(a+b*cos(x)^2)^-2*cos(x)*sin(x)+-3//8*b*a^-2*(a+b)^-2*(a+b*cos(x)^2)^-1*(b+2a)*cos(x)*sin(x))
@test integrate((1+cos(x)^2)^-1, x) == :((1/2)*x*2^(1/2)+-1//2*2^(1/2)*arctan((1+2^(1/2)+cos(x)^2)^-1*cos(x)*sin(x)))
@test integrate((1+cos(x)^2)^-2, x) == :(-3//8*2^(1/2)*arctan((1+2^(1/2)+cos(x)^2)^-1*cos(x)*sin(x))+3//8*x*2^(1/2)+-1*(4+4*cos(x)^2)^-1*cos(x)*sin(x))
@test integrate((1+cos(x)^2)^-3, x) == :(-19//64*2^(1/2)*arctan((1+2^(1/2)+cos(x)^2)^-1*cos(x)*sin(x))+19//64*x*2^(1/2)+-9*(32+32*cos(x)^2)^-1*cos(x)*sin(x)+-1//8*(1+cos(x)^2)^-2*cos(x)*sin(x))
@test integrate((1+-1*cos(x)^2)^(1/2), x) == :(-1*(sin(x)^2)^(1/2)*cot(x))
@test integrate((-1+cos(x)^2)^(1/2), x) == :(-1*(-1*sin(x)^2)^(1/2)*cot(x))
@test integrate((1+-1*cos(x)^2)^3//2, x) == :(-2//3*(sin(x)^2)^(1/2)*cot(x)+-1//3*(sin(x)^2)^3//2*cot(x))
@test integrate((-1+cos(x)^2)^3//2, x) == :(-1//3*(-1*sin(x)^2)^3//2*cot(x)+2//3*(-1*sin(x)^2)^(1/2)*cot(x))
@test integrate((1+-1*cos(x)^2)^-1//2, x) == :(-1*(sin(x)^2)^-1//2*arctanh(cos(x))*sin(x))
@test integrate((-1+cos(x)^2)^-1//2, x) == :(-1*(-1*sin(x)^2)^-1//2*arctanh(cos(x))*sin(x))
@test integrate((1+-1*cos(x)^2)^-3//2, x) == :(-1//2*(sin(x)^2)^-1//2*cot(x)+-1//2*(sin(x)^2)^-1//2*arctanh(cos(x))*sin(x))
@test integrate((-1+cos(x)^2)^-3//2, x) == :((1/2)*(-1*sin(x)^2)^-1//2*cot(x)+(1/2)*(-1*sin(x)^2)^-1//2*arctanh(cos(x))*sin(x))
@test integrate((1+cos(x)^2)^(1/2), x) == :(Elliptic.E(x+(1/2)*pi,-1))
@test integrate((-1+-1*cos(x)^2)^(1/2), x) == :((1+cos(x)^2)^-1//2*(-1+-1*cos(x)^2)^(1/2)*Elliptic.E(x+(1/2)*pi,-1))
@test integrate((a+b*cos(x)^2)^(1/2), x) == :((1+b*a^-1*cos(x)^2)^-1//2*(a+b*cos(x)^2)^(1/2)*Elliptic.E(x+(1/2)*pi,-1*b*a^-1))
@test integrate((1+cos(x)^2)^3//2, x) == :(2*Elliptic.E(x+(1/2)*pi,-1)+-2//3*Elliptic.F(x+(1/2)*pi,-1)+1//3*(1+cos(x)^2)^(1/2)*cos(x)*sin(x))
@test integrate((-1+-1*cos(x)^2)^3//2, x) == :(-2*(1+cos(x)^2)^-1//2*(-1+-1*cos(x)^2)^(1/2)*Elliptic.E(x+(1/2)*pi,-1)+-2//3*(1+cos(x)^2)^(1/2)*(-1+-1*cos(x)^2)^-1//2*Elliptic.F(x+(1/2)*pi,-1)+-1//3*(-1+-1*cos(x)^2)^(1/2)*cos(x)*sin(x))
@test integrate((a+b*cos(x)^2)^3//2, x) == :(1//3*b*(a+b*cos(x)^2)^(1/2)*cos(x)*sin(x)+1//3*(1+b*a^-1*cos(x)^2)^-1//2*(a+b*cos(x)^2)^(1/2)*(2b+4a)*Elliptic.E(x+(1/2)*pi,-1*b*a^-1)+-1//3*a*(1+b*a^-1*cos(x)^2)^(1/2)*(a+b*cos(x)^2)^-1//2*(a+b)*Elliptic.F(x+(1/2)*pi,-1*b*a^-1))
@test integrate((1+cos(x)^2)^-1//2, x) == :(Elliptic.F(x+(1/2)*pi,-1))
@test integrate((-1+-1*cos(x)^2)^-1//2, x) == :((1+cos(x)^2)^(1/2)*(-1+-1*cos(x)^2)^-1//2*Elliptic.F(x+(1/2)*pi,-1))
@test integrate((a+b*cos(x)^2)^-1//2, x) == :((1+b*a^-1*cos(x)^2)^(1/2)*(a+b*cos(x)^2)^-1//2*Elliptic.F(x+(1/2)*pi,-1*b*a^-1))
@test integrate((1+cos(x)^2)^-3//2, x) == :((1/2)*Elliptic.E(x+(1/2)*pi,-1)+-1//2*(1+cos(x)^2)^-1//2*cos(x)*sin(x))
@test integrate((-1+-1*cos(x)^2)^-3//2, x) == :((1/2)*(1+cos(x)^2)^-1//2*(-1+-1*cos(x)^2)^(1/2)*Elliptic.E(x+(1/2)*pi,-1)+(1/2)*(-1+-1*cos(x)^2)^-1//2*cos(x)*sin(x))
@test integrate((a+b*cos(x)^2)^-3//2, x) == :(a^-1*(1+b*a^-1*cos(x)^2)^-1//2*(a+b)^-1*(a+b*cos(x)^2)^(1/2)*Elliptic.E(x+(1/2)*pi,-1*b*a^-1)+-1*b*a^-1*(a+b)^-1*(a+b*cos(x)^2)^-1//2*cos(x)*sin(x))
@test integrate((1+cos(x)^2)^-1//2*cos(x), x) == :(arcsin((1/2)*2^(1/2)*sin(x)))
@test integrate((3+cos(5+3x)^2)^-1//2*cos(5+3x), x) == :(1//3*arcsin((1/2)*sin(5+3x)))
@test integrate((4+-1*cos(x)^2)^-1//2*cos(x), x) == :(arcsinh(1//3*3^(1/2)*sin(x)))
@test integrate((a+b*cos(x)^4)^-1, x) == :(-1//4*2^(1/2)*a^-3//4*(a+b)^-1//4*(a+b+a^(1/2)*(a+b)^(1/2))^-1//2*(a^(1/2)+(a+b)^(1/2))*arctan(a^-1//4*(a+b+a^(1/2)*(a+b)^(1/2))^-1//2*(a^1//4*(a+b+-1*a^(1/2)*(a+b)^(1/2))^(1/2)+2^(1/2)*(a+b)^3//4*cot(x)))+-1//8*2^(1/2)*a^-3//4*(a+b)^-1//4*(a+b+-1*a^(1/2)*(a+b)^(1/2))^-1//2*(a^(1/2)+-1*(a+b)^(1/2))*log(a^(1/2)*(a+b)^1//4+(a+b)^3//4*cot(x)^2+-1*2^(1/2)*a^1//4*(a+b+-1*a^(1/2)*(a+b)^(1/2))^(1/2)*cot(x))+1//4*2^(1/2)*a^-3//4*(a+b)^-1//4*(a+b+a^(1/2)*(a+b)^(1/2))^-1//2*(a^(1/2)+(a+b)^(1/2))*arctan(a^-1//4*(a+b+a^(1/2)*(a+b)^(1/2))^-1//2*(a^1//4*(a+b+-1*a^(1/2)*(a+b)^(1/2))^(1/2)+-1*2^(1/2)*(a+b)^3//4*cot(x)))+1//8*2^(1/2)*a^-3//4*(a+b)^-1//4*(a+b+-1*a^(1/2)*(a+b)^(1/2))^-1//2*(a^(1/2)+-1*(a+b)^(1/2))*log(a^(1/2)*(a+b)^1//4+(a+b)^3//4*cot(x)^2+2^(1/2)*a^1//4*(a+b+-1*a^(1/2)*(a+b)^(1/2))^(1/2)*cot(x)))
@test integrate((a+-1*b*cos(x)^4)^-1, x) == :(-1//2*a^-3//4*(a^(1/2)+b^(1/2))^-1//2*arctan(a^-1//4*(a^(1/2)+b^(1/2))^(1/2)*cot(x))+-1//2*a^-3//4*(a^(1/2)+-1*b^(1/2))^-1//2*arctan(a^-1//4*(a^(1/2)+-1*b^(1/2))^(1/2)*cot(x)))
@test integrate((1+cos(x)^4)^-1, x) == :((1/2)*x*(-1+2^(1/2))^-1//2+-1//8*(-1+2^(1/2))^(1/2)*log(1+2^(1/2)*cot(x)^2+(-2+2*2^(1/2))^(1/2)*cot(x))+1//4*(-1+2^(1/2))^-1//2*arctan((2+(1+2^(1/2))^(1/2)+sin(x)^2*(-2+2^(1/2))+-2*(-1+2^(1/2))^(1/2)*cos(x)*sin(x))^-1*((-1+2^(1/2))^(1/2)*(-1+2*sin(x)^2)+(-2+2^(1/2))*cos(x)*sin(x)))+1//4*(-1+2^(1/2))^-1//2*arctan((2+(1+2^(1/2))^(1/2)+sin(x)^2*(-2+2^(1/2))+2*(-1+2^(1/2))^(1/2)*cos(x)*sin(x))^-1*((-1+2^(1/2))^(1/2)*(1+-2*sin(x)^2)+(-2+2^(1/2))*cos(x)*sin(x)))+1//8*(-1+2^(1/2))^(1/2)*log(2^(1/2)+2*cot(x)^2+-2*(-1+2^(1/2))^(1/2)*cot(x)))
@test integrate((1+-1*cos(x)^4)^-1, x) == :(-1//2*cot(x)+-1//4*2^(1/2)*arctan((1+2^(1/2)+cos(x)^2)^-1*cos(x)*sin(x))+1//4*x*2^(1/2))
@test integrate((a+b*cos(x)^5)^-1, x) == :(2//5*a^-4//5*(a^1//5+b^1//5)^-1//2*(a^1//5+-1*b^1//5)^-1//2*arctan((a^1//5+b^1//5)^-1//2*(a^1//5+-1*b^1//5)^(1/2)*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^1//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^1//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^1//5*b^1//5)^(1/2)*(a^1//5+-1*(-1)^1//5*b^1//5)^-1//2*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^2//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^2//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^2//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^2//5*b^1//5)^(1/2)*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^3//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^3//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^3//5*b^1//5)^(1/2)*(a^1//5+-1*(-1)^3//5*b^1//5)^-1//2*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^4//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^4//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^4//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^4//5*b^1//5)^(1/2)*tan((1/2)*x)))
@test integrate((a+b*cos(x)^6)^-1, x) == :(-1//3*a^-5//6*(a^1//3+b^1//3)^-1//2*arctan(a^-1//6*(a^1//3+b^1//3)^(1/2)*cot(x))+-1//3*a^-5//6*(a^1//3+(-1)^2//3*b^1//3)^-1//2*arctan(a^-1//6*(a^1//3+(-1)^2//3*b^1//3)^(1/2)*cot(x))+-1//3*a^-5//6*(a^1//3+-1*(-1)^1//3*b^1//3)^-1//2*arctan(a^-1//6*(a^1//3+-1*(-1)^1//3*b^1//3)^(1/2)*cot(x)))
@test integrate((a+b*cos(x)^8)^-1, x) == :(1//4*(-1a)^-7//8*(b^1//4+(-1a)^1//4)^-1//2*arctan((-1a)^-1//8*(b^1//4+(-1a)^1//4)^(1/2)*cot(x))+1//4*(-1a)^-7//8*((-1a)^1//4+-1*b^1//4)^-1//2*arctan((-1a)^-1//8*((-1a)^1//4+-1*b^1//4)^(1/2)*cot(x))+1//4*(-1a)^-7//8*((-1a)^1//4+im*b^1//4)^-1//2*arctan((-1a)^-1//8*((-1a)^1//4+im*b^1//4)^(1/2)*cot(x))+1//4*(-1a)^-7//8*((-1a)^1//4+-1*im*b^1//4)^-1//2*arctan((-1a)^-1//8*((-1a)^1//4+-1*im*b^1//4)^(1/2)*cot(x)))
@test integrate((a+-1*b*cos(x)^5)^-1, x) == :(2//5*a^-4//5*(a^1//5+b^1//5)^-1//2*(a^1//5+-1*b^1//5)^-1//2*arctan((a^1//5+b^1//5)^(1/2)*(a^1//5+-1*b^1//5)^-1//2*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^1//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^1//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^1//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^1//5*b^1//5)^(1/2)*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^2//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^2//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^2//5*b^1//5)^(1/2)*(a^1//5+-1*(-1)^2//5*b^1//5)^-1//2*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^3//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^3//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^3//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^3//5*b^1//5)^(1/2)*tan((1/2)*x))+2//5*a^-4//5*(a^1//5+(-1)^4//5*b^1//5)^-1//2*(a^1//5+-1*(-1)^4//5*b^1//5)^-1//2*arctan((a^1//5+(-1)^4//5*b^1//5)^(1/2)*(a^1//5+-1*(-1)^4//5*b^1//5)^-1//2*tan((1/2)*x)))
@test integrate((a+-1*b*cos(x)^6)^-1, x) == :(-1//3*a^-5//6*(a^1//3+-1*b^1//3)^-1//2*arctan(a^-1//6*(a^1//3+-1*b^1//3)^(1/2)*cot(x))+-1//3*a^-5//6*(a^1//3+(-1)^1//3*b^1//3)^-1//2*arctan(a^-1//6*(a^1//3+(-1)^1//3*b^1//3)^(1/2)*cot(x))+-1//3*a^-5//6*(a^1//3+-1*(-1)^2//3*b^1//3)^-1//2*arctan(a^-1//6*(a^1//3+-1*(-1)^2//3*b^1//3)^(1/2)*cot(x)))
@test integrate((a+-1*b*cos(x)^8)^-1, x) == :(-1//4*a^-7//8*(a^1//4+b^1//4)^-1//2*arctan(a^-1//8*(a^1//4+b^1//4)^(1/2)*cot(x))+-1//4*a^-7//8*(a^1//4+-1*b^1//4)^-1//2*arctan(a^-1//8*(a^1//4+-1*b^1//4)^(1/2)*cot(x))+-1//4*a^-7//8*(a^1//4+im*b^1//4)^-1//2*arctan(a^-1//8*(a^1//4+im*b^1//4)^(1/2)*cot(x))+-1//4*a^-7//8*(a^1//4+-1*im*b^1//4)^-1//2*arctan(a^-1//8*(a^1//4+-1*im*b^1//4)^(1/2)*cot(x)))
@test integrate((1+cos(x)^5)^-1, x) == :((5+5*cos(x))^-1*sin(x)+-2//5*(-1+(-1)^2//5)^-1//2*arctanh((-1*(1+(-1)^1//5)^-1*(1+-1*(-1)^1//5))^-1//2*tan((1/2)*x))+2//5*(1+(-1)^3//5)^-1//2*arctan(((1+(-1)^4//5)^-1*(1+-1*(-1)^4//5))^(1/2)*tan((1/2)*x))+2//5*(1+-1*(-1)^4//5)^-1//2*arctan(((1+(-1)^2//5)^-1*(1+-1*(-1)^2//5))^(1/2)*tan((1/2)*x))+-2*(-1*(1+-1*(-1)^3//5)^-1*(1+(-1)^3//5))^(1/2)*(5+5*(-1)^3//5)^-1*arctanh((-1*(1+-1*(-1)^3//5)^-1*(1+(-1)^3//5))^(1/2)*tan((1/2)*x)))
@test integrate((1+cos(x)^8)^-1, x) == :(-1//4*(1+(-1)^1//4)^-1//2*arctan((1+(-1)^1//4)^(1/2)*cot(x))+-1//4*(1+(-1)^3//4)^-1//2*arctan((1+(-1)^3//4)^(1/2)*cot(x))+-1//4*(1+-1*(-1)^1//4)^-1//2*arctan((1+-1*(-1)^1//4)^(1/2)*cot(x))+-1//4*(1+-1*(-1)^3//4)^-1//2*arctan((1+-1*(-1)^3//4)^(1/2)*cot(x)))
@test integrate((1+-1*cos(x)^5)^-1, x) == :(-1*(5+-5*cos(x))^-1*sin(x)+-2//5*(-1+(-1)^4//5)^-1//2*arctanh((-1*(1+(-1)^2//5)^-1*(1+-1*(-1)^2//5))^-1//2*tan((1/2)*x))+2//5*(1+(-1)^1//5)^-1//2*arctan(((1+(-1)^3//5)^-1*(1+-1*(-1)^3//5))^(1/2)*tan((1/2)*x))+2//5*(1+-1*(-1)^2//5)^-1//2*arctan(((1+(-1)^1//5)^-1*(1+-1*(-1)^1//5))^(1/2)*tan((1/2)*x))+2//5*(-1+-1*(-1)^3//5)^-1//2*arctanh((-1*(1+-1*(-1)^4//5)^-1*(1+(-1)^4//5))^(1/2)*tan((1/2)*x)))
@test integrate((1+-1*cos(x)^6)^-1, x) == :(-1//3*cot(x)+-1//3*(1+(-1)^1//3)^-1//2*arctan((1+(-1)^1//3)^(1/2)*cot(x))+-1//3*(1+-1*(-1)^2//3)^-1//2*arctan((1+-1*(-1)^2//3)^(1/2)*cot(x)))
@test integrate((1+-1*cos(x)^8)^-1, x) == :(-1//4*cot(x)+-1//4*(1+im)^-1//2*arctan((1+im)^(1/2)*cot(x))+-1//4*(1+-1im)^-1//2*arctan((1+-1im)^(1/2)*cot(x))+-1//8*2^(1/2)*arctan((1+2^(1/2)+cos(x)^2)^-1*cos(x)*sin(x))+1//8*x*2^(1/2))
@test integrate((1+cos(x)^2)^-1*tan(x), x) == :((1/2)*log(1+cos(x)^2)+-1*log(cos(x)))
@test integrate((a+b*cos(x)^2)^(1/2)*tan(x), x) == :(-1*(a+b*cos(x)^2)^(1/2)+a^(1/2)*arctanh(a^-1//2*(a+b*cos(x)^2)^(1/2)))
@test integrate((1+-1*cos(x)^2)^(1/2)*tan(x), x) == :(-1*(sin(x)^2)^(1/2)+arctanh((sin(x)^2)^(1/2)))
@test integrate((a+b*cos(x)^2)^-1//2*tan(x), x) == :(a^-1//2*arctanh(a^-1//2*(a+b*cos(x)^2)^(1/2)))
@test integrate((1+cos(x)^2)^-1//2*tan(x), x) == :(arctanh((1+cos(x)^2)^(1/2)))
@test integrate((1+-1*cos(x)^2)^-1//2*tan(x), x) == :(arctanh((sin(x)^2)^(1/2)))
@test integrate((a+b*cos(x)^3)^-1*tan(x)^3, x) == :(a^-1*log(cos(x))+(1/2)*a^-1*sec(x)^2+-1//3*a^-1*log(a+b*cos(x)^3)+-1//6*a^-5//3*b^2//3*log(a^2//3+b^2//3*cos(x)^2+-1*a^1//3*b^1//3*cos(x))+1//3*a^-5//3*b^2//3*log(a^1//3+b^1//3*cos(x))+-1//3*3^(1/2)*a^-5//3*b^2//3*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*b^1//3*cos(x))))
@test integrate((a+b*cos(x)^3)^(1/2)*tan(x), x) == :(-2//3*(a+b*cos(x)^3)^(1/2)+2//3*a^(1/2)*arctanh(a^-1//2*(a+b*cos(x)^3)^(1/2)))
@test integrate((a+b*cos(x)^3)^-1//2*tan(x), x) == :(2//3*a^-1//2*arctanh(a^-1//2*(a+b*cos(x)^3)^(1/2)))
@test integrate((a+b*cos(x)^4)^(1/2)*tan(x), x) == :(-1//2*(a+b*cos(x)^4)^(1/2)+(1/2)*a^(1/2)*arctanh(a^-1//2*(a+b*cos(x)^4)^(1/2)))
@test integrate((a+b*cos(x)^4)^-1//2*tan(x), x) == :((1/2)*a^-1//2*arctanh(a^-1//2*(a+b*cos(x)^4)^(1/2)))
@test integrate((a+b*cos(x)^n)^(1/2)*tan(x), x) == :(-2*n^-1*(a+b*cos(x)^n)^(1/2)+2*a^(1/2)*n^-1*arctanh(a^-1//2*(a+b*cos(x)^n)^(1/2)))
@test integrate((a+b*cos(x)^n)^-1//2*tan(x), x) == :(2*a^-1//2*n^-1*arctanh(a^-1//2*(a+b*cos(x)^n)^(1/2)))

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.7 (d trig)^m (a+b (c cos)^n)^p.jl


text:
function _get_trades(url::String)

  res = HTTP.get(url)
  @assert !isnothing(res)

  rawData = JSON.parse(String(res.body))

  df = vcat(DataFrame.(rawData)...)
  df.time = TimeDate.(String.(chop.(df.time)))
  df.price = parse.(Float64, df.price)
  df.size = parse.(Float64, df.size)

  (df, get_pages(res))
end

function trades(ccy::String)
  url = ENDPOINT * "products/" * ccy * "/trades"
  _get_trades(url)
end

function trades_page(ccy::String, pageid::Number, beforeafter::String)
  @assert beforeafter in ["before", "after"]
  url = ENDPOINT * "products/" * ccy * "/trades?($beforeafter)=$(pageid)"
  _get_trades(url)
end

function lasttrades(ccy::String, n::Int64)
  nreqs = div(n, 100) + 1
  res = Array{DataFrame}(undef, nreqs)

  df, pages = trades(ccy)
  if nreqs == 1
    return df[1:n, :]
  else
    nextpage = pages[2]
    for i in eachindex(res)
      df2, pages2 = trades_page(ccy, nextpage, "after")
      nextpage = pages2[2]
      res[i] = df2
    end
  end
  allres = vcat(res...)
  allres[1:n, :]
end


text:
mutable struct SlopePerRanEf
    grouping_vars::Dict{String,Vector{String}}
    # inner constructors
    SlopePerRanEf() = new(Dict{String,Vector{String}}())
    SlopePerRanEf(d::Dict{String,Vector{String}}) = new(d)
end
isequal(x::SlopePerRanEf, y::SlopePerRanEf) = x.grouping_vars == y.grouping_vars
==(x::SlopePerRanEf, y::SlopePerRanEf) = x.grouping_vars == y.grouping_vars
length(x::SlopePerRanEf) = length(values(x.grouping_vars))

"""
    data_response(formula::FormulaTerm, data)

Constructs the response y vector.

Returns a `Vector` of the response variable in the `formula` and present inside `data`.
# Arguments
- `formula`: a `FormulaTerm` created by `@formula` macro.
- `data`:  a `data` object that satisfies the
[Tables.jl](https://github.com/JuliaData/Tables.jl) interface such as a DataFrame.
"""
function data_response(formula::FormulaTerm, data::D) where {D}
    y = response(formula, data)
    return y
end

"""
    data_fixed_effects(formula::FormulaTerm, data)

Constructs the matrix X of fixed-effects (a.k.a. population-level) predictors.

Returns a `Matrix` of the fixed-effects predictors variables in the `formula`
and present inside `data`.

# Arguments
- `formula`: a `FormulaTerm` created by `@formula` macro.
- `data`:  a `data` object that satisfies the
[Tables.jl](https://github.com/JuliaData/Tables.jl) interface such as a DataFrame.
"""
function data_fixed_effects(formula::FormulaTerm, data::D) where {D}
    if has_ranef(formula)
        X = MixedModels.modelmatrix(MixedModel(formula, data))
        X = X[:, 2:end]
    else
        X = StatsModels.modelmatrix(formula, data)
        if hasintercept(formula)
            X = X[:, 2:end]
        end
    end
    return X
end

"""
    data_random_effects(formula::FormulaTerm, data)

Constructs the vector(s)/matrix(ces) Z(s) of random-effects (a.k.a. group-level)
slope predictors.

Returns a `Dict{String, AbstractArray}` of `Vector`/`Matrix` as values of the random-effects
predictors slope variables (keys) in the `formula` and present inside `data`.

# Arguments
- `formula`: a `FormulaTerm` created by `@formula` macro.
- `data`:  a `data` object that satisfies the
[Tables.jl](https://github.com/JuliaData/Tables.jl) interface such as a DataFrame.
"""
function data_random_effects(formula::FormulaTerm, data::D) where {D}
    if !has_ranef(formula)
        return nothing
    end
    slopes = slope_per_ranef(ranef(formula))

    Z = Dict{String,AbstractArray}() # empty Dict
    if length(slopes) > 0
        # add the slopes to Z
        # this would need to create a vector from the column of the X matrix from the
        # slope term
        for slope in values(slopes.grouping_vars)
            if slope isa String
                Z["slope_" * slope] = get_var(term(slope), data)
            else
                for s in slope
                    Z["slope_" * s] = get_var(term(s), data)
                end
            end
        end
    else
        Z = nothing
    end
    return Z
end

"""
    has_ranef(formula::FormulaTerm)

Returns `true` if any of the terms in `formula` is a `FunctionTerm` or false
otherwise.
"""
function has_ranef(formula::FormulaTerm)
    return any(t -> t isa FunctionTerm{typeof(|)}, formula.rhs)
end

"""
    ranef(formula::FormulaTerm)

Returns a tuple of the `FunctionTerm`s parsed as `RandomEffectsTerm`s in `formula`.
If there are no `FunctionTerm`s in `formula` returns `nothing`.
"""
function ranef(formula::FormulaTerm)
    if has_ranef(formula)
        terms = filter(t -> t isa FunctionTerm{typeof(|)}, formula.rhs)
        terms = map(terms) do t
            lhs, rhs = first(t.args_parsed), last(t.args_parsed)
            RandomEffectsTerm(lhs, rhs)
        end
        return terms
    else
        return nothing
    end
end

"""
    n_ranef(formula::FormulaTerm)

Returns the number of `RandomEffectsTerm`s in `formula`.
"""
function n_ranef(formula::FormulaTerm)
    terms = ranef(formula)
    if terms === nothing
        return 0
    elseif any(t -> t isa RandomEffectsTerm, terms)
        counter = 0
        for ts in terms
            if ts.lhs isa Tuple
                counter += tuple_length(ts.lhs)
            else
                counter += 1
            end
        end
        return counter
    else
        # fallback
        return 1
    end
end

"""
    intercept_per_ranef(terms::Tuple{RandomEffectsTerm})

Returns a vector of `String`s where the entries are the grouping variables that have
a group-level intercept.
"""
function intercept_per_ranef(terms::Tuple)
    vec_intercepts = Vector{String}()
    for ts in terms
        if ts.lhs isa ConstantTerm
            push!(vec_intercepts, string(ts.rhs))
        elseif ts.lhs isa Tuple && any(t -> t isa ConstantTerm, ts.lhs)
            push!(vec_intercepts, string(ts.rhs))
        end
    end
    return vec_intercepts
end

"""
    slope_per_ranef(terms::Tuple{RandomEffectsTerm})

Returns a `SlopePerRanEf` object where the entries are the grouping variables that have
a group-level slope.
"""
function slope_per_ranef(terms::Tuple)
    slopes = SlopePerRanEf()
    for ts in terms
        if ts.lhs isa Term
            slopes.grouping_vars[string(ts.rhs)] = [string(ts.rhs)]
        elseif ts.lhs isa Tuple && any(t -> t isa Term, ts.lhs)
            # empty first
            slopes.grouping_vars[string(ts.rhs)] = Vector{String}()
            # now populate
            for tleft in ts.lhs
                if tleft isa Term
                    push!(slopes.grouping_vars[string(ts.rhs)], string(tleft))
                end
            end
        end
    end
    return slopes
end

"""
    get_idx(term::Term, data)

Returns a tuple with the first element as the ID vector of `Int`s that represent
group membership for a specific random-effect intercept group `t` of observations
present in `data`. The second element of the tuple is a `Dict` specifying which string is
which integer in the ID vector.
"""
function get_idx(t::Term, data::D) where {D}
    col = Symbol(t)
    idx = Tables.getcolumn(data, col)
    return convert_str_to_indices(idx)
end

"""
    get_var(term::Term, data)

Returns the corresponding vector of column in `data` for the a specific
random-effect slope `term` of observations present in `data`.
"""
function get_var(t::Term, data::D) where {D}
    col = Symbol(t)
    return Tables.getcolumn(data, col)
end


text:
@testset "NumField/NfRel" begin
  include("NfRel/Conjugates.jl")
  include("NfRel/ComplexEmbeddings.jl")
  include("NfRel/Simplify.jl")
end


text:
# https://github.com/JuliaPlots/Makie.jl
# http://juliaplots.org/MakieReferenceImages/
# http://juliaplots.org/MakieReferenceImages/gallery/index.html

# Need to restart Julia in order to switch the backend
glmakie    = false
cairomakie = false
wglmakie   = false

glmakie    = true
#cairomakie = true
#wglmakie = true

# GLMakie.jl
# GLMakie.jl is in a sub-directory of Makie.jl.  The old GLMakie.jl package
# is deprecated.
# add GLMakie
if glmakie
    using GLMakie
end

# CairoMakie.jl
if cairomakie
    using CairoMakie
end

if wglmakie
    using WGLMakie
end

#using AbstractPlotting # not necessary
using LinearAlgebra
import GeometryBasics
gb = GeometryBasics
using Colors

function learn_basic_lines()
    x = range(0, 10, length=100)
    y = sin.(x)
    #fig = lines(x, y)
    lines(x, y)
    # @show typeof(fig)
    # if glmakie
    # # fig isa AbstractPlotting.FigureAxisPlot
    #     display(fig)
    # end
    # if cairomakie
    # # fig isa Makie.FigureAxisPlot
    #     display(fig)
    # end

    # fig
end

function learn_basic_scatter()
    x = range(0, 10, length=100)
    y = sin.(x)
    scatter(x, y)
end

function learn_basic_multiple_plots()
    x = range(0, 10, length=100)
    y1 = sin.(x)
    y2 = cos.(x)

    lines(x, y1)
    lines!(x, y2)
    current_figure()
end

function learn_basic_attributes()
    x = range(0, 10, length=100)
    y1 = sin.(x)
    y2 = cos.(x)

    lines(x, y1, color = colorant"darkred", linewidth=5)
    lines!(x, y2, color = colorant"darkblue")
    current_figure()
end

function learn_basic_attributes_2()
    x = range(0, 10, length=100)
    y1 = sin.(x)
    y2 = cos.(x)

    scatter(x, y1, color = colorant"darkred", markersize=5)
    plt = scatter!(x, y2, color = colorant"darkblue", markersize=10)
    plt.color = colorant"darkgreen"

    current_figure()
end
function learn_basic_attributes_3()
    x = range(0, 10, length=100)
    y1 = sin.(x)
    y2 = cos.(x)

    scatter(x, y1, color = :red, markersize = range(5, 15, length=100))
    sc = scatter!(x, y2, color = range(0, 1, length=100), colormap = :thermal)
    sc.colorrange = (0.25, 0.75)
    current_figure()
end

function learn_basic_attributes_4()
    x = range(0, 10, length=100)
    y = sin.(x)

    colors = repeat([:crimson, :dodgerblue, :slateblue1, :sienna1, :orchid1], 20)

    scatter(x, y, color = colors, markersize = 20)
end

function learn_basic_simple_legend()
    x = range(0, 10, length=100)
    y1 = sin.(x)
    y2 = cos.(x)

    lines(x, y1, color = :red, label = "sin")
    lines!(x, y2, color = :blue, label = "cos")
    axislegend()
    current_figure()
end
"""

Construct `Figure` manually.
"""
function learn_basic_subplots()
    #x = LinRange(0, 10, 100)
    x = range(0, 10, length=100)
    y = sin.(x)

    fig = Figure()
    lines(fig[1, 1], x, y, color = :red)
    lines(fig[1, 2], x, y, color = :blue)
    lines(fig[2, 1:2], x, y, color = :green)
    display(fig)
end

function learn_basic_constructing_axes()
    fig = Figure()
    ax1 = Axis(fig[1,1])
    ax2 = Axis(fig[1,2])
    ax3 = Axis(fig[2, 1:2])
    display(fig)

    # .. from IntervalSets.jl
    lines!(ax1, -pi..pi, sin)
    lines!(ax2, 0..10, cos)
    lines!(ax3, 0..10, sqrt)

    ax1.title = "sin"
    ax2.title = "cos"
    ax3.title = "sqrt"

    ax1.ylabel = "amplitude"
    ax3.ylabel = "amplitude"
    ax3.xlabel = "time"
end

function learn_basic_legend()
    fig = Figure()
    ax1, l1 = lines(fig[1,1], 0..10, sin, color=:darkred)
    ax2, l2 = lines(fig[2,1], 0..10, cos, color=:lightblue)
    Legend(fig[1:2,2], [l1,l2], ["sin", "cos"])

    display(fig)
end

function learn_basic_colorbar()
    fig, ax, hm = heatmap(randn(20,20))

    Colorbar(fig[1,2], hm)

    #@show fig # = (backgroundcolor = :pink)
    #fig.attributes[:backgroundcolor] = :pink
    fig.scene.backgroundcolor = :pink
    ax.aspect = 1
    display(fig)

    # does not work
    # heatmap(
    #     randn(20,20),
    #     figure = (backgroundcolor = :pink)
    # )
end

function learn_plotting_arrows()
    x2y = 2/3
    fig = Figure(resolution = (500, 500/x2y))
    ax = Axis(fig[1,1], backgroundcolor=:black)
    ax.aspect = x2y

    xs = range(0, 2pi, length=20)
    ys = range(0, 2pi/x2y, length=20)

    us = [sin(x)*cos(y) for x in xs, y in ys]
    vs = [-cos(x)*sin(y) for x in xs, y in ys]
    # strength = vec(
    #     @. sqrt(us^2 + vs^2)
    # )
    strength = vec(sqrt.(us .^ 2 .+ vs .^ 2))

    ar = arrows!(xs, ys, us, vs, arrowsize=0.1, lengthscale=0.3, arrowcolor=strength, linecolor=strength)

    display(fig)
end

function learn_plotting_arrows_2()
    ps = [Point3(x,y,z) for x in -5:2:5 for y in -5:2:5 for z in -5:2:5]
    ns = map(p -> 0.1 * Vec3(p[2], p[3], p[1]), ps)
    lengths = norm.(ns)
    fig, ax, ar = arrows(
        ps, ns, fxaa=true, # anti_aliasing
        #linecolor = :gray, arrowcolor = :black,
        linecolor = lengths, arrowcolor = lengths,
        linewidth = 0.1, #arrowsize = Vec3(0.3,0.3,0.4),
        align = :center, axis = (type=Axis3,),
        quality = 50
    )
    #di = DataInspector(fig)
    display(fig)
end

# broken
function sec_animated_surface_wireframe()
    scene = Scene();
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r)/r)
    end

    r = range(-2, stop = 2, length = 50)
    surf_func(i) = [Float32(xy_data(x*i, y*i)) for x = r, y = r]
    z = surf_func(20)
    surf = surface!(scene, r, r, z)[end]

    wf = wireframe!(scene, r, r, lift(x-> x .+ 1.0, surf[3]),
        linewidth = 2f0, color = lift(x-> to_colormap(x)[5], surf[:colormap])
    )
    N = 150
    scene
    record(scene, "output.mp4", range(5, stop = 40, length = N)) do i
        surf[3] = surf_func(i)
    end
end

function sec_arrows_3d()
    function SphericalToCartesian(r::T,θ::T,ϕ::T) where T<:AbstractArray
        x = @.r*sin(θ)*cos(ϕ)
        y = @.r*sin(θ)*sin(ϕ)
        z = @.r*cos(θ)
        Point3f0.(x, y, z)
    end
    n = 100^2 #number of points to generate
    r = ones(n);
    θ = acos.(1 .- 2 .* rand(n))
    φ = 2π * rand(n)
    pts = SphericalToCartesian(r,θ,φ)
    arrows(pts, (normalize.(pts) .* 0.1f0), arrowsize = 0.02, linecolor = :green, arrowcolor = :darkblue)
end

function sec_axis_surface()
    vx = -1:0.01:1
    vy = -1:0.01:1

    f(x, y) = (sin(x*10) + cos(y*10)) / 4
    scene = Scene(resolution = (500, 500))
    # One way to style the axis is to pass a nested dictionary / named tuple to it.
    surface!(scene, vx, vy, f, axis = (frame = (linewidth = 2.0,),))
    psurf = scene[end] # the surface we last plotted to scene
    # One can also directly get the axis object and manipulate it
    axis = scene[Axis] # get axis

    # You can access nested attributes likes this:
    axis[:names, :axisnames] = ("\\bf{ℜ}[u]", "\\bf{𝕴}[u]", " OK\n\\bf{δ}\n γ")
    tstyle = axis[:names] # or just get the nested attributes and work directly with them

    tstyle[:textsize] = 10
    tstyle[:textcolor] = (:red, :green, :black)
    tstyle[:font] = "helvetica"

    psurf[:colormap] = :RdYlBu
    wh = widths(scene)
    t = text!(
        campixel(scene),
        "Multipole Representation of first resonances of U-238",
        position = (wh[1] / 2.0, wh[2] - 20.0),
        align = (:center,  :center),
        textsize = 20,
        font = "helvetica",
        raw = :true
    )
    c = lines!(scene, Circle(Point2f0(0.1, 0.5), 0.1f0), color = :red, offset = Vec3f0(0, 0, 1))
    scene
    #update surface
    # TODO explain and improve the situation here
    psurf.converted[3][] = f.(vx .+ 0.5, (vy .+ 0.5)')
    scene
end
function sec_merged_color_mesh()
    function colormesh((geometry, color))
        mesh1 = gb.normal_mesh(geometry)
        npoints = length(GeometryBasics.coordinates(mesh1))
        return gb.pointmeta(mesh1; color=fill(color, npoints))
    end
    # create an array of differently colored boxes in the direction of the 3 axes
    x = Vec3f0(0); baselen = 0.2f0; dirlen = 1f0
    rectangles = [
        (Rect(Vec3f0(x), Vec3f0(dirlen, baselen, baselen)), RGBAf0(1,0,0,1)),
        (Rect(Vec3f0(x), Vec3f0(baselen, dirlen, baselen)), RGBAf0(0,1,0,1)),
        (Rect(Vec3f0(x), Vec3f0(baselen, baselen, dirlen)), RGBAf0(0,0,1,1))
    ]

    meshes = map(colormesh, rectangles)
    #@show meshes
    mesh(merge(meshes))
end
function sec_shading()
    #mesh(gb.Sphere(Point3f0(0), 1f0), color = :orange, shading = true)
    mesh(gb.Sphere(gb.Point3(0.0), 1.0), color = :orange, shading = true)
end

function sec_volume()
    volume(rand(32, 32, 32), algorithm = :mip)
end

function test_rect()
    rect1 = gb.Rect(0.,0., 0.,  1.,2., 3.)
    mesh(rect1, color = :purple, shading = true)
    rect2 = gb.Rect(1.,1., 1.,  1.,2., 3.)
    mesh!(rect2, color = :orange, shading = true)

    #mesh([rect1,rect2], color = [:purple, :orange], shading = true)
    current_figure()
end

function test_mesh()
    vertices = [
        0.0 0.0 1.0;
        1.0 0.0 1.0;
        1.0 1.0 1.0;
        0.0 1.0 0.0;
    ]
    faces = [
       1 2 3;
     #  3 4 1;
    ]
    scene = mesh(vertices, faces, color = :orange)
    faces = [
       3 4 1;
    ]
    mesh!(vertices, faces, color = :purple)

    # mesh(
    # [(0.0, 0.0, 1.0), (0.5, 1.0, 0.0), (1.0, 0.0, 0.0)], color = :orange,
    # shading = false
    # )

    current_figure()
end

## Basic Tutorial ----------------------------------------------- Basic Tutorial
#learn_basic_lines()
#learn_basic_scatter()
#learn_basic_multiple_plots()

## http://makie.juliaplots.org/dev/basic-tutorial.html#Array-attributes
#learn_basic_attributes()
#learn_basic_attributes_2()
#learn_basic_attributes_3()
#learn_basic_attributes_4()

## http://makie.juliaplots.org/dev/basic-tutorial.html#Simple-legend
#learn_basic_simple_legend()
#learn_basic_subplots()
#learn_basic_constructing_axes()
#learn_basic_legend()
#learn_basic_colorbar()

## Layout Tutorial --------------------------------------------- Layout Tutorial

## Plotting Functions --------------------------------------- Plotting Functions
#learn_plotting_arrows()
learn_plotting_arrows_2()

#sec_animated_surface_wireframe()
#sec_arrows_3d()
#sec_axis_surface()
#sec_merged_color_mesh()
#sec_shading()
#sec_volume()

#test_rect()
#test_mesh()


text:
# This file was generated, do not modify it. # hide
tm = TunedModel(model=dtc, ranges=r, operation=predict_mode,
                measure=misclassification_rate)
m = machine(tm, X, y)
fit!(m)
fitted_params(m).best_model.max_depth


text:
@polyvar x y
@variable(model, t)
p1 = t*(1 + x*y)^2 - x*y + (1 - y)^2
p2 = (1 - x*y)^2 + x*y + t*(1 + y)^2
@constraint(model, p1 >= 0)
@constraint(model, p2 >= 0)
@objective(model, Min, t)
optimize!(model)
println("Solution: $(value(t))")
# Returns t ~ 0.25

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
@testset "$TEST $G" begin

g3 = PathGraph(5, G)
g4 = PathDiGraph(5, DG)

d1 = EdgeMap(g3, float([ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
d2 = EdgeMap(g3, sparse(float([ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])))

@test a_star(g3, 1, 4, d1) ==
    a_star(g4, 1, 4, d1) ==
    a_star(g3, 1, 4, d2)

@test a_star(g4, 4, 1) == Vector{Edge{V}}()

end # testset


text:
using Sigma

PI = 3.1415926535897
# arc of a circle
X = uniform(-2,2)
Y = uniform(-2,2)
Z = uniform(-2,2)
cond1 = X*X + Y*Y == 0.1
# cond2 = in(asin(X/Y), (-1/2PI,1/2PI))
cond2 = asin(X/Y) > 1/2PI
xy = RandArray(Any[X,Y])
samples = rand(xy, cond1 & cond2, 100)

# x = [sample[1] for sample in samples]
# y = [sample[2] for sample in samples]
#
# plot(x=x,y=y,Geom.point)

# Ellipsoid
## =========
a = 2
b = 1
c = 0.5
ellip_cond = ((X*X)/(a*a) + (Y*Y)/(b*b) + (Z*Z)/(c*c)) ==  1.0
xyz = RandArray(Any[X,Y,Z])
samples_ellipsoid = rand(xyz, ellip_cond, 100)

## Fractal
## =======
Zsqr = Z*Z
Ysqr = Y*Y
Xsqr = X*X

xcond = (((3Zsqr - Xsqr - Ysqr)*X*(Xsqr-3Ysqr))/(Xsqr+Ysqr)) == X
ycond = (((3Zsqr - Xsqr - Ysqr) * Y * (3*Xsqr - Ysqr))/(Xsqr+Ysqr)) == Y
zcond = (Z*(Zsqr-3*Xsqr - 3*Ysqr)) == Z
samples_fractal = rand(xyz, xcond & ycond & zcond, 10)


text:
using DataIndexer
using Test

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


text:
"""
    NoiseHandling(ασ = 1.1, callback = s -> s > 0,
                  r = .3, ϵ = 0., θ = .5, c = .3)

The standard settings may work well for noisy objective functions. To avoid
premature convergence due to too fast decrease of sigma, there is the option
`noise_handling = CMAEvolutionStrategy.NoiseHandler(ασ = 1.1, callback = s -> s > 0)`.
Choose `ασ` such that sigma decreases slowly (and does not diverge). The callback
function can be used to change the objective function, e.g. increase the
measurement duration, if this leads to smaller noise. The variable `s` indicates
whether CMA-ES can handle the current level of noise: `s > 0` indicates that the
noise level is too high. Whenever the callback returns `true`, sigma gets multiplied by
`ασ`, which is the case when `s > 0`, with the default callback.

For details on noise handling see [Hansen 2009](http://dx.doi.org/10.1109/TEVC.2008.924423).
"""
Base.@kwdef mutable struct NoiseHandling{F}
    r::Float64 = .3
    ϵ::Float64 = 0.
    θ::Float64 = .5
    c::Float64 = .3
    ασ::Float64 = 1.1
    sbar::Float64 = 0.
    fevals::Int = 0
    callback::F = s -> s > 0
end
"""
    NoiseHandling(n; kwargs...) = NoiseHandling(; ασ = 1 + 2/(n + 10), kwargs...)
"""
NoiseHandling(n; kwargs...) = NoiseHandling(; ασ = 1 + 2/(n + 10), kwargs...)
noisefevals(::Any) = 0
noisefevals(n::NoiseHandling) = n.fevals
noise_handling!(::Any, ::Any, ::Any, ::Any, ::Any) = nothing
mutate(y, ϵ) = ϵ == 0 ? y : error("Not implemented.")
function rank_changes(x, y, λreev)
    λ = length(x)
    rank = sortperm(sortperm([x; y]))
    [rank[i] - rank[λ + i] - sign(rank[i] - rank[λ + i]) for i in 1:λreev], rank
end
Δlim(R, θ, λ) = θ * median([abs(i - R) for i in 1:2λ-1])
function noise_handling!(p::Parameters{<:Any,<:NoiseHandling}, f, y, fvals, perm)
    n = p.noise_handling
    rate = p.λ * n.r
    λreev = floor(Int, rate) + (rand(p.rng) < rate - floor(rate))
    n.fevals += λreev
    yreev = mutate(@view(y[:, 1:λreev]), n.ϵ)
    new_fvals = evaluate(p, f, compute_input(p, yreev))
    Δ, joint_perm = rank_changes(fvals, [new_fvals; fvals[λreev+1:end]], λreev)
    s = mean(2 * abs(Δ[i]) -
             Δlim(joint_perm[i] - (fvals[i] > fvals[λreev + i]), n.θ, p.λ) -
             Δlim(joint_perm[λreev + i] - (fvals[λreev + i] > fvals[i]), n.θ, p.λ)
             for i in 1:λreev)
    n.sbar = (1 - n.c) * n.sbar + n.c * s
    if n.callback(n.sbar)
        p.sigma.σ *= n.ασ
    end
    @views fvals[1:λreev] .= .5 * fvals[1:λreev] .+ .5 * new_fvals
    perm .= sortperm(fvals)
    nothing
end


text:
function init(start_time, m::JuniperProblem; incumbent = nothing)
    hash_val = string(hash(hcat(m.l_var,m.u_var)))
    node = BnBNode(1, 1, m.l_var, m.u_var, m.relaxation_solution, 0, :Branch, m.relaxation_status, m.relaxation_objval,[],hash_val)
    obj_gain_m = zeros(m.num_disc_var)
    obj_gain_p = zeros(m.num_disc_var)
    obj_gain_mc = zeros(Int64, m.num_disc_var)
    obj_gain_pc = zeros(Int64, m.num_disc_var)
    obj_gain_inf = zeros(Int64, m.num_disc_var)
    obj_gain = GainObj(obj_gain_m, obj_gain_p, obj_gain_mc, obj_gain_pc, obj_gain_inf)
    disc2var_idx = zeros(m.num_disc_var)
    var2disc_idx = zeros(m.num_var)
    disc_i = 1
    for i=1:m.num_var
        if m.var_type[i] != :Cont
            disc2var_idx[disc_i] = i
            var2disc_idx[i] = disc_i
            disc_i += 1
        end
    end
    factor = 1
    if m.obj_sense == :Min
        factor = -1
    end
    bnbTree = BnBTreeObj()
    bnbTree.m           = m
    bnbTree.obj_gain    = obj_gain
    bnbTree.disc2var_idx = disc2var_idx
    bnbTree.var2disc_idx = var2disc_idx
    bnbTree.options     = m.options
    bnbTree.obj_fac     = factor
    bnbTree.start_time  = start_time
    bnbTree.nsolutions  = 0
    bnbTree.branch_nodes = [node]
    bnbTree.best_bound  = m.relaxation_objval
    bnbTree.limit = :None
    bnbTree.global_solver = is_global_status(m.status)

    if incumbent !== nothing
        bnbTree.incumbent = incumbent
        bnbTree.nsolutions += 1
        if m.options.incumbent_constr
            add_incumbent_constr(m,bnbTree.incumbent)
            m.ncuts += 1
        end
    end

    return bnbTree
end

function new_default_node(idx, level, l_var, u_var, solution;
                            var_idx=0, state=:Solve, relaxation_state=MOI.OPTIMIZE_NOT_CALLED, best_bound=NaN, path=[])

    l_var = copy(l_var)
    u_var = copy(u_var)
    solution = copy(solution)
    hash_val = string(hash(hcat(l_var,u_var)))
    return BnBNode(idx, level, l_var, u_var, solution, var_idx, state, relaxation_state, best_bound, path, hash_val)
end

function new_left_node(node, u_var;
                       var_idx=0, state=:Solve, relaxation_state=MOI.OPTIMIZE_NOT_CALLED, best_bound=NaN, path=[])
    l_var = copy(node.l_var)
    u_var = copy(u_var)
    solution = NaN*ones(length(node.solution))
    hash_val = string(hash(hcat(l_var,u_var)))
    return BnBNode(node.idx*2, node.level+1, l_var, u_var, solution, var_idx, state, relaxation_state, best_bound, path, hash_val)
end

function new_right_node(node, l_var;
                       var_idx=0, state=:Solve, relaxation_state=MOI.OPTIMIZE_NOT_CALLED, best_bound=NaN, path=[])
    l_var = copy(l_var)
    u_var = copy(node.u_var)
    solution = NaN*ones(length(node.solution))
    hash_val = string(hash(hcat(l_var,u_var)))
    return BnBNode(node.idx*2+1, node.level+1, l_var, u_var, solution, var_idx, state, relaxation_state, best_bound, path, hash_val)
end

function init_time_obj()
    return TimeObj(0.0,0.0,0.0,0.0,0.0)
end

function new_default_step_obj(m,node)
    gains_m = zeros(m.num_disc_var)
    gains_mc = zeros(Int64, m.num_disc_var)
    gains_p = zeros(m.num_disc_var)
    gains_pc = zeros(Int64, m.num_disc_var)
    gains_inf = zeros(Int64, m.num_disc_var)
    gains = GainObj(gains_m, gains_p, gains_mc, gains_pc, gains_inf)
    idx_time = 0.0
    node_idx_time = 0.0
    upd_gains_time = 0.0
    node_branch_time = 0.0
    branch_time = 0.0
    step_obj = StepObj()
    step_obj.node             = node
    step_obj.var_idx          = 0
    step_obj.state            = :None
    step_obj.nrestarts        = 0
    step_obj.gain_gap         = 0.0
    step_obj.obj_gain         = gains
    step_obj.idx_time         = idx_time
    step_obj.node_idx_time    = node_idx_time
    step_obj.upd_gains_time   = upd_gains_time
    step_obj.node_branch_time = node_branch_time
    step_obj.branch_time      = branch_time
    step_obj.integral         = []
    step_obj.branch           = []
    step_obj.counter          = 0
    step_obj.upd_gains        = :None
    step_obj.strong_disc_vars = Int64[]
    step_obj.branch_strategy  = :None
    step_obj.strong_branching = nothing
    return step_obj
end

_ci(vi::MOI.VariableIndex, S::Type) = MOI.ConstraintIndex{MOI.VariableIndex, S}(vi.value)
function init_juniper_problem!(jp::JuniperProblem, model::MOI.AbstractOptimizer)
    jp.num_var = MOI.get(model, MOI.NumberOfVariables())

    jp.status = MOI.OPTIMIZE_NOT_CALLED
    jp.relaxation_status = MOI.OPTIMIZE_NOT_CALLED
    jp.nlp_data = nothing
    jp.objective = nothing
    for attr in MOI.get(model, MOI.ListOfModelAttributesSet())
        if attr isa MOI.NLPBlock
            jp.nlp_data = MOI.get(model, MOI.NLPBlock())
        elseif attr isa MOI.ObjectiveFunction
            jp.objective = MOI.get(model, attr)
        end
    end

    jp.objval = NaN
    jp.best_bound = NaN
    jp.solution = fill(NaN, jp.num_var)
    jp.nsolutions = 0
    jp.solutions = []
    jp.num_disc_var = 0
    jp.nintvars = 0
    jp.nbinvars = 0
    jp.nnodes = 1 # is set to one for the root node
    jp.ncuts = 0
    jp.nbranches = 0
    jp.nlevels = 1
    jp.relaxation_time = 0.0

    jp.start_time = time()

    jp.nl_solver = model.options.nl_solver

    if model.options.mip_solver !== nothing
        jp.mip_solver = model.options.mip_solver
    end
    jp.options = model.options
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
        jp.obj_sense = :Min
    else
        jp.obj_sense = :Max
    end
    vis = MOI.get(model, MOI.ListOfVariableIndices())
    jp.l_var = map(vis) do vi
        lb = MOIU.get_bounds(model, Float64, vi)[1]
        MOI.is_valid(model, _ci(vi, MOI.ZeroOne)) ? max(lb, 0.0) : lb
    end
    jp.u_var = map(vis) do vi
        ub = MOIU.get_bounds(model, Float64, vi)[2]
        MOI.is_valid(model, _ci(vi, MOI.ZeroOne)) ? min(ub, 1.0) : ub
    end
    integer_bool_arr = map(vis) do vi
        MOI.is_valid(model, _ci(vi, MOI.Integer))
    end
    binary_bool_arr = map(vis) do vi
        MOI.is_valid(model, _ci(vi, MOI.ZeroOne))
    end
    primal_start_arr = map(vis) do vi
        start = MOI.get(model, MOI.VariablePrimalStart(), vi)
        start === nothing ? 0.0 : start
    end
    jp.primal_start = primal_start_arr
    jp.nintvars = sum(integer_bool_arr)
    jp.nbinvars = sum(binary_bool_arr)
    jp.num_disc_var = sum(integer_bool_arr)+sum(binary_bool_arr)
    jp.var_type = [:Cont for i in 1:jp.num_var]
    jp.var_type[integer_bool_arr .== true] .= :Int
    jp.var_type[binary_bool_arr .== true] .= :Bin
    jp.disc2var_idx = zeros(jp.num_disc_var)
    jp.var2disc_idx = zeros(jp.num_var)
    int_i = 1
    for i=1:jp.num_var
        if jp.var_type[i] != :Cont
            jp.disc2var_idx[int_i] = i
            jp.var2disc_idx[i] = int_i
            int_i += 1
        end
    end
end

function new_default_strong_branch_step_obj(var_idx)
    strong_step_obj = StrongBranchStep()
    strong_step_obj.var_idx                 = var_idx
    strong_step_obj.l_relaxation_state      = MOI.OPTIMIZE_NOT_CALLED
    strong_step_obj.r_relaxation_state      = MOI.OPTIMIZE_NOT_CALLED
    strong_step_obj.init_restart            = false
    return strong_step_obj
end


text:
# Automatically generated using Clang.jl

const CSFML_VERSION_MAJOR = 2
const CSFML_VERSION_MINOR = 5
const CSFML_VERSION_PATCH = 0

# Skipping MacroDefinition: CSFML_API_EXPORT extern "C" __attribute__ ( ( __visibility__ ( "default" ) ) )
# Skipping MacroDefinition: CSFML_API_IMPORT CSFML_EXTERN_C __attribute__ ( ( __visibility__ ( "default" ) ) )
# Skipping MacroDefinition: CSFML_DEPRECATED __attribute__ ( ( deprecated ) )

const sfFalse = 0
const sfTrue = 1
const sfBool = Cint
const sfInt8 = UInt8
const sfUint8 = Cuchar
const sfInt16 = Int16
const sfUint16 = UInt16
const sfInt32 = Cint
const sfUint32 = UInt32
const sfInt64 = Clonglong
const sfUint64 = Culonglong


text:
# please swith it off if you do not use CUDA
using Fire
#using Yao, Yao.ConstGate, BitBasis
#using QuantumPEPS
using Statistics: var, mean, std
using Printf
#using Flux
#using BenchmarkTools

include("data/decoder.jl")

@main function show_exact(nx=4, ny=4)
    model = J1J2(nx, ny; J2=0.5, periodic=false)
    @show ground_state(model)[1]
end

@main function decode(nx::Int, ny::Int; depth::Int=5, nv::Int=1)
    decode(nx, ny, depth, nv)
end

@main function vargrad(nx::Int=4, ny::Int=4;
                    depth::Int=5, nv::Int=1,
                    nbatch::Int=1024, maxiter::Int=20,
                    use_mps::Bool=false,
                    )
    gradients = readdlm("data/gradients-nx$nx-ny$ny-nv$nv-d$depth-B$nbatch-iter$maxiter$(use_mps ? "mps" : "").dat")
    gradients = vec(gradients)
    println("Std-Gradient = ", @sprintf "%.3f" mean(std(gradients)))
    #println("Mean-Abs-Gradient = ", @sprintf "%.3f" mean(abs.(gradients)))

    gradients = readdlm("data/fixparam-gradients-nx$nx-ny$ny-nv$nv-d$depth-B$nbatch-iter$maxiter$(use_mps ? "mps" : "").dat")
    println("Std-Sampling = ", @sprintf "%.3f" mean(std(gradients, dims=2)))
    #println("Std-Std-Sampling = ",@sprintf "%.3f" std(std(gradients, dims=2)))
end

@main function grad_analyse()
    for use_mps in [true, false]
        #for nbatch in [512, 1024, 2048, 4096]
        for nbatch in [4096]
            for nx in [4, 6]
                println("MPS: $use_mps, Size: $nx, Batch: $nbatch")
                maxiter = nx==4 ? 50 : 20
                depth = use_mps ? (nx==4 ? 3 : 2) : 5
                nv = use_mps ? nx+1 : 1
                @show (nx, nx, nbatch, maxiter, nv, use_mps, depth)
                vargrad(nx, nx, nbatch=nbatch, maxiter=maxiter, nv=nv, use_mps=use_mps, depth=depth)
            end
        end
    end
end


text:
using Test, Pkg
using MPI

if haskey(Pkg.installed(), "CuArrays")
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
prodrank = prod(1:rank+1)

# Not possible to PROD a Char (and neither Int8 with OpenMPI)
for T in setdiff([Base.uniontypes(MPI.MPIDatatype)...], [Char, Int8, UInt8])
    A = ArrayType{T}(fill(T(rank+1), 4))
    B = similar(A)
    MPI.Scan!(A, B, *, comm)
    @test B == ArrayType{T}(fill(T(prodrank), 4))

    B = MPI.Scan(A, *, comm)
    @test B isa ArrayType{T}
    @test B == ArrayType{T}(fill(T(prodrank), 4))

    B = MPI.Scan(T(rank+1), *, comm)
    @test B[1] == prodrank
end

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


text:
#ENV["GUROBI_HOME"] = "/Library/gurobi910/mac64"

using ClutteredEnvPathOpt
using LightGraphs
using Pipe
using Plots
using JuMP, Gurobi

function plot_scatter(num_obstacles, num_samples)
    io = open("points.txt", "w")
    close(io)

    X = []
    Y = []
    N = []

    for i in 1:num_samples
        println("SAMPLE: ", i)

        x, y, n = gen_point(num_obstacles)

        push!(X, x)
        push!(Y, y)
        push!(N, n)
    end

    scatter(X, Y, title="Computed Biclique Cover Sizes vs. Optimal")
    xlabel!("Computed Biclique Cover Sizes")
    ylabel!("Optimal (via ILP)")
    Plots.savefig("scatter")

    io = open("points.txt", "a")
    println(io, "X=", X)
    println(io, "Y=", Y)
    println(io, "N=", N)
    close(io)

    return (X, Y, N)
end

function gen_point(num_obstacles)
    obstacles, points, g, faces = ClutteredEnvPathOpt.plot_new(num_obstacles)
    skeleton = LabeledGraph(g)

    cover = ClutteredEnvPathOpt.find_biclique_cover(skeleton, faces)
    x = length(cover)

    feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(faces))
    E = @pipe complement(feg.graph) |> incidence_matrix(_)

    lower = Int(ceil(log2(length(faces))))
    upper = x

    while upper - lower > 1
        t = (upper + lower) ÷ 2 # integer division
        is_feasible = is_biclique_feasible(E, t)

        if is_feasible
            upper = t
        else
            lower = t
        end
    end

    y = lower
    plot_optimal(E, obstacles, points, lower, upper)

    io = open("points.txt", "a")
    println(io, "x=",x,";y=",y,";")
    close(io)

    return (x, y, length(points))
end

function is_biclique_feasible(E, t)
    println("CHECKING FEASIBILITY: t = ", t)

    n = size(E)[1]
    J = 1:n

    model = JuMP.Model(Gurobi.Optimizer)
    JuMP.@variable(model, x[1:t,1:n], Bin)
    JuMP.@variable(model, y[1:t,1:n], Bin)
    JuMP.@variable(model, z[1:t,1:n,1:n], Bin)

    for j in 1:t
        for r in J, s in J
            if r >= s
                continue
            end
            JuMP.@constraints(model, begin
                z[j,r,s] <= x[j,r] + x[j,s]
                z[j,r,s] <= x[j,r] + y[j,r]
                z[j,r,s] <= x[j,s] + y[j,s]
                z[j,r,s] <= y[j,r] + y[j,s]
                z[j,r,s] >= x[j,r] + y[j,s] - 1
                z[j,r,s] >= x[j,s] + y[j,r] - 1
            end)
        end
        for r in J
            JuMP.@constraint(model, x[j,r] + y[j,r] <= 1)
        end
    end

    for r in J, s in J
        if r >= s
            continue
        end
        if E[r,s] == 1
            JuMP.@constraint(model, sum(z[j,r,s] for j in 1:t) == 0)
        else
            JuMP.@constraint(model, sum(z[j,r,s] for j in 1:t) >= 1)
        end
    end

    # Maybe throw in an objctive to maximize spacial logic?
    # JuMP.@objective(model, Min, sum(x) + sum(y))

    stat = JuMP.optimize!(model)

    return termination_status(model) == MOI.OPTIMAL
end

function plot_optimal(E, obstacles, points, lower, upper)
    while upper - lower > 1
        t = (upper + lower) ÷ 2 # integer division
        is_feasible = is_biclique_feasible(E, t)

        if is_feasible
            upper = t
        else
            lower = t
        end
    end

    t = upper
    n = size(E)[1]
    J = 1:n

    model = JuMP.Model(Gurobi.Optimizer)
    JuMP.@variable(model, x[1:t,1:n], Bin)
    JuMP.@variable(model, y[1:t,1:n], Bin)
    JuMP.@variable(model, z[1:t,1:n,1:n], Bin)

    for j in 1:t
        for r in J, s in J
            if r >= s
                continue
            end
            JuMP.@constraints(model, begin
                z[j,r,s] <= x[j,r] + x[j,s]
                z[j,r,s] <= x[j,r] + y[j,r]
                z[j,r,s] <= x[j,s] + y[j,s]
                z[j,r,s] <= y[j,r] + y[j,s]
                z[j,r,s] >= x[j,r] + y[j,s] - 1
                z[j,r,s] >= x[j,s] + y[j,r] - 1
            end)
        end
        for r in J
            JuMP.@constraint(model, x[j,r] + y[j,r] <= 1)
        end
    end

    for r in J, s in J
        if r >= s
            continue
        end
        if E[r,s] == 1
            JuMP.@constraint(model, sum(z[j,r,s] for j in 1:t) == 0)
        else
            JuMP.@constraint(model, sum(z[j,r,s] for j in 1:t) >= 1)
        end
    end

    # Maybe throw in an objctive to maximize spacial logic?
    # JuMP.@objective(model, Min, sum(x) + sum(y))

    stat = JuMP.optimize!(model)

    # x[i, j] <- Element i is in left of biclique j
    # y[i, j] <- Element i is in right of biclique j
    if (termination_status(model) == MOI.OPTIMAL)
        for i in 1:size(x)[1]
            Plots.scatter()

            for j in 1:size(x)[2]
                if (value(x[i, j]) == 1)
                    point = points[j]
                    Plots.scatter!([point.first], [point.second], color="red", lims=(-0.1, 1.1), series_annotations=[Plots.text(string(j), :right, 6, "courier")])
                end

                if (value(y[i, j]) == 1)
                    point = points[j]
                    Plots.scatter!([point.first], [point.second], color="blue", lims=(-0.1, 1.1), series_annotations=[Plots.text(string(j), :right, 6, "courier")])
                end
            end

            ClutteredEnvPathOpt.plot_field(obstacles)
            ClutteredEnvPathOpt.plot_lines(obstacles)

            Plots.savefig(string("optimal", i))
        end
    end

    return termination_status(model) == MOI.OPTIMAL
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Libuuid_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Libuuid")
JLLWrappers.@generate_main_file("Libuuid", UUID("38a345b3-de98-5d2b-a5d3-14cd9215e700"))
end  # module Libuuid_jll


text:
type Request
  fixturename::Union{Void, Symbol}
  param::Any
end

function set_fixturename!(r::Request, name)
  r.fixturename = name
  nothing
end

function set_param!(r::Request, param)
  r.param = param
  nothing
end

Request() = Request(nothing, nothing)

@fixture function request()
  Request()
end


text:
# First Faber Schauder Basis evaluated at time `t`
Λ(t, T::Float64) = sqrt(T)*0.5 - abs((t % T)/sqrt(T) - sqrt(T)*0.5)

# Rescaled Faber Schauder Basis evaluated at time `t`
Λ(t, l⁻::Int64, T::Float64) = Λ(t*(1<<l⁻), T)/sqrt(1<<l⁻)

"""
    dotψ(ξ, s, L, T, u, v)

Given the truncated FS expansion with truncation level `L` and
coefficients `ξ`, output the value of the diffuion bridge at time `s` (`r`)
with initial value `u` at time 0 and final value `v` at `T`.
"""
function dotψ(ξ, s, L, T)
    0 <= s < T || error("out of bounds")
    r = s/sqrt(T)*ξ[end] + sqrt(T)*(1 - s/T)*ξ[1]
    for i in 0:L
        j = floor(Int, s/T * (1 << (L - i)))*(2 << i) + (1 << i) + 1
        r += ξ[j]*Λ(s, L-i, T)
    end
    r
end

"""
    dotψmoving(t, ξ, θ, t′, s, F, L, T, u, v)

Jointly updates the coefficeints (locally) and evaluates the diffuion bridge at time `s`.

Given the truncated FS expansion with truncation level `L` and
coefficients `ξ` and velocities `θ`, move first the coefficients required for
the evaluation of the diffuion bridge at time `s` up to time `t′` (according
to the dynamics of the sampler `F`) and output the
value of diffuion bridge at time `s` with initial value `u` at time 0 and final value `v` at `T`.
"""
function dotψmoving(t, ξ, θ, t′, s, F, L, T)
    0 <= s <= T || error("out of bounds")
    r = s/sqrt(T)*ξ[end] + sqrt(T)*(1 - s/T)*ξ[1]
    for i in 0:L
        j = floor(Int, s/T*(1 << (L - i)))*(2 << i) + (1 << i) + 1
        ZigZagBoomerang.smove_forward!(j, t, ξ, θ, t′, F)
        r += ξ[j]*Λ(s, L-i, T)
    end
    r
end

# find level of index i
"""
    lvl0(i)

Level of element i in the Faber-Schauder base for functions pinned
down `f(0) = f(T) = 0` with elements
ordered according to their midpoint,

```
lvl0.(1:17)
[0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0]
```
"""
function lvl0(i)
    l = 0
    while (i & 1) == 0
        l += 1
        i = i >> 1
    end
    l
end

"""
    lvl(i, L)

Level of element i in the Faber-Schauder base with elements
ordered according to their midpoint.

```
lvl.(1:17, 3)
[3, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 3]
```
"""
function lvl(i, L)
    if i == 1 || i == (2 << L) + 1
        return L
    else
        return lvl0(i-1)
    end
end

"""
    ∇ϕ(ξ, i, K, L, T)

Unbiased estimate for the `i`th partial derivative of the potential function.
The variance of the estimate can be reduced by averaging over `K` independent realization.
`ξ` is the current position of the coefficients, `L` the truncation level.
The bridge has initial value `u` at time 0 and final value `v` at `T`.
"""
function ∇ϕ(ξ, i, K, L, T) # formula (17)
    if i == (2 << L) + 1    # final point
        s = T*(rand())
        x = dotψmoving(t, ξ, θ, t′, s, F, L,  T)
        return -b(ξ[end]*sqrt(T))*sqrt(T) + 0.5*sqrt(T)*s*(2b(x)*b′(x) + b″(x)) + ξ[i] - ξ[1]
    elseif i == 1   # initial point
        s = T*(rand())
        x = dotψmoving(t, ξ, θ, t′, s, F, L,  T)
        return b(ξ[1]*sqrt(T))*sqrt(T) + 0.5*T^(1.5)*(1 - s/T)*(2b(x)*b′(x) + b″(x)) + ξ[i]
    else
        l = lvl(i, L)
        k = (i - 1) ÷ (2 << l)
        δ = T/(1 << (L-l)) # T/(2^(L-l))
        r = 0.0
        for _ in 1:K
            s = δ*(k + rand())
            x = dotψ(ξ, s, L,  T)
            r += 0.5*δ*Λ(s, L-l, T)*(2b(x)*b′(x) + b″(x)) + ξ[i]
        end
        return r/K
    end
end

"""
    ∇ϕmoving(t, ξ, θ, i, t′, F, L, T, u, v)

Jointly updates the coefficeints (locally) and estimates
the `i`th partial derivative of the potential function.
The bridge has initial value `u` at time 0 and final value `v` at `T`.
"""
function ∇ϕmoving(t, ξ, θ, i, t′, F, L, T) # formula (17)
    if i == (2 << L) + 1    # final point
        s = T * (rand())
        x = dotψmoving(t, ξ, θ, t′, s, F, L, T)
        return -b(ξ[end]*sqrt(T))*sqrt(T) + 0.5 * sqrt(T) * s * (2b(x) * b′(x) + b″(x)) + ξ[i] - ξ[1]
    elseif i == 1   # initial point
        s = T * (rand())
        x = dotψmoving(t, ξ, θ, t′, s, F, L, T)
        return 0.5 * T^(1.5) * (1 - s / T) * (2b(x) * b′(x) + b″(x)) + ξ[1]
    else
        l = lvl(i, L)
        k = (i - 1) ÷ (2 << l)
        δ = T / (1 << (L - l))
        s = δ * (k + rand())
        x = dotψmoving(t, ξ, θ, t′, s, F, L, T)
        return 0.5 * δ * Λ(s, L - l, T) * (2b(x) * b′(x) + b″(x)) + ξ[i]
    end
end

"""
    ∇ϕ!(y, ξ, k, L,  T, u, v)

In-place evaluation of the gradient of the potential function.
`ξ` is the current position, `k` is the number of MC realization,
`L` is the truncation level. The bridge has initial value `u` at time 0 and final value `v` at `T`.
"""
function ∇ϕ!(y, ξ, k, L, T)
    for i in eachindex(ξ)
        y[i] = ∇ϕ(ξ, i, k, L, T)
    end
    y
end

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
include("utils.jl")
using StatsBase: AbstractWeights, pweights
using Random: GLOBAL_RNG, seed!
import ManifoldsBase:
    manifold_dimension, exp!, log!, inner, zero_vector!, decorated_manifold, base_manifold
using Manifolds:
    AbstractEstimationMethod,
    CyclicProximalPointEstimation,
    GeodesicInterpolation,
    GeodesicInterpolationWithinRadius,
    GradientDescentEstimation,
    WeiszfeldEstimation
import Manifolds: mean!, median!, var, mean_and_var

struct TestStatsSphere{N} <: AbstractManifold{ℝ} end
TestStatsSphere(N) = TestStatsSphere{N}()
manifold_dimension(::TestStatsSphere{N}) where {N} = manifold_dimension(Sphere(N))
function exp!(::TestStatsSphere{N}, w, x, y; kwargs...) where {N}
    return exp!(Sphere(N), w, x, y; kwargs...)
end
function log!(::TestStatsSphere{N}, y, x, v; kwargs...) where {N}
    return log!(Sphere(N), y, x, v; kwargs...)
end
function inner(::TestStatsSphere{N}, x, v, w; kwargs...) where {N}
    return inner(Sphere(N), x, v, w; kwargs...)
end
function zero_vector!(::TestStatsSphere{N}, v, x; kwargs...) where {N}
    return zero_vector!(Sphere(N), v, x; kwargs...)
end

struct TestStatsEuclidean{N} <: AbstractManifold{ℝ} end
TestStatsEuclidean(N) = TestStatsEuclidean{N}()
manifold_dimension(::TestStatsEuclidean{N}) where {N} = manifold_dimension(Euclidean(N))
function exp!(::TestStatsEuclidean{N}, y, x, v; kwargs...) where {N}
    return exp!(Euclidean(N), y, x, v; kwargs...)
end
function log!(::TestStatsEuclidean{N}, w, x, y; kwargs...) where {N}
    return log!(Euclidean(N), w, x, y; kwargs...)
end
function inner(::TestStatsEuclidean{N}, x, v, w; kwargs...) where {N}
    return inner(Euclidean(N), x, v, w; kwargs...)
end
function zero_vector!(::TestStatsEuclidean{N}, v, x; kwargs...) where {N}
    return zero_vector!(Euclidean(N), v, x; kwargs...)
end

struct TestStatsNotImplementedEmbeddedManifold <:
       AbstractEmbeddedManifold{ℝ,TransparentIsometricEmbedding} end
decorated_manifold(::TestStatsNotImplementedEmbeddedManifold) = Sphere(2)
base_manifold(::TestStatsNotImplementedEmbeddedManifold) = Sphere(2)

struct TestStatsNotImplementedEmbeddedManifold2 <:
       AbstractEmbeddedManifold{ℝ,DefaultIsometricEmbeddingType} end
decorated_manifold(::TestStatsNotImplementedEmbeddedManifold2) = Sphere(2)
base_manifold(::TestStatsNotImplementedEmbeddedManifold2) = Sphere(2)

struct TestStatsNotImplementedEmbeddedManifold3 <:
       AbstractEmbeddedManifold{ℝ,DefaultEmbeddingType} end
decorated_manifold(::TestStatsNotImplementedEmbeddedManifold3) = Sphere(2)
base_manifold(::TestStatsNotImplementedEmbeddedManifold3) = Sphere(2)

function test_mean(M, x, yexp=nothing, method...; kwargs...)
    @testset "mean unweighted" begin
        y = mean(M, x; kwargs...)
        @test is_point(M, y; atol=10^-9)
        if yexp !== nothing
            @test isapprox(M, y, yexp; atol=10^-7)
        end
        y2, _ = mean_and_var(M, x; kwargs...)
        @test isapprox(M, y2, y; atol=10^-7)
        y3, _ = mean_and_std(M, x; kwargs...)
        @test isapprox(M, y3, y; atol=10^-7)
    end

    @testset "mean weighted" begin
        n = length(x)
        w1 = pweights(ones(n) / n)
        w2 = pweights(ones(n))
        w3 = pweights(2 * ones(n))
        y = mean(M, x; kwargs...)
        for w in (w1, w2, w3)
            @test is_point(M, mean(M, x, w; kwargs...); atol=10^-9)
            @test isapprox(M, mean(M, x, w; kwargs...), y)

            @test isapprox(M, mean_and_var(M, x, w; kwargs...)[1], y; atol=10^-7)
            @test isapprox(M, mean_and_std(M, x, w; kwargs...)[1], y; atol=10^-7)
        end
        @test_throws DimensionMismatch mean(M, x, pweights(ones(n + 1)); kwargs...)
    end
    return nothing
end

function test_median(
    M,
    x,
    yexp=nothing;
    method::Union{Nothing,AbstractEstimationMethod}=nothing,
    kwargs...,
)
    @testset "median unweighted$(!isnothing(method) ? " ($method)" : "")" begin
        y = isnothing(method) ? median(M, x; kwargs...) : median(M, x, method; kwargs...)
        @test is_point(M, y; atol=10^-9)
        if yexp !== nothing
            @test isapprox(M, y, yexp; atol=10^-5)
        end
    end

    @testset "median weighted$(!isnothing(method) ? " ($method)" : "")" begin
        n = length(x)
        w1 = pweights(ones(n) / n)
        w2 = pweights(ones(n))
        w3 = pweights(2 * ones(n))
        y = median(M, x; kwargs...)
        for w in (w1, w2, w3)
            if isnothing(method)
                @test is_point(M, median(M, x, w; kwargs...); atol=10^-9)
            else
                @test is_point(M, median(M, x, w, method; kwargs...); atol=10^-9)
            end
            @test isapprox(M, median(M, x, w; kwargs...), y; atol=10^-4)
        end
        if isnothing(method)
            @test_throws Exception median(M, x, pweights(ones(n + 1)); kwargs...)
        else
            @test_throws Exception median(M, x, pweights(ones(n + 1)), method; kwargs...)
        end
    end
    return nothing
end

function test_var(M, x, vexp=nothing; kwargs...)
    n = length(x)
    @testset "var unweighted" begin
        v = var(M, x; kwargs...)
        if vexp !== nothing
            @test v ≈ vexp
        end
        @test v ≈ var(M, x; corrected=true, kwargs...)
        _, v2 = mean_and_var(M, x; kwargs...)
        @test v2 ≈ v
        m = mean(M, x; kwargs...)
        @test var(M, x, m; kwargs...) ≈ var(M, x; kwargs...)
        @test var(M, x; corrected=false, kwargs...) ≈
              var(M, x, m; corrected=false, kwargs...)
        @test var(M, x; corrected=false, kwargs...) ≈ var(M, x; kwargs...) * (n - 1) / n
        @test var(M, x, m; corrected=false, kwargs...) ≈
              var(M, x, m; kwargs...) * (n - 1) / n
    end

    @testset "var weighted" begin
        n = length(x)
        w1 = pweights(ones(n) / n)
        w2 = pweights(ones(n))
        w3 = pweights(2 * ones(n))

        v = var(M, x; corrected=false, kwargs...)
        for w in (w1, w2, w3)
            @test var(M, x, w; kwargs...) ≈ v
            @test var(M, x, w; corrected=false, kwargs...) ≈ v
            @test mean_and_var(M, x, w; kwargs...)[2] ≈ var(M, x, w; kwargs...)
            m = mean(M, x, w; kwargs...)
            @test var(M, x, w, m; kwargs...) ≈ var(M, x, w; kwargs...)
            @test var(M, x, w; corrected=true, kwargs...) ≈
                  var(M, x, w, m; corrected=true, kwargs...)
            @test var(M, x, w; corrected=true, kwargs...) ≈
                  var(M, x, w; kwargs...) * n / (n - 1)
            @test var(M, x, w, m; corrected=true, kwargs...) ≈
                  var(M, x, w, m; kwargs...) * n / (n - 1)
        end
        @test_throws DimensionMismatch var(M, x, pweights(ones(n + 1)); kwargs...)
    end
    return nothing
end

function test_std(M, x, sexp=nothing; kwargs...)
    n = length(x)
    @testset "std unweighted" begin
        s = std(M, x; kwargs...)
        if sexp !== nothing
            @test s ≈ sexp
        end
        @test s ≈ std(M, x; corrected=true, kwargs...)
        @test s ≈ √var(M, x; corrected=true, kwargs...)
        _, s2 = mean_and_std(M, x; kwargs...)
        @test s2 ≈ s
        m = mean(M, x; kwargs...)
        @test std(M, x, m; kwargs...) ≈ std(M, x; kwargs...)
        @test std(M, x; corrected=false, kwargs...) ≈
              std(M, x, m; corrected=false, kwargs...)
        @test std(M, x; corrected=false, kwargs...) ≈
              std(M, x; kwargs...) * sqrt((n - 1) / n)
        @test std(M, x, m; corrected=false, kwargs...) ≈
              std(M, x, m; kwargs...) * sqrt((n - 1) / n)
    end

    @testset "std weighted" begin
        n = length(x)
        w1 = pweights(ones(n) / n)
        w2 = pweights(ones(n))
        w3 = pweights(2 * ones(n))

        s = std(M, x; corrected=false, kwargs...)
        for w in (w1, w2, w3)
            @test std(M, x, w; kwargs...) ≈ s
            @test std(M, x, w; corrected=false, kwargs...) ≈ s
            @test mean_and_std(M, x, w; kwargs...)[2] ≈ std(M, x, w; kwargs...)
            m = mean(M, x, w; kwargs...)
            @test std(M, x, w, m; kwargs...) ≈ std(M, x, w; kwargs...)
            @test std(M, x, w; corrected=true, kwargs...) ≈
                  std(M, x, w, m; corrected=true, kwargs...)
            @test std(M, x, w; corrected=true, kwargs...) ≈
                  std(M, x, w; kwargs...) * sqrt(n / (n - 1))
            @test std(M, x, w, m; corrected=true, kwargs...) ≈
                  std(M, x, w, m; kwargs...) * sqrt(n / (n - 1))
        end
        @test_throws DimensionMismatch std(M, x, pweights(ones(n + 1)); kwargs...)
    end
    return nothing
end

function test_moments(M, x)
    n = length(x)
    @testset "moments unweighted" begin
        m = mean(M, x)
        for i in 1:5
            @test moment(M, x, i) ≈ mean(distance.(Ref(M), Ref(m), x) .^ i)
            @test moment(M, x, i, m) ≈ moment(M, x, i)
        end
        @test moment(M, x, 2) ≈ var(M, x; corrected=false)
        @test skewness(M, x) ≈ moment(M, x, 3) / moment(M, x, 2)^(3 / 2)
        @test kurtosis(M, x) ≈ moment(M, x, 4) / moment(M, x, 2)^2 - 3

        @test moment(M, x, 2, m) ≈ var(M, x; corrected=false)
        @test skewness(M, x, m) ≈ moment(M, x, 3) / moment(M, x, 2)^(3 / 2)
        @test kurtosis(M, x, m) ≈ moment(M, x, 4) / moment(M, x, 2)^2 - 3
    end

    @testset "moments weighted" begin
        n = length(x)
        w1 = pweights(ones(n) / n)
        w2 = pweights(ones(n))
        w3 = pweights(2 * ones(n))

        for w in (w1, w2, w3)
            m = mean(M, x, w)
            for i in 1:5
                @test moment(M, x, i, w) ≈ mean(distance.(Ref(M), Ref(m), x) .^ i, w)
                @test moment(M, x, i, w, m) ≈ moment(M, x, i, w)
            end
            @test moment(M, x, 2, w) ≈ var(M, x, w; corrected=false)
            @test skewness(M, x, w) ≈ moment(M, x, 3, w) / moment(M, x, 2, w)^(3 / 2)
            @test kurtosis(M, x, w) ≈ moment(M, x, 4, w) / moment(M, x, 2, w)^2 - 3

            @test moment(M, x, 2, w, m) ≈ var(M, x, w; corrected=false)
            @test skewness(M, x, w, m) ≈ moment(M, x, 3, w) / moment(M, x, 2, w)^(3 / 2)
            @test kurtosis(M, x, w, m) ≈ moment(M, x, 4, w) / moment(M, x, 2, w)^2 - 3
        end
        @test_throws DimensionMismatch moment(M, x, 3, pweights(ones(n + 1)))
        @test_throws DimensionMismatch skewness(M, x, pweights(ones(n + 1)))
        @test_throws DimensionMismatch kurtosis(M, x, pweights(ones(n + 1)))
    end
    return nothing
end

struct TestStatsOverload1 <: AbstractManifold{ℝ} end
struct TestStatsOverload2 <: AbstractManifold{ℝ} end
struct TestStatsOverload3 <: AbstractManifold{ℝ} end
struct TestStatsMethod1 <: AbstractEstimationMethod end

function mean!(
    ::TestStatsOverload1,
    y,
    ::AbstractVector,
    ::AbstractWeights,
    ::GradientDescentEstimation,
)
    return fill!(y, 3)
end
mean!(::TestStatsOverload2, y, ::AbstractVector, ::AbstractWeights) = fill!(y, 4)
function mean!(
    ::TestStatsOverload2,
    y,
    ::AbstractVector,
    ::AbstractWeights,
    ::GradientDescentEstimation,
)
    return fill!(y, 3)
end
function mean!(
    ::TestStatsOverload3,
    y,
    ::AbstractVector,
    ::AbstractWeights,
    ::TestStatsMethod1=TestStatsMethod1(),
)
    return fill!(y, 5)
end

function median!(
    ::TestStatsOverload1,
    y,
    ::AbstractVector,
    ::AbstractWeights,
    ::CyclicProximalPointEstimation,
)
    return fill!(y, 3)
end
median!(::TestStatsOverload2, y, ::AbstractVector, ::AbstractWeights) = fill!(y, 4)
function median!(
    ::TestStatsOverload2,
    y,
    ::AbstractVector,
    ::AbstractWeights,
    ::CyclicProximalPointEstimation,
)
    return fill!(y, 3)
end
function median!(
    ::TestStatsOverload3,
    y,
    ::AbstractVector,
    ::AbstractWeights,
    ::TestStatsMethod1=TestStatsMethod1(),
)
    return fill!(y, 5)
end

function var(::TestStatsOverload1, ::AbstractVector, ::AbstractWeights, m; corrected=false)
    return 4 + 5 * corrected
end
function mean_and_var(
    ::TestStatsOverload1,
    ::AbstractVector,
    ::AbstractWeights;
    corrected=false,
    kwargs...,
)
    return [4.0], 4 + 5 * corrected
end
function mean_and_var(
    ::TestStatsOverload1,
    ::AbstractVector,
    ::AbstractWeights,
    ::TestStatsMethod1;
    corrected=false,
    kwargs...,
)
    return [5.0], 9 + 7 * corrected
end

@testset "Statistics" begin
    @testset "defaults and overloads" begin
        w = pweights([2.0])
        x = [[0.0]]
        @testset "mean" begin
            M = TestStatsOverload1()
            @test mean(M, x) == [3.0]
            @test mean(M, x, w) == [3.0]
            @test mean(M, x, w, GradientDescentEstimation()) == [3.0]
            @test mean(M, x, GradientDescentEstimation()) == [3.0]
            M = TestStatsOverload2()
            @test mean(M, x) == [4.0]
            @test mean(M, x, w) == [4.0]
            @test mean(M, x, w, GradientDescentEstimation()) == [3.0]
            @test mean(M, x, GradientDescentEstimation()) == [3.0]
            M = TestStatsOverload3()
            @test mean(M, x) == [5.0]
            @test mean(M, x, w) == [5.0]
        end

        @testset "median" begin
            M = TestStatsOverload1()
            @test median(M, x) == [3.0]
            @test median(M, x, w) == [3.0]
            @test median(M, x, w, CyclicProximalPointEstimation()) == [3.0]
            @test median(M, x, CyclicProximalPointEstimation()) == [3.0]
            M = TestStatsOverload2()
            @test median(M, x) == [4.0]
            @test median(M, x, w) == [4.0]
            @test median(M, x, w, CyclicProximalPointEstimation()) == [3.0]
            @test median(M, x, CyclicProximalPointEstimation()) == [3.0]
            M = TestStatsOverload3()
            @test median(M, x) == [5.0]
            @test median(M, x, w) == [5.0]
        end

        @testset "var" begin
            M = TestStatsOverload1()
            @test mean_and_var(M, x) == ([4.0], 9)
            @test mean_and_var(M, x, w) == ([4.0], 4)
            @test mean_and_std(M, x) == ([4.0], 3.0)
            @test mean_and_std(M, x, w) == ([4.0], 2.0)
            @test var(M, x) == 9
            @test var(M, x, 2) == 9
            @test var(M, x, w) == 4
            @test var(M, x, w, 2) == 4
            @test std(M, x) == 3.0
            @test std(M, x, 2) == 3.0
            @test std(M, x, w) == 2.0
            @test std(M, x, w, 2) == 2.0

            @test mean_and_var(M, x, TestStatsMethod1()) == ([5.0], 16)
            @test mean_and_var(M, x, w, TestStatsMethod1()) == ([5.0], 9)
            @test mean_and_std(M, x, TestStatsMethod1()) == ([5.0], 4.0)
            @test mean_and_std(M, x, w, TestStatsMethod1()) == ([5.0], 3.0)
        end
    end

    @testset "decorator dispatch" begin
        ps = [normalize([1, 0, 0] .+ 0.1 .* randn(3)) for _ in 1:3]
        M1 = TestStatsNotImplementedEmbeddedManifold()
        @test mean!(M1, similar(ps[1]), ps) == mean!(Sphere(2), similar(ps[1]), ps)
        @test mean(M1, ps) == mean(Sphere(2), ps)
        @test median!(M1, similar(ps[1]), ps) == median!(Sphere(2), similar(ps[1]), ps)
        @test median(M1, ps) == median(Sphere(2), ps)

        M2 = TestStatsNotImplementedEmbeddedManifold2()
        @test_throws ErrorException mean(M2, ps)
        @test_throws ErrorException mean!(M2, similar(ps[1]), ps)
        @test_throws ErrorException median(M2, ps)
        @test_throws ErrorException median!(M2, similar(ps[1]), ps)

        M3 = TestStatsNotImplementedEmbeddedManifold3()
        @test_throws ErrorException mean(M3, ps)
        @test_throws ErrorException mean!(M3, similar(ps[1]), ps)
        @test_throws ErrorException median(M3, ps)
        @test_throws ErrorException median!(M3, similar(ps[1]), ps)
    end

    @testset "TestStatsSphere" begin
        M = TestStatsSphere(2)

        @testset "consistency" begin
            p = [0.0, 0.0, 1.0]
            n = 3
            x = [
                exp(M, p, π / 6 * [cos(α), sin(α), 0.0]) for
                α in range(0, 2 * π - 2 * π / n, length=n)
            ]
            test_mean(M, x)
            test_median(M, x; atol=1e-12)
            test_median(M, x; method=CyclicProximalPointEstimation(), atol=1e-12)
            test_median(M, x; method=WeiszfeldEstimation())
            test_var(M, x)
            test_std(M, x)
            test_moments(M, x)
        end

        @testset "zero variance" begin
            x = fill([0.0, 0.0, 1.0], 5)
            test_mean(M, x, [0.0, 0.0, 1.0])
            test_median(M, x, [0.0, 0.0, 1.0]; atol=10^-12)
            test_var(M, x, 0.0)
            test_std(M, x, 0.0)
        end

        @testset "three points" begin
            x = [[1.0, 0.0, 0.0], [1.0, 1.0, 0.0] / √2, [0.0, 1.0, 0.0]]
            θ = π / 4
            @test isapprox(M, mean(M, x), geodesic(M, x[1], x[3], θ))
            test_mean(M, x, [1.0, 1.0, 0.0] / √2)
            test_median(M, x, [1.0, 1.0, 0.0] / √2; atol=10^-12)
            test_var(M, x, θ^2)
            test_std(M, x, θ)
            test_moments(M, x)
        end
    end

    @testset "TestStatsEuclidean{N}" begin
        @testset "N=1" begin
            M = TestStatsEuclidean(1)

            rng = MersenneTwister(42)
            x = [randn(rng, 1) for _ in 1:10]
            vx = vcat(x...)

            test_mean(M, x, mean(x))
            test_median(M, x; atol=10^-12)
            test_var(M, x, var(vx))
            test_std(M, x, std(vx))
            test_moments(M, x)

            w = pweights(rand(rng, 10))
            @test mean(M, x, w) ≈ [mean(vx, w)]
            @test var(M, x, w) ≈ var(vx, w)
            @test std(M, x, w) ≈ std(vx, w)
        end

        @testset "N=5" begin
            M = TestStatsEuclidean(5)

            rng = MersenneTwister(42)
            x = [randn(rng, 5) for _ in 1:10]
            test_mean(M, x, mean(x))
            test_var(M, x, sum(var(x)))
            test_std(M, x, sqrt(sum(std(x) .^ 2)))
            test_moments(M, x)

            w = pweights(rand(rng, 10))
            ax = hcat(x...)
            @test mean(M, x, w) ≈ mean(ax, w; dims=2)
            @test var(M, x, w) ≈ sum(var(ax, w, 2))
            @test std(M, x, w) ≈ sqrt(sum(std(ax, w, 2) .^ 2))
        end
    end

    @testset "Euclidean statistics" begin
        @testset "N=1" begin
            M = Euclidean()
            @testset "scalar" begin
                x = [1.0, 2.0, 3.0, 4.0]
                w = pweights(ones(length(x)) / length(x))
                @test mean(M, x) ≈ mean(x)
                @test mean(M, x, w) ≈ mean(x, w)
                @test median(M, x; rng=MersenneTwister(1212), atol=10^-12) ≈ median(x)
                @test median(M, x, w; rng=MersenneTwister(1212), atol=10^-12) ≈ median(x, w)
                @test var(M, x) ≈ var(x)
                @test var(M, x, w) ≈ var(x, w)
                @test std(M, x) ≈ std(x)
                @test std(M, x, w) ≈ std(x, w)

                test_mean(M, x)
                test_median(M, x; atol=10^-12)
                test_var(M, x)
                test_std(M, x)
                test_moments(M, x)
            end

            @testset "vector" begin
                x = [[1.0], [2.0], [3.0], [4.0]]
                vx = vcat(x...)
                w = pweights(ones(length(x)) / length(x))
                @test mean(M, x) ≈ mean(x)
                @test mean(M, x, w) ≈ [mean(vx, w)]
                @test var(M, x) ≈ var(vx)
                @test var(M, x, w) ≈ var(vx, w)
                @test std(M, x) ≈ std(vx)
                @test std(M, x, w) ≈ std(vx, w)

                test_mean(M, x)
                test_median(M, x; atol=10^-12)
                test_var(M, x)
                test_std(M, x)
                test_moments(M, x)
            end
        end
    end

    @testset "GeodesicInterpolation" begin
        @testset "mean" begin
            rng = MersenneTwister(1212)
            @testset "exact for Euclidean" begin
                M = Euclidean(2, 2)
                x = [randn(rng, 2, 2) for _ in 1:10]
                w = pweights([rand(rng) for _ in 1:10])
                @test mean(M, x, GeodesicInterpolation()) ≈ mean(x)
                @test mean(M, x, w, GeodesicInterpolation()) ≈
                      mean(M, x, w, GradientDescentEstimation())
            end

            @testset "three points" begin
                S = Sphere(2)
                x = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0] / √2]
                @test mean(S, [x[1], x[1]], GeodesicInterpolation()) ≈ x[1]
                @test mean(S, x[1:2], GeodesicInterpolation()) ≈ x[3]
                @test mean(S, x[1:2], pweights([1, 0]), GeodesicInterpolation()) ≈ x[1]
                @test mean(S, x[1:2], pweights([0, 1]), GeodesicInterpolation()) ≈ x[2]
                @test mean(S, x[1:2], pweights([1, 2]), GeodesicInterpolation()) ≈
                      shortest_geodesic(S, x[1], x[2], 2 / 3)
                m = mean(S, x, pweights([1, 2, 3]), GeodesicInterpolation())
                @test m ≈ shortest_geodesic(
                    S,
                    shortest_geodesic(S, x[1], x[2], 2 / 3),
                    x[3],
                    1 / 2,
                )
            end

            @testset "resumable" begin
                S = Sphere(2)
                x = [normalize(randn(rng, 3)) for _ in 1:10]
                w = pweights([rand(rng) for _ in 1:10])
                ypart = mean(S, x[1:5], pweights(w[1:5]), GeodesicInterpolation())
                yfull = mean(S, x, w, GeodesicInterpolation())
                x2 = [[ypart]; x[6:end]]
                w2 = pweights([sum(w[1:5]); w[6:end]])
                @test mean(S, x2, w2, GeodesicInterpolation()) ≈ yfull
            end
        end

        @testset "welford" begin
            rng = MersenneTwister(56)

            @testset "exact for Euclidean" begin
                M = TestStatsEuclidean(4)
                x = [randn(rng, 4) for _ in 1:10]
                w = pweights([rand(rng) for _ in 1:10])
                m1, v1 = mean_and_var(M, x, GeodesicInterpolation())
                @test m1 ≈ mean(x)
                @test v1 ≈ sum(var(x))
                m2, v2 = mean_and_var(M, x, w, GeodesicInterpolation())
                @test m2 ≈ mean(M, x, w, GradientDescentEstimation())
                @test v2 ≈ var(M, x, w)
            end

            @testset "three points" begin
                S = Sphere(2)
                x = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0] / √2]
                m1, v1 = mean_and_var(S, [x[1], x[1]], GeodesicInterpolation())
                @test m1 ≈ x[1]
                @test v1 ≈ 0 atol = 1e-6
                m2, v2 = mean_and_var(S, x[1:2], GeodesicInterpolation())
                @test m2 ≈ x[3]
                @test v2 ≈ (π / 4)^2 * 2 / (2 - 1)
                m3, v3 = mean_and_var(S, x[1:2], pweights([1, 0]), GeodesicInterpolation())
                @test m3 ≈ x[1]
                @test v3 ≈ 0 atol = 1e-6
                m4, v4 = mean_and_var(S, x[1:2], pweights([0, 1]), GeodesicInterpolation())
                @test m4 ≈ x[2]
                @test v4 ≈ 0 atol = 1e-6
                m5, v5 = mean_and_var(S, x, pweights([1, 2, 3]), GeodesicInterpolation())
                @test m5 ≈ shortest_geodesic(
                    S,
                    shortest_geodesic(S, x[1], x[2], 2 / 3),
                    x[3],
                    1 / 2,
                )
                @test v5 ≈ var(S, x, pweights([1, 2, 3]), m5)
            end

            @testset "within radius" begin
                rng = MersenneTwister(47)
                @test repr(GeodesicInterpolationWithinRadius(Inf)) ==
                      "GeodesicInterpolationWithinRadius(Inf)"
                @test_throws DomainError GeodesicInterpolationWithinRadius(-0.1)
                S = Sphere(2)
                x = [[1.0, 0, 0], [0, 1.0, 0]]
                m = mean(S, x, GeodesicInterpolation())
                mg = mean(S, x, GradientDescentEstimation(); p0=m)
                vg = var(S, x, mg)

                @test mean(S, x, GeodesicInterpolationWithinRadius(Inf)) == m
                @test mean(S, x, GeodesicInterpolationWithinRadius(π)) == m
                @test mean(S, x, GeodesicInterpolationWithinRadius(π / 8)) != m
                @test mean(S, x, GeodesicInterpolationWithinRadius(π / 8)) == mg

                m, v = mean_and_var(S, x, GeodesicInterpolation())
                @test mean_and_var(S, x, GeodesicInterpolationWithinRadius(Inf)) == (m, v)
                @test mean_and_var(S, x, GeodesicInterpolationWithinRadius(π)) == (m, v)
                @test mean_and_var(S, x, GeodesicInterpolationWithinRadius(π / 8)) != (m, v)
                @test mean_and_var(S, x, GeodesicInterpolationWithinRadius(π / 8)) ==
                      (mg, vg)
            end
        end

        @testset "SymmetricPositiveDefinite default" begin
            rng = MersenneTwister(36)
            P1 = SymmetricPositiveDefinite(3)
            P2 = MetricManifold(P1, LinearAffineMetric())
            @testset "$P" for P in [P1, P2]
                p0 = collect(exp(Symmetric(randn(rng, 3, 3) * 0.1)))
                x = [exp(P, p0, Symmetric(randn(rng, 3, 3) * 0.1)) for _ in 1:10]
                w = pweights([rand(rng) for _ in 1:length(x)])
                m = mean(P, x, w)
                mg = mean(P, x, w, GeodesicInterpolation())
                mf = mean(P, x, w, GradientDescentEstimation())
                @test m == mg
                @test m != mf
            end
        end

        @testset "Sphere default" begin
            rng = MersenneTwister(47)
            S = Sphere(2)
            p0 = [1.0, 0, 0]
            x = [normalize(randn(rng, 3)) for _ in 1:10]
            x = [x; -x]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(S, x, w)
            mg = mean(S, x, w, GeodesicInterpolation())
            mf = mean(S, x, w, GradientDescentEstimation(); p0=mg)
            @test m != mg
            @test m == mf

            μ = randn(rng, 3) .* 10
            x = [normalize(randn(rng, 3) .+ μ) for _ in 1:10]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(S, x, w)
            mg = mean(S, x, w, GeodesicInterpolation())
            @test m == mg
        end

        @testset "ProjectiveSpace default" begin
            rng = MersenneTwister(47)
            M = ProjectiveSpace(2)
            p0 = [1.0, 0, 0]
            x = [normalize(randn(rng, 3)) for _ in 1:10]
            x = [x; -x]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(M, x, w)
            mg = mean(M, x, w, GeodesicInterpolation())
            mf = mean(M, x, w, GradientDescentEstimation(); p0=mg)
            @test !isapprox(M, m, mg)
            @test isapprox(M, m, mf)

            μ = randn(rng, 3) .* 10
            x = [normalize(randn(rng, 3) .+ μ) for _ in 1:10]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(M, x, w)
            mg = mean(M, x, w, GeodesicInterpolation())
            @test m == mg
        end

        @testset "Rotations default" begin
            rng = MersenneTwister(47)
            R = Manifolds.Rotations(3)
            p0 = collect(Diagonal(ones(3)))
            v1 = hat(R, p0, [1.0, 0.0, 0.0])
            v2 = hat(R, p0, [0.0, 1.0, 0.0])
            x = [geodesic(R, p0, v1, π / 2 * (1:4)); geodesic(R, p0, v2, π / 2 * (1:4))]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(R, x, w)
            mg = mean(R, x, w, GeodesicInterpolation())
            mf = mean(R, x, w, GradientDescentEstimation(); p0=mg)
            @test m != mg
            @test m == mf

            μ = project(R, randn(3, 3))
            d = normal_tvector_distribution(R, μ, 0.1)
            x = [exp(R, μ, rand(rng, d)) for _ in 1:10]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(R, x, w)
            mg = mean(R, x, w, GeodesicInterpolation())
            @test m == mg
        end

        @testset "Grassmann default" begin
            rng = MersenneTwister(47)
            G = Manifolds.Grassmann(3, 2)
            p0 = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            x = [exp(G, p0, project(G, p0, randn(rng, 3, 2) * 10)) for _ in 1:10]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(G, x, w)
            mg = mean(G, x, w, GeodesicInterpolation())
            mf = mean(G, x, w, GradientDescentEstimation(); p0=mg)
            @test m != mg
            @test m == mf

            x = [exp(G, p0, project(G, p0, randn(rng, 3, 2) * 0.01)) for _ in 1:10]
            w = pweights([rand(rng) for _ in 1:length(x)])
            m = mean(G, x, w)
            mg = mean(G, x, w, GeodesicInterpolation())
            @test m == mg
        end
    end

    @testset "Extrinsic mean" begin
        rng = MersenneTwister(47)
        S = Sphere(2)
        x = [normalize(randn(rng, 3)) for _ in 1:10]
        w = pweights([rand(rng) for _ in 1:length(x)])
        m = normalize(mean(reduce(hcat, x), w; dims=2)[:, 1])
        mg = mean(S, x, w, ExtrinsicEstimation())
        @test isapprox(S, m, mg)
    end

    @testset "Extrinsic median" begin
        rng = MersenneTwister(47)
        S = Sphere(2)
        x = [normalize(randn(rng, 3)) for _ in 1:10]
        w = pweights([rand(rng) for _ in 1:length(x)])
        m = normalize(median(Euclidean(3), x, w))
        mg = median(S, x, w, ExtrinsicEstimation())
        @test isapprox(S, m, mg)
    end

    @testset "Covariance matrix, Euclidean" begin
        rng = MersenneTwister(47)
        M = Euclidean(3)
        x = [randn(rng, 3) for _ in 1:10]
        @test isapprox(cov(M, x), cov(x))
        covest = SimpleCovariance(; corrected=false)
        @test isapprox(cov(M, x; tangent_space_covariance_estimator=covest), cov(covest, x))
    end
    @testset "Covariance matrix, sphere" begin
        rng = MersenneTwister(47)
        S = Sphere(2)
        x = [normalize(randn(rng, 3)) for _ in 1:10]
        covm = cov(S, x)
        @test size(covm) == (2, 2)
        @test isposdef(covm)
        @test issymmetric(covm)
    end
end


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using Cxx

cxxinclude("TFitResultPtr.h")

import Base: show

const TFitResultPtr = cxxt"TFitResultPtr"

function show(io::IO, x::TFitResultPtr)
    print(io, "TFitResultPtr(")
    show(io, reinterpret(Culong, @cxx x->Get()))
    print(io, ")")
end