text:
using LinearAlgebra
# Testing the following distribution: WishartGaussian, DirichletMultinomial, GammaNormal, NormalNormal, BetaBernoulli

@testset "testing posterior paramters" begin

  @testset "Normal-Inverse-Wishart" begin

    # prior
    μ0 = zeros(2)
    κ0 = 5.0
    ν0 = 9.0
    Σ0 = Matrix{Float64}(I, 2, 2)

    # data
    x = [0.2 0.1; 0.1 0.0; 0.2 0.05]
    (N, D) = size(x)

    # distribution
    d = WishartGaussian(μ0, κ0, ν0, Σ0)

    # test prior
    (μ, κ, ν, Σ) = BayesianNonparametrics.posteriorParameters(d)

    @test κ == κ0
    @test ν == ν0
    @test all( μ .== μ0 )
    @test all( Σ .== Σ0 )

    # add data
    BayesianNonparametrics.add!(d, x)

    # test posterior parameters
    (μ, κ, ν, Σ) = BayesianNonparametrics.posteriorParameters(d)

    mu = vec(mean(x, dims=1))
    SS = x' * x

    @test κ == κ0 + N
    @test ν == ν0 + N
    @test all( μ .≈ (κ0 * μ0 .+ mu * N) ./ κ )
    @test all( Σ .== Σ0 + SS - κ * (μ * μ') + κ0 * (μ0 * μ0') )

    # remove data
    BayesianNonparametrics.remove!(d, x)

    # test prior
    (μ, κ, ν, Σ) = BayesianNonparametrics.posteriorParameters(d)

    @test κ == κ0
    @test ν == ν0
    @test all( μ .== μ0 )
    @test all( Σ .== Σ0 )
  end

  @testset "Normal-Gamma" begin

    # prior
    μ0 = 0.0
    λ0 = 1.0
    α0 = 1.0
    β0 = 1.0

    # data
    x = [0.1, 0.0, 0.05]
    N = length(x)

    # distribution
    d = GammaNormal(μ0 = μ0, λ0 = λ0, α0 = α0, β0 = β0)

    # test prior
    (μ, λ, α, β) = BayesianNonparametrics.posteriorParameters(d)

    @test μ == μ0
    @test λ == λ0
    @test α == α0
    @test β == β0

    # add data
    BayesianNonparametrics.add!(d, x)

    # test posterior paramters
    # see: https://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf page 8.
    (μ, λ, α, β) = BayesianNonparametrics.posteriorParameters(d)

    @test λ == λ0 + N
    @test α == α0 + (N / 2)
    @test μ == (λ0 * μ0 + N * mean(x)) / (λ0 + N)
    @test β == β0 + 1/2 * sum( (x .- mean(x)).^2 ) + ( λ0 * N * (mean(x) - μ0)^2 ) / (2 * λ)

    # remove data
    BayesianNonparametrics.remove!(d, x)

    # test prior
    (μ, λ, α, β) = BayesianNonparametrics.posteriorParameters(d)

    @test μ == μ0
    @test λ == λ0
    @test α == α0
    @test β == β0
  end

  @testset "Beta-Bernoulli" begin

    # prior
    α0 = 1.0
    β0 = 1.0

    # data
    x = Bool[0, 1, 0, 0, 1]
    N = length(x)

    # distribution
    d = BetaBernoulli(α0 = α0, β0 = β0)

    # test prior
    (α, β) = BayesianNonparametrics.posteriorParameters(d)

    @test α == α0
    @test β == β0

    # add data
    BayesianNonparametrics.add!(d, x)

    # test posterior paramters
    (α, β) = BayesianNonparametrics.posteriorParameters(d)

    @test α == α0 + sum(x)
    @test β == β0 + (N - sum(x))

    # remove data
    BayesianNonparametrics.remove!(d, x)

    # test prior
    (α, β) = BayesianNonparametrics.posteriorParameters(d)

    @test α == α0
    @test β == β0
  end

end


text:
#///////////////////////////////////////
#// File Name: state_transition.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/15
#// Description: State transition model for Risk Sensitive Stochastic SAC
#///////////////////////////////////////

using LinearAlgebra
using Random

# Instantaneous Transition.
function transition(e_state::RobotState,
                    u::Vector{Float64})
    @assert length(u) == 2 "Invalid control input dimension!"
    ev = get_velocity(e_state);
    return [ev; u]
end

function transition_jacobian(ex_vec::Vector{Float64},
                             u::Vector{Float64})
    @assert length(ex_vec) == 4 "Invalid robot state dimension!"
    @assert length(u) == 2 "Invalid control input dimension!"
    A = [0.0 0.0 1.0 0.0;
         0.0 0.0 0.0 1.0;
         0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0];
    return A
end

function transition_jacobian(e_state::RobotState,
                             u::Vector{Float64})
    return transition_jacobian(e_state.x, u)
end

function transition_control_coeff(e_state::RobotState)
    H = [0.0 0.0;
         0.0 0.0;
         1.0 0.0;
         0.0 1.0];
    return H
end

# Forward Euler Integration of Instataneous Transition
function transition(e_state::RobotState,
                    u::Vector{Float64},
                    dt::Real)
    x_new = e_state.x + transition(e_state, u)*dt;
    return RobotState(x_new, e_state.t + Duration(dt))
end

ext: jl


text:
using GageCompuscope
using Interact
using Blink
using NamedTupleTools
using Interact: hbox, vbox, vline, hline, pad
gage = GageCard();

win = Window();

oAcq = Observable(gage.acquisition_config)
oCh1 = Observable(gage.channel_config[1])
oCh2 = Observable(gage.channel_config[2])

Acq=gage.acquisition_config
#-------------------------------------------------------------------------------
# Acquistion Configuation.
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# Channel 1
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# Channel 2
#-------------------------------------------------------------------------------

body!(win, hbox(acq_ctl, vline(), c1Ctl, c2Ctl, vline()));

max_stars_repo_name: pacific-zitter/GingerGUI

max_issues_repo_name: pacific-zitter/GingerGUI

max_forks_repo_name: pacific-zitter/GingerGUI


text:
using Documenter, SellmeierFit

makedocs(
    modules = [SellmeierFit],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Wonseok Shin",
    sitename = "SellmeierFit.jl",
    pages = Any[
        "Home" => "index.md"
        "Installation" => "install.md"
        "Tutorials" => "tutorials.md"
        "Types and Functions" => "types_funs.md"
    ]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/wsshin/SellmeierFit.jl.git",
    devbranch = "main",  # see https://github.com/JuliaDocs/Documenter.jl/issues/1443#issuecomment-735905544
    push_preview = true
)


text:
#ДАНО: Робот - в юго-западном углу поля, на котором расставлено некоторое количество маркеров

#РЕЗУЛЬТАТ: Функция вернула значение средней температуры всех замаркированных клеток

function sum_temp(r::Robot)
    sum = 0
    count = 0
    side = Ost
    while isborder(r,Nord) == false
        s, c=moves!(r,side)
        sum+=s
        count+=c
        side = inverse(side)
    end
    s, c=moves!(r,side)
    sum+=s
    count+=c
    return (sum/count)
end
function moves!(r::Robot,side::HorizonSide)
    sum = 0
    count = 0
    while isborder(r,side) == false
        if ismarker(r) == true
            count+=1
            sum+= temperature(r)
        end
        move!(r,side)
    end
    if isborder(r,Nord) == false
        move!(r,Nord)
    end
    return sum, count
end
inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))


text:
using Distributed
@everywhere using Combinatorics
@everywhere using Random
@everywhere include("./Grid.jl")
@everywhere using .Grid
@everywhere struct Extremes
    min_share::Float64
    max_share::Float64
    min_voters::Set{Array{Int}}
    max_voters::Set{Array{Int}}
    voters_hash::UInt64
end
include("./dedupe.jl")

const batch_size = 250000
const channel_size = 1
const n_voters = 12

@everywhere function eval_voter_dists(jobs, results)
    while true
        voter_dists = take!(jobs)
        max_share = 0
        min_share = 1
        max_voters = Set([])
        min_voters = Set([])
        v_hash = hash(Set(v for v in eachcol(voter_dists)))
        for i in 1:size(voter_dists, 2)
            voters = voter_dists[:, i]
            share = expected_seat_share(voters)
            if share > max_share
                max_share = share
                max_voters = Set([voters])
            elseif share == max_share
                push!(max_voters, voters)
            end
            if share < min_share
                min_share = share
                min_voters = Set([voters])
            elseif share == min_share
                push!(min_voters, voters)
            end
        end
        extremes = Extremes(min_share, max_share, min_voters, max_voters)#, v_hash)
        put!(results, extremes)
    end
end

function add_combos(jobs::RemoteChannel, meta::RemoteChannel, batch_size::Int, n_voters::Int)
    batch_count = 0
    batch_idx = 1
    batch = zeros(Int, n_voters, batch_size)
    for c in combinations(1:ensemble.size, n_voters)
        if !combo_duplicated(c, ensemble.width)
            batch[:, batch_idx] = c
            batch_idx += 1
            if batch_idx > batch_size
                put!(jobs, copy(batch))
                batch_idx = 1
                batch_count += 1
            end
        end
    end
    # Handle remaining partial batch
    put!(jobs, copy(batch[:, 1:batch_idx - 1]))
    batch_count += 1
    put!(meta, batch_count)
end

function main()
    jobs = RemoteChannel(()->Channel{Array{Int, 2}}(channel_size))
    meta = RemoteChannel(()->Channel{Int}(1))
    results = RemoteChannel(()->Channel{Extremes}(channel_size))
    # Start workers.
    for p in workers()
        remote_do(eval_voter_dists, p, jobs, results)
    end

    max_share = 0
    min_share = 1
    max_voters = Set([])
    min_voters = Set([])
    @async add_combos(jobs, meta, batch_size, n_voters)
    i = 0
    n_batches = missing
    while true
        batch_extremes = take!(results)
        println("[batch ", i, "]")
        println("\tmin share: ", batch_extremes.min_share)
        println("\tmax share: ", batch_extremes.max_share)
        println("\tvoters hash: ", batch_extremes.voters_hash)
        println("\tmin voter dists:")
        for v in batch_extremes.min_voters
            println("\t\t", v)
        end
        println("\tmax voter dists:")
        for v in batch_extremes.max_voters
            println("\t\t", v)
        end

        if batch_extremes.min_share < min_share
            min_share = batch_extremes.min_share
            min_voters = batch_extremes.min_voters
        elseif batch_extremes.min_share == min_share
            min_voters = union(min_voters, batch_extremes.min_voters)
        end
        if batch_extremes.max_share > max_share
            max_share = batch_extremes.max_share
            max_voters = batch_extremes.max_voters
        elseif batch_extremes.max_share == max_share
            max_voters = union(max_voters, batch_extremes.max_voters)
        end
        i += 1
       if isready(meta)
            n_batches = take!(meta)
            println("got number of batches: ", n_batches)
        end
        if !ismissing(n_batches) && n_batches == i
            break  # done!
        end
    end
    println("min share: ", min_share)
    println("max share: ", max_share)
    println("min voter dists:")
    for v in min_voters
        println(v)
    end
    println("max voter dists:")
    for v in max_voters
        println(v)
    end

    # Kill workers.
    for wid in workers()
        if wid != 1
            rmprocs(wid)
        end
    end
end

main()


text:
using LinearAlgebra
using Plots; pyplot()
using JLD, FileIO
using DataFrames, CSV
using LsqFit
using LaTeXStrings

include("fftFilter2d.jl")

function tau(z; dl=5)
    """
    returns Lp/L, averaged over every dl lines
    """
    L =size(z,1)
    Lp  = 0.0
    count=0
    for j=1:dl:L # eval Lp for each line
        count+=1
        for i=1:L-1
            v1 = [1; 0; z[j,i+1] - z[j,i]]
            Lp += norm(v1)
        end
    end
    Lp /=count
    return Lp/L
end

#sa = make2DSurfs(; maxlevel=10, hursts=0.3:0.05:0.85, samples=1:10, atol=0.01)

files = ["level8.jld", "level9.jld", "level10.jld"]
ss = [:l8, :l9, :l10]
df = DataFrame()
df2 = DataFrame()
for (n,file) in enumerate(files)
    sa=load(file)
    hurst = sa["hursts"]
    if n==1
        df[:h] = hurst
        df2[:h] = hurst
    end

    τ = []
    σ = []
    nsamples = size(sa["samples"],1)
    for h in hurst
        τav = 0.0
        sav = 0.0
        for s in sa["samples"]
            z = sa["$h"]["$s"]
            τav += tau(z)
            sav += std(z)
        end
        push!(τ, τav/nsamples)
        push!(σ, sav/nsamples)
    end
    df[!,ss[n]] = τ
    df2[!,ss[n]] = σ
end
df

p1 = plot(size=(400,300),
              legendfontsize=8,
              legend=:inside,
              bottommargin=50,
              fg_legend=:lightgray,
              dpi=150, frame=:box,
              xaxis=(L"H", font(11)),
              yaxis=(L"τ", font(11)),
              grid=false)
@. model(x,p) = sqrt(1.0 + (p[1]^2)*(4*p[2]/p[3])^(2*x))
x = df.h
dx = []
c = []
N = [1, 2^8, 2^9,  2^10]
for (n,col) in enumerate([1,df.l8, df.l9, df.l10])
    if n>2
        y = col
        ℓ =  N[n]
        p=curve_fit(model,x,y, [5.0, 1.0, N[n]], lower=[0.0, 1.0, 0.999999*ℓ], upper=[5, 2, 1.0*ℓ])
        println(N[n], ",", p.param)
        push!(dx, p.param[2])
        push!(c, p.param[1])

        xfit = range(0.29, stop=0.9, length=100)
        yfit = model(xfit, p.param)

        scatter!(x, y, marker=:auto, markersize=8, label="$ℓ")
        plot!(xfit, yfit, linestyle=:dash, color=:black, label="")

        xfit2 = df2.h
        yfit2 = sqrt.(1 .+ (df2.l9).^2 .* ((4/N[n]).^(2*df.h)))
        plot!(xfit2, yfit2, color="#aaaaaa", label="")

    end
end
p1
savefig("/home/humberto/OneD/Documentos/Petrobras/figs/tau.png")

df2.l9


text:
function throw_error_message(wasm_error::Ptr{wasmtime_error_t})
    message_buf = WasmByteVec()
    wasmtime_error_message(
        wasm_error,
        message_buf,
    )
    error_msg = unsafe_string(message_buf.data, message_buf.size)
    @ccall libwasmtime.wasmtime_error_delete(wasm_error::Ptr{wasmtime_error_t})::Cvoid

    error(error_msg)
end

macro wt_check(call)
    quote
        wt_err = $(esc(call))
        if wt_err != C_NULL
            throw_error_message(wt_err)
        end
    end
end


text:
# configure logger; see: https://docs.julialang.org/en/v1/stdlib/Logging/index.html
# and: https://github.com/oxinabox/LoggingExtras.jl
logger = FileLogger("/var/tmp/stb_log.txt")
global_logger(logger)
LogLevel( Logging.Debug )


text:
immutable Chisq <: ContinuousUnivariateDistribution
    df::Float64 # non-integer degrees of freedom are meaningful
    function Chisq(d::Real)
        d > zero(d) || error("df must be positive")
        new(float64(d))
    end
end

@_jl_dist_1p Chisq chisq

@continuous_distr_support Chisq 0.0 Inf

function entropy(d::Chisq)
    x = d.df / 2.0 + log(2.0) + lgamma(d.df / 2.0)
    x + (1.0 - d.df / 2.0) * digamma(d.df / 2.0)
end

kurtosis(d::Chisq) = 12.0 / d.df

mean(d::Chisq) = d.df

# TODO: Switch to using quantile?
function median(d::Chisq)
    k = d.df
    k * (1.0 - 2.0 / (9.0 * k))^3
end

function mgf(d::Chisq, t::Real)
    k = d.df
    (1.0 - 2.0 * t)^(-k / 2.0)
end

cf(d::Chisq, t::Real) = (1.0 - 2.0 * im * t)^(-d.df / 2.0)

mode(d::Chisq) = d.df > 2.0 ? d.df - 2.0 : 0.0

function gradloglik(d::Chisq, x::Float64)
  insupport(Chisq, x) ? (d.df / 2.0 - 1) / x - 0.5 : 0.0
end

# rand - the distribution chi^2(df) is 2 * gamma(df / 2)
# for integer n, a chi^2(n) is the sum of n squared standard normals
function rand(d::Chisq)
    d.df == 1 ? randn()^2 : 2.0 * rand(Gamma(d.df / 2.0))
end

function rand!(d::Chisq, A::Array{Float64})
    if d.df == 1
        for i = 1:length(A)
            @inbounds A[i] = randn()^2
        end
    else
        s = GammaSampler(d.df / 2.0)
        for i = 1:length(A)
            @inbounds A[i] = 2.0 * rand(s)
        end
    end
    return A
end

skewness(d::Chisq) = sqrt(8.0 / d.df)

var(d::Chisq) = 2.0 * d.df


text:
module XFails

using Test

struct XFailTestSet <: Test.AbstractTestSet
    ts::Test.AbstractTestSet

    XFailTestSet(ts::Test.AbstractTestSet) = new(ts)
end

function XFailTestSet(args...; kwargs...)
    ts = Test.get_testset()
    if ts isa XFailTestSet
        ts = ts.ts
    end
    TestSetType = typeof(ts)
    @assert !(TestSetType isa XFailTestSet)
    XFailTestSet(TestSetType(args...; kwargs...))
end

function Test.finish(ts::XFailTestSet)
    if Test.get_testset_depth() != 0
        parent_ts = Test.get_testset()
        if parent_ts isa XFailTestSet
            Test.record(parent_ts.ts, ts.ts)
            return ts.ts
        end
    end
    return Test.finish(ts.ts)
end

Test.record(ts::XFailTestSet, result::Test.Pass) = Test.record(
    ts.ts,
    Test.Fail(
        Symbol(:xfail_, result.test_type),
        result.orig_expr,
        result.data,
        result.value,
        LineNumberNode(0, Symbol("<unknown>")),
    ),
)
Test.record(ts::XFailTestSet, result::Test.Fail) = Test.record(
    ts.ts,
    Test.Pass(
        Symbol(:xfail_, result.test_type),
        result.orig_expr,
        result.data,
        result.value,
    ),
)
Test.record(ts::XFailTestSet, result::Test.Error) = Test.record(ts.ts, result)

end  # module


text:
# Run tests on Beyns contour integral method

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

push!(LOAD_PATH, string(@__DIR__, "/../src/gallery_extra"))
using GalleryNLEVP

@testset "GUN (NLEVP interface)" begin

    nep_org = nep_gallery(NLEVP_NEP,"gun")
    nep1 = nep_gallery("nlevp_native_gun")
    n=size(nep1,1);
    tol=1e-11;
    @bench @testset "Running alg" begin
        λ1,v1=quasinewton(nep1,λ=150^2+1im,v=ones(n),logger=displaylevel,tol=tol,maxit=500);
        λ1 = λ1[1]
        v1 = vec(v1)
        normalize!(v1)

        @test norm(compute_Mlincomb(nep1,λ1,v1))<tol*100
        @test norm(compute_Mder(nep1,λ1)*v1)<tol*100

        @test norm(compute_Mlincomb(nep_org,λ1,v1))<tol*100
        @test norm(compute_Mder(nep_org,λ1)*v1)<tol*100

    end
    @bench @testset "Compute derivatives" begin
        # Test compute_Mlincomb:
        λ=150^2+2im;
        V=randn(n,2);    v=V[:,1];
        z1=compute_Mlincomb(nep1,λ,V[:,1], [1.0], 1)
        z2=compute_Mlincomb(nep_org,λ,V[:,1], [1.0], 1)
        # Compare with divided difference
        ee=1e-4
        z3=(compute_Mlincomb(nep_org,λ+ee,v)-compute_Mlincomb(nep_org,λ-ee,v))/(2*ee);

        @test norm(z3-z1)<(ee^2)*1000  # Compare native and FD
        @test norm(z3-z2)<(ee^2)*1000 # Compare MATLAB and FD. This fails on NLEVP 3.0 due to a bug in gun.m in NLEVP
        @test norm(z1-z2)<sqrt(eps()) # Compare MATLAB and native.  This fails on NLEVP 3.0 due to a bug in gun.m in NLEVP
    end

    @bench @testset "Compare SPMF and Native in alg" begin
        # Check that two steps of quasinewton always gives the same result
        λ_org=0
        try
            quasinewton(nep_org,maxit=2,λ=150^2+1im,v=ones(n),logger=displaylevel)
        catch e
            λ_org=e.λ
        end

        λ1=0
        try
            quasinewton(nep1,maxit=2,λ=150^2+1im,v=ones(n),logger=displaylevel)
        catch e
            λ1=e.λ
        end

        @test abs(λ1-λ_org)<sqrt(eps())

    end

    @bench @testset "Compare MATLAB loaded vs MAT-loaded" begin
        nep2=nep_gallery("nlevp_native_gun");
        z=ones(n); λ=150^2;
        @test (norm(compute_Mlincomb(nep2,λ,z)-compute_Mlincomb(nep1,λ,z))/norm(compute_Mlincomb(nep1,λ,z)) < 1e-14)
    end

end


text:
using FillArrays

"""
    ADI.SDIAlgorithm <: ADI.ADIAlgorithm

Spectral differential imaging (SDI) algorithms. These work on 4-D SDI tensors. To use these algorithms, simply treat them like functions

```julia
(::SDIAlgorithm)(data::AbstractArray{T,4}, angles, scales; kwargs...)
(::SDIAlgorithm)(data::AbstractArray{T,4}, angles, data_ref, scales; kwargs...)
```

The data is expected to be laid out in `(nλ, nf, ny, nx)` format, so you may need to `permutedims` before processing the data. The `scales` correspond to the relative wavelength scales for each spectrum, and can be retrieved with `HCIToolbox.scale_list`.

# Algorithms

The current SDI implementations are
* [`SingleSDI`](@ref)
* [`DoubleSDI`](@ref)
* [`SliceSDI`](@ref)
"""
abstract type SDIAlgorithm <: ADIAlgorithm end

"""
    SingleSDI(alg)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in a single pass. This means that each channel will be scaled and then concatenated together, so an SDI tensor `(nλ, nf, y, x)` becomes a stack `(nλ * nf, y, x)` which is processed by the underlying `alg` like ADI data.

!!! tip
    `SingleSDI` is the default SDI mode. This means instead of writing
    ```julia
    SingleSDI(PCA(15))(data, angles, scales)
    ```
    you can just write
    ```julia
    PCA(15)(data, angles, scales)
    ```
"""
struct SingleSDI{ALG<:ADIAlgorithm} <: SDIAlgorithm
    alg::ALG
end

(alg::ADIAlgorithm)(spcube::AbstractArray{T,4}, angles, scales; kwargs...) where {T} =
    SingleSDI(alg)(spcube, angles, scales; kwargs...)

(alg::ADIAlgorithm)(spcube::AbstractArray{T,4}, angles, scales, spcube_ref; kwargs...) where {T} =
    SingleSDI(alg)(spcube, angles, scales, spcube_ref; kwargs...)

function (sdi::SingleSDI)(spcube::AbstractArray{T,4}, angles, scales; kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    big_cube = scale_and_stack(spcube, scales)

    # do single-pass reconstruction
    S = reconstruct(sdi.alg, big_cube, repeat(angles, inner=nλ); kwargs...)
    big_resid_cube = big_cube .- S
    # bin across spectral dim
    resid_cube = invscale_and_collapse(big_resid_cube, scales, frame_size)
    # derotate and combine
    return collapse!(resid_cube, angles; kwargs...)
end

function (sdi::SingleSDI)(spcube::AbstractArray{T,4}, angles, scales, spcube_ref; kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    big_cube = scale_and_stack(spcube, scales)
    big_cube_ref = scale_and_stack(spcube_ref, scales)

    # do single-pass reconstruction
    S = reconstruct(sdi.alg, big_cube, repeat(angles, inner=nλ), big_cube_ref; kwargs...)
    big_resid_cube = big_cube .- S
    # bin across spectral dim
    resid_cube = invscale_and_collapse(big_resid_cube, scales, frame_size)
    # derotate and combine
    return collapse!(resid_cube, angles; kwargs...)
end

"""
    DoubleSDI(alg)
    DoubleSDI(alg_spec, alg_temp)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in two passes. The first pass uses `alg_spec` to reduce each spectral cube slice in the SDI tensor. Then, the spectral residual frames will be reduced using `alg_temp`, which will include the derotation and final combination.

The difference between [`DoubleSDI`](@ref) and [`SliceSDI`](@ref) is that `DoubleSDI` does its first pass in the spectral slice, effectively collapsing the slice before performing ADI on the residual cube. `SliceSDI` does its first pass in the temporal slice, collapsing it first before performing ADI on the residual cube.

# Examples

```julia
julia> data, angles, scales = # load data...

# Median subtraction for each spectral slice,
# GreeDS{PCA} subtraction on spectral residual cube
julia> res = DoubleSDI(Median(), GreeDS(15))(data, angles, scales)
```
"""
struct DoubleSDI{ALG<:ADIAlgorithm, ALG2<:ADIAlgorithm} <: SDIAlgorithm
    alg_spec::ALG
    alg_temp::ALG2
end

DoubleSDI(alg) = DoubleSDI(alg, alg)

function (sdi::DoubleSDI)(spcube::AbstractArray{T,4}, angles, scales; kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    spec_resids = similar(spcube, n, ny, nx)
    # do first pass in spectral domain
    Threads.@threads for n in axes(spcube, 2)
        cube = @view spcube[:, n, :, :]
        scaled_cube = scale(cube, scales)
        S = reconstruct(sdi.alg_spec, scaled_cube, Fill(angles[n], nλ); kwargs...)
        spec_resid = invscale(scaled_cube .- S, scales, frame_size)
        spec_resids[n, :, :] .= collapse(spec_resid)
    end
    # do second pass in temporal domain
    return sdi.alg_temp(spec_resids, angles; kwargs...)
end

function (sdi::DoubleSDI)(spcube::AbstractArray{T,4}, angles, scales, spcube_ref; kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    spec_resids = similar(spcube, n, ny, nx)
    # do first pass in spectral domain
    Threads.@threads for n in axes(spcube, 2)
        cube = @view spcube[:, n, :, :]
        cube_ref = @view spcube_ref[:, n, :, :]
        scaled_cube = scale(cube, scales)
        scaled_cube_ref = scale(cube_ref, scales)
        S = reconstruct(sdi.alg_spec, scaled_cube, Fill(angles[n], nλ), scaled_cube_ref; kwargs...)
        spec_resid = invscale(scaled_cube .- S, scales, frame_size)
        spec_resids[n, :, :] .= collapse(spec_resid)
    end
    # do second pass in temporal domain
    return sdi.alg_temp(spec_resids, angles; kwargs...)
end

"""
    SliceSDI(alg)
    SliceSDI(alg_spec, alg_temp)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in two passes. The first pass uses `alg_temp` to reduce each temporal cube slice in the SDI tensor. These residuals will be rescaled and stacked into a new cube. Then, the temporal residual frames will be reduced using `alg_spec`, which will include the derotation and final combination.

The difference between [`SliceSDI`](@ref) and [`DoubleSDI`](@ref) is that `DoubleSDI` does its first pass in the spectral slice, effectively collapsing the slice before performing ADI on the residual cube. `SliceSDI` does its first pass in the temporal slice, collapsing it first before performing ADI on the residual cube.

# Examples

```julia
julia> data, angles, scales = # load data...

# Median subtraction for each spectral slice,
# GreeDS{PCA} subtraction on spectral residual cube
julia> res = SliceSDI(Median(), GreeDS(15))(data, angles, scales)
```
"""
struct SliceSDI{ALG<:ADIAlgorithm, ALG2<:ADIAlgorithm} <: SDIAlgorithm
    alg_spec::ALG
    alg_temp::ALG2
end

SliceSDI(alg) = SliceSDI(alg, alg)

function (sdi::SliceSDI)(spcube::AbstractArray{T,4}, angles, scales; kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    temp_resids = similar(spcube, nλ, ny, nx)
    # do first pass in temporal domain
    Threads.@threads for n in axes(spcube, 1)
        cube = spcube[n, :, :, :]
        temp_resids[n, :, :] .= sdi.alg_temp(cube, angles; kwargs...)
    end
    # do second pass in temporal domain
    scaled_resid_cube = scale(temp_resids, scales)
    resid = sdi.alg_spec(scaled_resid_cube, angles; kwargs...)
    return invscale(resid, maximum(scales))
end

function (sdi::SliceSDI)(spcube::AbstractArray{T,4}, angles, scales, spcube_ref; kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    temp_resids = similar(spcube, nλ, ny, nx)
    # do first pass in temporal domain
    Threads.@threads for n in axes(spcube, 1)
        cube = spcube[n, :, :, :]
        cube_ref = spcube_ref[n, :, :, :]
        temp_resids[n, :, :] .= sdi.alg_temp(cube, angles, cube_ref; kwargs...)
    end
    # do second pass in temporal domain
    scaled_resid_cube = scale(temp_resids, scales)
    resid = sdi.alg_spec(scaled_resid_cube, angles; kwargs...)
    return invscale(resid, maximum(scales))
end


text:
using Tables

function _construct_pandas_from_tables(source)
    if(!Tables.istable(source))
        return nothing
    end
    source_columns = Tables.columns(source)
    source_as_dict = Dict(column => Tables.getcolumn(source_columns, column) for column in Tables.columnnames(source_columns))
    return invoke(DataFrame, Tuple{Vararg{Any}}, source_as_dict)
end

Tables.columnaccess(::DataFrame) = true
Tables.rowaccess(::DataFrame) = true
Tables.istable(::DataFrame) = true


text:
import ProgressMeter
using Statistics: mean

mutable struct BasinInfo{B, IF, RF, UF, D, T, Q, K}
    basin::Array{Int16, B}
    grid_steps::SVector{B, Float64}
    grid_maxima::SVector{B, Float64}
    grid_minima::SVector{B, Float64}
    iter_f!::IF
    complete_and_reinit!::RF
    get_projected_state::UF
    state::Symbol
    current_att_color::Int
    current_bas_color::Int
    consecutive_match::Int
    consecutive_lost::Int
    prev_clr::Int
    attractors::Dict{Int16, Dataset{D, T}}
    visited::Q
    search_trees::K
end

"""
This is the low level function that creates & computes the basins of attraction, and is
agnostic of the dynamical system. `integ` is an integrator, `iter_f!` a function that
steps the integrator, `complete_and_reinit!` a function that re-inits the integrator
at a new full state, given the state on the grid.
"""
function draw_basin!(
        grid::Tuple, integ, iter_f!::Function, complete_and_reinit!, get_projected_state::Function;
        show_progress = true, attractors = nothing, kwargs...,
    )
    B = length(grid)
    D = length(get_state(integ)) # dimension of the full state space
    complete = false
    trees = if isnothing(attractors)
        nothing
    else
        Dict(k => searchstructure(KDTree, att, Euclidean()) for (k, att) in attractors)
    end
    grid_steps = step.(grid)
    grid_maxima = maximum.(grid)
    grid_minima = minimum.(grid)
    bsn_nfo = BasinInfo(
        ones(Int16, map(length, grid)),
        SVector{B, Float64}(grid_steps),
        SVector{B, Float64}(grid_maxima),
        SVector{B, Float64}(grid_minima),
        iter_f!,
        complete_and_reinit!,
        get_projected_state,
        :att_search,
        2,4,0,1,1,
        Dict{Int16,Dataset{D,eltype(get_state(integ))}}(),
        Vector{CartesianIndex{B}}(),
        trees
    )
    reset_basin_counters!(bsn_nfo)
    I = CartesianIndices(bsn_nfo.basin)
    j = 1
    progress = ProgressMeter.Progress(
        length(bsn_nfo.basin); desc = "Basins of attraction: ", dt = 1.0
    )

    while !complete
        ind, complete, j = next_uncolored_cell(bsn_nfo, j, I)
        show_progress && ProgressMeter.update!(progress, j)
        complete && break
        # Tentatively assign a color: odd is for basins, even for attractors.
        # First color is 2 for attractor and 3 for basins
        bsn_nfo.basin[ind] = bsn_nfo.current_bas_color
        y0 = generate_ic_on_grid(grid, ind)
        bsn_nfo.basin[ind] = get_color_point!(bsn_nfo, integ, y0; kwargs...)
    end
    # remove attractors and rescale from 1 to max number of attractors
    ind = iseven.(bsn_nfo.basin)
    bsn_nfo.basin[ind] .+= 1
    bsn_nfo.basin .= (bsn_nfo.basin .- 1) .÷ 2
    return bsn_nfo
end

function next_uncolored_cell(bsn_nfo, j, I)
    @inbounds for k in j:length(bsn_nfo.basin)
        if bsn_nfo.basin[I[k]] == 1
            j = k
            ind = I[k]
            return ind, false, j
        end
    end
    return I[1], true, length(bsn_nfo.basin)
end

function get_color_point!(bsn_nfo::BasinInfo, integ, y0; kwargs...)
    # This routine identifies the attractor using the previously defined basin.
    # reinitialize integrator
    bsn_nfo.complete_and_reinit!(integ, y0)
    reset_basin_counters!(bsn_nfo)
    cellcolor = 0

    while cellcolor == 0
        bsn_nfo.iter_f!(integ)
        new_y = bsn_nfo.get_projected_state(integ)
        n = basin_cell_index(new_y, bsn_nfo)
        u_att = get_state(integ) # in case we need the full state to save the attractor
        cellcolor = _identify_basin_of_cell!(bsn_nfo, n, u_att; kwargs...)
    end
    return cellcolor
end

"""
Main procedure motivated by Nusse & Yorke for the grid cell `n`. The algorithm can be
though as a finite state machine with five states: :att_hit, :att_search, :att_found,
:bas_hit, :lost. The transition between states depends on the current number in the
cell being visited by the trajectory of the dynamical systems. When the automata switches
to a state a counter starts and one of the following action can happen depending on the internal
counter and the input number:
* The automata stays in the same state and the counter increases.
* The automata switches to another state.
* The automata has decided upon the basin of the initial condition and returns a code corresponding
  to this basin.

The basins and attractors are coded in the array with odd numbers for the basins and even numbers
for the attractors. The attractor `2n` has the corresponding basin `2n+1`. This codification
is changed when the basins and attractors are returned to the user. Diverging trajectories
and the trajectories staying outside the grid are coded with -1.
"""
function _identify_basin_of_cell!(
        bsn_nfo::BasinInfo, n::CartesianIndex, u_full_state;
        mx_chk_att = 2, mx_chk_hit_bas = 10, mx_chk_fnd_att = 100, mx_chk_loc_att = 100,
        horizon_limit = 1e6, ε = 1e-3,
        mx_chk_lost = isnothing(bsn_nfo.search_trees) ? 20 : 1000,
    )

    #if n[1]==-1 means we are outside the grid
    nxt_clr = (n[1]==-1  || isnan(u_full_state[1])) ? -1 : bsn_nfo.basin[n]

    # search attractors directly
    if !isnothing(bsn_nfo.search_trees)
        for (k, t) in bsn_nfo.search_trees # this is a `Dict`
            idxs = isearch(t, u_full_state, WithinRange(ε))
            if !isempty(idxs)
                nxt_clr = 2*k + 1
                break
            end
        end
    end

    check_next_state!(bsn_nfo,nxt_clr)

    if bsn_nfo.state == :att_hit
        if nxt_clr == bsn_nfo.prev_clr
             bsn_nfo.consecutive_match += 1
        end
        if bsn_nfo.consecutive_match ≥ mx_chk_att
            # Wait if we hit the attractor a mx_chk_att times in a row just
            # to check if it is not a nearby trajectory
            hit_att = nxt_clr + 1
            recolor_visited_cell!(bsn_nfo, bsn_nfo.current_bas_color, 1)
            reset_basin_counters!(bsn_nfo)
            return hit_att
         end
         bsn_nfo.prev_clr = nxt_clr
         return 0
    end

    if bsn_nfo.state == :att_search
        if nxt_clr == 1
            # uncolored box, color it with current odd color and reset counter
            bsn_nfo.basin[n] = bsn_nfo.current_bas_color
            push!(bsn_nfo.visited,n) # keep track of visited cells
            bsn_nfo.consecutive_match = 1
        elseif nxt_clr == bsn_nfo.current_bas_color
            # hit a previously visited box with the current color, possible attractor?
            bsn_nfo.consecutive_match += 1
        end

        if bsn_nfo.consecutive_match >= mx_chk_fnd_att
            bsn_nfo.basin[n] = bsn_nfo.current_att_color
            store_attractor!(bsn_nfo, u_full_state)
            bsn_nfo.state = :att_found
            bsn_nfo.consecutive_match = 1
        end
        bsn_nfo.prev_clr = nxt_clr
        return 0
    end

    if bsn_nfo.state == :att_found
        if nxt_clr == 1 || nxt_clr == bsn_nfo.current_bas_color
            # Maybe chaotic attractor, perodic or long recursion.
            # Color this box as part of an attractor
            bsn_nfo.basin[n] = bsn_nfo.current_att_color
            bsn_nfo.consecutive_match = 1
            store_attractor!(bsn_nfo, u_full_state)
        elseif iseven(nxt_clr) && (bsn_nfo.consecutive_match <  mx_chk_loc_att)
            # We make sure we hit the attractor another mx_chk_loc_att consecutive times
            # just to be sure that we have the complete attractor
            bsn_nfo.consecutive_match += 1
        elseif iseven(nxt_clr) && bsn_nfo.consecutive_match >= mx_chk_loc_att
            # We have checked the presence of an attractor: tidy up everything
            # and get a new cell
            recolor_visited_cell!(bsn_nfo, bsn_nfo.current_bas_color, 1)
            # pick the next color for coloring the basin.
            bsn_nfo.current_bas_color += 2
            bsn_nfo.current_att_color += 2
            reset_basin_counters!(bsn_nfo)
            return nxt_clr + 1;
        end
        return 0
    end

    if bsn_nfo.state == :bas_hit
        # hit a colored basin point of the wrong basin, happens all the time,
        # we check if it happens mx_chk_hit_bas times in a row or if it happens
        # N times along the trajectory whether to decide if it is another basin.
        if bsn_nfo.prev_clr == nxt_clr
            bsn_nfo.consecutive_match += 1
        else
            bsn_nfo.consecutive_match = 1
        end
        if  bsn_nfo.consecutive_match > mx_chk_hit_bas
            recolor_visited_cell!(bsn_nfo, bsn_nfo.current_bas_color, 1)
            reset_basin_counters!(bsn_nfo)
            return nxt_clr
        end
        bsn_nfo.prev_clr = nxt_clr
        return 0
    end

    if bsn_nfo.state == :lost
        #grid_mid_point = (bsn_nfo.grid_maxima - bsn_nfo.grid_minima) ./2 + bsn_nfo.grid_minima
        bsn_nfo.consecutive_lost += 1
        if   bsn_nfo.consecutive_lost > mx_chk_lost || norm(u_full_state) > horizon_limit
            recolor_visited_cell!(bsn_nfo, bsn_nfo.current_bas_color, 1)
            reset_basin_counters!(bsn_nfo)
            # problematic IC : diverges or wanders outside the defined grid
            return -1
        end
        bsn_nfo.prev_clr = nxt_clr
        return 0
    end
end

function store_attractor!(bsn_nfo::BasinInfo{B, IF, RF, UF, D, T, Q},
    u_full_state) where {B, IF, RF, UF, D, T, Q}
    # bsn_nfo.current_att_color is the number of the attractor multiplied by two
    attractor_id = bsn_nfo.current_att_color ÷ 2
    V = SVector{D, T}
    if haskey(bsn_nfo.attractors, attractor_id)
        push!(bsn_nfo.attractors[attractor_id], V(u_full_state))
    else
        # initialize container for new attractor
        bsn_nfo.attractors[attractor_id] = Dataset([V(u_full_state)])
    end
end

function recolor_visited_cell!(bsn_nfo::BasinInfo, old_c, new_c)
    while !isempty(bsn_nfo.visited)
        ind = pop!(bsn_nfo.visited)
        if bsn_nfo.basin[ind] == old_c
            bsn_nfo.basin[ind] = new_c
        end
    end
end

function basin_cell_index(y_grid_state, bsn_nfo::BasinInfo{B}) where {B}
    iswithingrid = true
    @inbounds for i in 1:length(bsn_nfo.grid_minima)
        if !(bsn_nfo.grid_minima[i] ≤ y_grid_state[i] ≤ bsn_nfo.grid_maxima[i])
            iswithingrid = false
            break
        end
    end
    if iswithingrid
        # Snap point to grid
        ind = @. round(Int, (y_grid_state - bsn_nfo.grid_minima)/bsn_nfo.grid_steps) + 1
        return CartesianIndex{B}(ind...)
    else
        return CartesianIndex{B}(-1)
    end
end

function reset_basin_counters!(bsn_nfo::BasinInfo)
    bsn_nfo.consecutive_match = 0
    bsn_nfo.consecutive_lost = 0
    bsn_nfo.prev_clr = 1
    bsn_nfo.state = :att_search
end

function check_next_state!(bsn_nfo, nxt_clr)
    next_state = :undef
    current_state = bsn_nfo.state
    if current_state == :att_found
        # this is a terminal state, once reached you don't get out
        return
    end

    if nxt_clr == 1 || nxt_clr == bsn_nfo.current_bas_color
        # uncolored box or previously visited box with the current color
        next_state = :att_search
    elseif iseven(nxt_clr)
        # hit an attractor box
        next_state = :att_hit
    elseif nxt_clr == -1
        # out of the grid we do not reset the counter of other state
        # since the trajectory can follow an attractor that spans outside the grid
        bsn_nfo.state = :lost
        return
    elseif isodd(nxt_clr)
        # hit an basin box
        next_state = :bas_hit
    end

    if next_state != current_state
        # reset counter except in lost state (the counter freezes in this case)
        if current_state == :lost
            bsn_nfo.consecutive_lost = 1
        else
            bsn_nfo.consecutive_match = 1
        end
    end
    bsn_nfo.state = next_state
end


text:
module TableInterface
import Compose: compose, context, fill, rectangle, UnitBox
import Formatting: format, generate_formatter
import Base: contains, convert, done, getindex,
  intersect, length, next, size, start, writemime, ==

if Base.VERSION < v"0.4"
  using Docile
  import Base: get, isequal, eltype
  include("backport/nullable.jl")
  export Nullable, get, isequal, isnull, eltype
end

export
  done, getindex, length, lengths, next, size, start, ==,
  Box, Point,
  AbstractDataElement, ExpandDataElement, CartesianDataElement,
    isvertical,
  CellRange,
  Layout,
  AbstractStyleElement, AbstractPresentationStyleElement,
    CartesianStyleElement, ExpandStyleElement, UniformStyleElement,
    style, styles,
  Table,
  AbstractCell, AbstractCells, render, convert, writerange!, tcell

include("Box.jl")
include("DataElement.jl")
include("Range.jl")
include("Layout.jl")
include("StyleElement.jl")
include("Style.jl")
include("Table.jl")
include("render.jl")

end # module

module Tables
module HTMLFormat
  import Base: convert, writemime
  import Formatting: generate_formatter
  using TableInterface
  import TableInterface: style, styles, tcell, renderformatspecific

  if Base.VERSION < v"0.4"
    export Nullable, get, isequal, isnull, eltype
  end

  export HTMLTable, HTMLCell, HTMLStyleElement, style, styles, writemime
  include("format/converters.jl")
end
end


text:
using StructHelpers: @batteries, StructHelpers
const SH = StructHelpers
using Test

struct SVanilla
    a
    b
end
struct SBatteries
    a
    b
end
@batteries SBatteries
struct SNoHash
    a
    b
end
@batteries SNoHash hash=false
struct Skw
    a
    b
end
@batteries Skw kwconstructor=true kwshow=true

struct Empty1 end
struct Empty2 end
@batteries Empty1
@batteries Empty2

struct Salt1 end
struct Salt1b end
struct Salt2 end
@batteries Salt1 typesalt = 1
@batteries Salt1b typesalt = 1
@batteries Salt2 typesalt = 2
struct NoSalt end
@batteries NoSalt

struct SErrors;a;b;c;end

@testset "@batteries" begin
    @test SBatteries(1,2) == SBatteries(1,2)
    @test SBatteries(1,[]) == SBatteries(1,[])
    @test Skw(1,[]) == Skw(1,[])
    @test SNoHash(1,[]) == SNoHash(1,[])
    @test SVanilla(1,[]) != SVanilla(1,[])

    @test SBatteries(2,1) != SBatteries(1,2)
    @test Skw(2,1) != Skw(1,2)
    @test SNoHash(2,1) != SNoHash(1,2)

    @test hash(SBatteries(1,[])) == hash(SBatteries(1,[]))
    @test hash(SVanilla(1,[])) != hash(SVanilla(1,[]))
    @test hash(SNoHash(1,[])) != hash(SNoHash(1,[]))

    @test hash(SBatteries(2,[])) != hash(SBatteries(1,[]))
    @test hash(SNoHash(2,[])) != hash(SNoHash(1,[]))
    @test hash(Skw(2,[])) != hash(Skw(1,[]))

    @test SH.getproperties(SBatteries(1,2)) === (a=1, b=2)
    @test SH.setproperties(SBatteries(1,2), (a=10, b=20)) === SBatteries(10, 20)

    @test Skw(a=1, b=2) === Skw(1,2)
    @test_throws MethodError SBatteries(a=1, b=2)

    s = sprint(show, Skw(a=1, b=2))
    @test occursin("=", s)

    s = sprint(show, SBatteries(1,2))
    @test !occursin("=", s)

    @test Empty1() !== Empty2()
    @test Empty1() != Empty2()
    @test hash(Empty1()) != hash(Empty2())

    @test_throws Exception @macroexpand @batteries SErrors kwconstructor="true"
    @test_throws Exception @macroexpand @batteries SErrors nonsense=true
    @macroexpand @batteries SErrors kwconstructor=true

    @test hash(Salt1()) === hash(Salt1b())
    @test hash(Salt1()) != hash(NoSalt())
    @test hash(Salt1()) != hash(Salt2())
end


text:
include("../common.jl")

version = v"10.0"

# Collection of sources required to build attr
sources = [GitSource(repo, "ab5e12aa84fd95d7b0cd7b7b3f27e6ea5eaba05e")]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="LLVM_full_jll", version=v"10.0.1")),
    #Dependency(PackageSpec(name="libLLVM_jll", version=v"10.0.1"))
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"7")

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil


text:
using HOODESolver
using DifferentialEquations
using LinearAlgebra
using Plots
using Random

function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (0, 0)
    for i=1:size(tab,1)
        for j=(i+1):size(tab,1)
            nm = norm(tab[i]-tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret, nmmin
end

function fctmain(n_tau, prec)

    setprecision(prec)
    Random.seed!(5612)
    u0=rand(BigFloat,4)

    t_max = big"1.0"
    epsilon=big"1e-3"
    println("epsilon=$epsilon")
    nbmaxtest=10
    ordmax=12
    debord=3
    pasord=1
    y = ones(Float64, nbmaxtest, div(ordmax-debord,pasord)+1 )
    fct = u -> [ u[2]^2,1/(1+u[3]^2),0.5-u[4],0.8u[1] ]
    res_gen = Array{ Array{BigFloat,1}, 2}(undef, nbmaxtest, div(ordmax-debord,pasord)+1 )
    x=zeros(Float64,nbmaxtest)
    A = [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    parphi = PreparePhi(epsilon, n_tau, A, fct)
    nm = NaN
    ordmax += 1
    ordprep=undef
    nb = 20*2^(nbmaxtest)
    solref=undef
    while isnan(nm)
        ordmax -= 1
        ordprep = ordmax+2
    	@time par_u0 = PrepareU0(parphi, ordprep, u0)
        @time pargen = PrepareTwoScalePureAB(nb, t_max, ordmax, par_u0)
        @time solref= twoscales_pure_ab(
pargen,
only_end=true,
diff_fft=false
)
	    nm = norm(solref, Inf)
    end

    tabsol = Array{Array{BigFloat,1},1}(undef,1)

    tabsol[1] = solref

    println("ordmax=$ordmax")

    println("solref=$solref")
    solref_gen = solref
    indref = 1

    ind=1
    for order=debord:pasord:ordmax
        ordprep=order+2
        println("eps=$epsilon solRef=$solref order=$order")
        nb = 10
        indc = 1
        labels=Array{String,2}(undef, 1, order-debord+1)
        resnorm=0
        resnormprec=1
        sol =undef
        println("preparation ordre $order + 2")
        @time par_u0 = PrepareU0(parphi, ordprep, u0)
        while indc <= nbmaxtest
            @time pargen = PrepareTwoScalePureAB(nb, t_max, order, par_u0)
            @time sol = twoscales_pure_ab(
    pargen,
    only_end=true,
    diff_fft=false
)
            push!(tabsol, sol)
            res_gen[indc,ind] = sol
            (a, b), nm = getmindif(tabsol)
            if a != indref
                println("New solref !!!! a=$a, b=$b nm=$nm")
                indref = a
                solref = tabsol[a]
                for i=1:ind
                    borne = (i <ind) ? size(res_gen,1) : indc
                    for j = 1:borne
                        nm2 = min( norm(res_gen[j,i] - solref, Inf), 1.1)
                        y[j,i] = nm2 == 0 ? nm : nm2
                   end
                end
            else
                diff=solref-sol
                y[indc,ind] = min(norm(diff,Inf), 1.1)
            end
            x[indc] = 1.0/nb
            println("result=$y")
            println("log2(y)=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i=debord:pasord:order
            labels[1,(i-debord)÷pasord+1] = " eps,order=$(convert(Float32,epsilon)),$i "
        end
        p=Plots.plot(
                        x,
                        view(y,:,1:ind),
                        xlabel="delta t",
                        xaxis=:log,
                        ylabel="error",
                        yaxis=:log,
                        legend=:topleft,
                        label=labels,
                        marker=2
                    )
        prec_v = precision(BigFloat)
        eps_v = convert(Float32,epsilon)
        Plots.savefig(p,"out/r2_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_fct3.pdf")
        if resnorm > resnormprec
            break
        end
        resnormprec = resnorm
        ind+= 1
    end
end

# testODESolver()

# for i=3:9
#     fctMain(2^i)
# end
fctmain(256,512)


text:
# Test the getvalue function for each variable type
@test isa(getvalue(T), Vector)
@test isa(getvalue(Qc), Float64)
@test isa(getvalue(x), Float64)


text:
import FiniteDifferences

function ngradient(f, args...)
    fdm = FiniteDifferences.central_fdm(5, 1)
    FiniteDifferences.grad(fdm, f, args...)
end

function gradcheck(f, args...; atol=1e-5, rtol=1e-5)
    y_grads = grad(f, args...)[2]
    # don't check gradient w.r.t. function since ngradient can't do it
    y_grads = y_grads[2:end]
    n_grad = ngradient(f, args...)
    results = []
    for n in 1:length(args)
        if y_grads[n] isa NoTangent
            push!(results, true)
        else
            push!(results, isapprox(y_grads[n], n_grad[n], rtol=rtol, atol=atol))
        end
    end
    return all(results)
end


text:
function branch_basekv(network_data::Dict{String, Any})
    basekv = Dict{String, Float64}()
    for branch_idx in keys(network_data["branch"])
        br = network_data["branch"][branch_idx]
        f, t = br["f_bus"], br["t_bus"]
        bf, bt = network_data["bus"][string(f)], network_data["bus"][string(t)]
        basekv[branch_idx] = max(bf["base_kv"], bt["base_kv"])
    end
    return basekv
end


text:
using Test, LinearMaps, LinearAlgebra, Quaternions

@testset "noncommutative number type" begin
    x = Quaternion.(rand(10), rand(10), rand(10), rand(10))
    v = rand(10)
    A = Quaternion.(rand(10,10), rand(10,10), rand(10,10), rand(10,10))
    α = UniformScaling(Quaternion.(rand(4)...))
    β = UniformScaling(Quaternion.(rand(4)...))
    L = LinearMap(A)
    @test Array(L) == A
    @test Array(L') == A'
    @test Array(transpose(L)) == transpose(A)
    @test Array(α * L) == α * A
    @test Array(L * α) == A * α
    @test Array(α * L) == α * A
    @test Array(L * α ) == A * α
    @test Array(α * L') == α * A'
    @test Array((α * L')') ≈ (α * A')' ≈ A * conj(α)
    @test L * x ≈ A * x
    @test L' * x ≈ A' * x
    @test α * (L * x) ≈ α * (A * x)
    @test α * L * x ≈ α * A * x
    @test (α * L') * x ≈ (α * A') * x
    @test (α * L')' * x ≈ (α * A')' * x
    @test (α * L')' * v ≈ (α * A')' * v
    @test Array(@inferred adjoint(α * L * β)) ≈ conj(β) * A' * conj(α)
    @test Array(@inferred transpose(α * L * β)) ≈ β * transpose(A) * α
end

@testset "nonassociative number type" begin
    x = Octonion.(rand(10), rand(10), rand(10), rand(10),rand(10), rand(10), rand(10), rand(10))
    v = rand(10)
    A = Octonion.(rand(10,10), rand(10,10), rand(10,10), rand(10,10),rand(10,10), rand(10,10), rand(10,10), rand(10,10))
    α = UniformScaling(Octonion.(rand(8)...))
    β = UniformScaling(Octonion.(rand(8)...))
    L = LinearMap(A)
    @test Array(L) == A
    @test Array(α * L) == α * A
    @test Array(L * α) == A * α
    @test Array(α * L) == α * A
    @test Array(L * α) == A * α
    @test (α * L')' * v ≈ (α * A')' * v
end


text:
using Distributions
using GaussianMixtures

struct DataGenGMM
    dimensions::Vector{Int}
    num_gaussians::Vector{Int}
    num_repetitions::Int
    num_data::Int
    outlier_threshold::Float64
    num_gaussians_train::Vector{Int}
end

function DataGenGMM(dimensions::Vector{Int}, num_gaussians::Vector{Int}, num_repetitions::Int, num_data::Int, outlier_threshold::Float64)
    return DataGenGMM(dimensions, num_gaussians, num_repetitions, num_data, outlier_threshold, num_gaussians)
end

function rand_gmm(; num_gaussians=1, d=2)
    μ = rand(Normal(0.5, 0.1), num_gaussians, d)
    Σ = Array{eltype(GaussianMixtures.FullCov{Float64})}(undef, num_gaussians)
    for i=1:num_gaussians
        T = rand(Normal(0, 0.1), d, d)
        Σ[i] = GaussianMixtures.cholinv(T' * T / d)
    end
    w = ones(num_gaussians) / num_gaussians
    hist = History("custom")
    return GMM(w, μ, Σ, [hist], 0)
end

function generate_gmm_test_data(gmm_sample::GMM, gmm_test::GMM, threshold::Float64, n_inlier::Int, n_outlier::Int)
    range_check(x) = vec(all(x .>= 0, dims=1)) .& vec(all(x .<= 1, dims=1))
    dist = MixtureModel(gmm_test)
    data_inliers, data_outliers = [], []
    while length(data_inliers) < n_inlier
        candidates = copy(transpose(rand(gmm_sample, n_inlier * 2)))
        mask = range_check(candidates) .& (pdf(dist, candidates) .>= threshold)
        candidates = candidates[:, mask]
        i = 1
        while i <= size(candidates, 2) && length(data_inliers) < n_inlier
            push!(data_inliers, candidates[:, i])
            i += 1
        end
    end
    @info "Inlier generation done."
    while length(data_outliers) < n_outlier
        candidates = rand(gmm_sample.d, n_outlier * 2)
        mask = range_check(candidates) .& (pdf(dist, candidates) .< threshold)
        candidates = candidates[:, mask]
        i = 1
        while i <= size(candidates, 2) && length(data_outliers) < n_outlier
            push!(data_outliers, candidates[:, i])
            i += 1
        end
    end
    @info "Outlier generation done."
    data = hcat(data_inliers..., data_outliers...)
    labels = vcat(fill(:inlier, n_inlier), fill(:outlier, n_outlier))
    return data, labels
end
generate_gmm_test_data(gmm::GMM, threshold::Float64, n_inlier::Int, n_outlier::Int) = generate_gmm_test_data(gmm, gmm, threshold, n_inlier, n_outlier)
generate_gmm_test_data(gmm::GMM, threshold::Float64, n::Int) = generate_gmm_test_data(gmm, threshold, n, n)

function deconstruct_gmm(g::GMM, components::Vector{Int})
    all([x ∈ Set(1:length(g.w)) for x in components]) || throw(ArgumentError("GMM does not contain all components: $components."))
    return GMM(ones(length(components)) / length(components), g.μ[components, :], g.Σ[components], g.hist, g.nx)
end

function gmm_to_string(g)
    s = join([g.n, g.d], ',') * '\n'
    s *= join(g.w, ',') * '\n'
    for i in 1:g.n
        s *= join(g.μ[i, :], ',') * '\n'
        s *= join(vcat(transpose(GaussianMixtures.covar(g.Σ[i]))...), ',') * '\n'
    end
    return s
end

max_stars_repo_name: englhardt/des-evaluation

max_issues_repo_name: englhardt/des-evaluation

max_forks_repo_name: englhardt/des-evaluation


text:
export
    isperfectsquare,
    factor, eulerphi,
    divisorcount, divisorsigma,
    isperfect, isdeficient, isabundant,
    primefactors, factors,
    least_number_with_d_divisors

isperfectsquare(n::Integer) = n == isqrt(n)^2

Nemo.factor(n::Int128) = (factor ∘ big)(n)
Nemo.factor(n::T) where {T<:Integer} = factor(n, Dict{T, Int})
Nemo.factor(n::Integer, out::Type{O}) where {OP<:Integer, OK<:Integer, O<:Associative{OP, OK}} = begin
    factorization = (factor ∘ fmpz)(n)
    pairs = (convert(OP, p) => convert(OK, k) for (p, k) in factorization)
    out(pairs)
end
Nemo.factor(n::Integer, out::Type{Nemo.Fac}) = (factor ∘ fmpz)(n)

##  Euler's Phi (or: totient) function
Nemo.eulerphi(n::T) where {T<:Integer} = convert(T, (eulerphi ∘ fmpz)(n))
Nemo.eulerphi(f::Associative{T, <:Integer}) where {T<:Integer} = begin
    phi = one(T)
    for (p, k) in f
        phi *= p^(k-1) * (p - 1)
    end
    phi
end

# https://oeis.org/wiki/Divisor_function
divisorcount(n::T) where {T<:Integer} = convert(T, sigma(fmpz(n), 0))
divisorcount(f::Associative{T, <:Integer}) where {T<:Integer} = begin
    c = one(T)
    for k in values(f)
        c *= k + 1
    end
    c
end

divisorsigma(n::T, s = 1) where {T<:Integer} = convert(T, sigma(fmpz(n), s))
divisorsigma(f::Associative{T, <:Integer}, s = 1) where {T<:Integer} = begin
    s >= 0 || error("Argument 's' must be an integer greater or equal 0")
    s == 0 && return divisorcount(f)

    σ = one(T)
    for (p, k) in f
        σ *= (p^((k + 1) * s) - 1) ÷ (p^s - 1)
    end
    σ
end

isperfect(n::Integer) = divisorsigma(n, 1) - n == n
isdeficient(n::Integer) = divisorsigma(n, 1) - n < n
isabundant(n::Integer) = divisorsigma(n, 1) - n > n

primefactors(n::Integer) = (sort! ∘ collect ∘ keys ∘ factor)(n)

# http://rosettacode.org/wiki/Factors_of_an_integer
factors(n::T, negative::Bool = false) where {T<:Integer} = begin
    n > 0 || error("Argument 'n' must be an integer greater 0")

    f = [one(n)]
    for (p, k) in factor(n)
        f = reduce(vcat, f, [f * p^j for j in 1:k])::Array{T,1}
    end

    if length(f) == 1
        f = [one(n), n]
    else
        sort!(f)
    end

    negative ? flatten([f -f]') : f
end

indexfactorization2number(x::Array{T,1}) where {T<:Integer} =
    prod(big(nthprime(i))^k for (i, k) in enumerate(x))

# http://www.primepuzzles.net/problems/prob_019.htm
least_number_with_d_divisors(d::Integer) =
    minimum(indexfactorization2number(e) for e in least_number_with_d_divisors_exponents(d))

least_number_with_d_divisors_exponents(d::T, i::Int = 1, prevn::T = zero(T)) where {T<:Integer} = begin
    d <= 1 && return Vector{T}[T[]]

    f = factor(d, SortedDict{T, Int})
    pmax = last(f)[1]
    k = sum(values(f))

    p = nthprime(k+i-1)
    p_i = nthprime(i)
    m = floor(Integer, log(p) / log(p_i))

    c = [pmax]
    for b in 2:m
        !(b in keys(f)) && continue

        for a in b:m
            a*b <= pmax && continue
            d % (a*b) != 0 && continue
            first(primefactors(a)) < b && continue
            i > 1 && a*b > prevn && continue

            push!(c, a*b)
        end
    end

    ans = Vector{T}[]
    for ni in c
        for tailn in least_number_with_d_divisors_exponents(d ÷ ni, i+1, ni)
            push!(tailn, ni-1)
            push!(ans, tailn)
        end
    end

    i == 1 ? [reverse(x) for x in ans] : ans
end


text:
module RandomBattles

if isdefined(Base, :Experimental) && isdefined(Base.Experimental,
      Symbol("@optlevel"))
      @eval Base.Experimental.@optlevel 3
end

# grab data and Pokemon structs so that it can be used
include("data.jl")
include("stats.jl")
include("moves.jl")

# team battles
include("team/pokemon.jl")
include("team/team.jl")
include("team/state.jl")
include("team/mechanics.jl")
include("team/decisions.jl")
include("team/battle_logic.jl")
include("team/diff.jl")
include("team/strategy.jl")
include("team/nash.jl")
include("team/mcts.jl")
#include("team/alpha_beta_pruning.jl")

# individual battles
include("individual/state.jl")
include("individual/mechanics.jl")
include("individual/decisions.jl")
include("individual/battle_logic.jl")
#include("individual/diff.jl")
include("individual/strategy.jl")
include("individual/nash.jl")

# higher level abstractions
include("pokemon_meta.jl")
include("rankings.jl")

# export useful variables and functions
export DynamicState, StaticState, DynamicTeam, StaticTeam, DynamicPokemon,
      StaticPokemon, DynamicIndividualState, StaticIndividualState,
      get_possible_decisions, battle_score, battle_scores, play_turn,
      play_battle, gamemaster,
      diff, reflect,
      Strategy, IndividualStrategy, plot_strategy, rank,
      PokemonMeta,
      SM, solve_battle,
      MCTSNode, update_tree_MCTS!, select_decisions_MCTS

end # module


text:
include("style.jl")

# Uncomment the following to install the necessary packages for plotting.
# Pkg.add(["Conda", "PyCall"])
# using Conda
# Conda.add("cartopy")

using PyCall, PyPlot, ClimateBase
ccrs = pyimport("cartopy.crs")
LONLAT = ccrs.PlateCarree()
DEFPROJ = ccrs.Mollweide()
Time = ClimateBase.Ti

"""
    earthscatter(A::ClimArray, projection = ccrs.Mollweide(); kwargs...) → fig, ax, cb
Plot a `ClimArray` with space type `GaussianEqualArea` as a scatter plot with the color
of the points being the value of `A` at these points. This requires that `A` has only
one dimension, the coordinates.

Keyword values are propagated to `scatter` and keys of interest are `cmap, vmin, vmax`.
"""
function earthscatter(A::ClimArray, projection = DEFPROJ; kwargs...)
    coords = dims(A, Coord)
    @assert length(dims(A)) == 1 "Input ClimArray must only have one `Coord` dimension."
    eqlon = [l[1] for l in coords]
    eqlat = [l[2] for l in coords]
    fig = figure()
    fig.tight_layout()
    ax = subplot(111, projection=projection)
    ax, cb = earthscatter!(ax, eqlon, eqlat, A.data, projection; kwargs...)
    if A.name ≠ Symbol("")
        ax.set_title(string(A.name))
    end
    return fig, ax, cb
end

function earthscatter!(ax, A, projection = DEFPROJ; kwargs...)
    coords = dims(A, Coord)
    lon = [l[1] for l in coords]
    lat = [l[2] for l in coords]
    earthscatter!(ax, lon, lat, A, projection; kwargs...)
end

function earthscatter!(ax, lon, lat, A, projection = DEFPROJ;
    coast = true, vmin = minimum(A), vmax = maximum(A), levels = 11,
    ticks = nothing, kwargs...)
    lon = Array(lon)
    lat = Array(lat)
    lvls = range(vmin, vmax, length = levels)
    cmap = matplotlib.cm.get_cmap(get(kwargs, :cmap, :viridis), length(lvls)-1)
    norm = matplotlib.colors.Normalize(vmin=lvls[1], vmax=lvls[end])
    cax, kw = matplotlib.colorbar.make_axes(ax,location="right",pad=0.02,shrink=0.8)
    cb = matplotlib.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm)

    ax.scatter(lon, lat; transform = LONLAT, c=A, kwargs..., cmap = cmap, vmin = vmin, vmax = vmax)
    coast && ax.coastlines()
    if projection == LONLAT
        ax.set_xticks(-180:60:180, crs=LONLAT)
        ax.set_yticks(-90:30:90, crs=LONLAT)
    end
    ax.gridlines(alpha = 0.25)
    cb.set_ticks(isnothing(ticks) ? lvls[1:max(1, levels÷5):end] : ticks)
    return ax, cb
end

"""
    earthsurface(A::ClimArray, projection = ccrs.Mollweide(); kwargs...) → fig, ax, cb
Plot a `ClimArray` with space type `LonLatGrid` as a surface plot with the color
of the surface being the value of `A`. This requires that `A` has exactly two dimensions,
`Lon, Lat`.

Keyword values are `coast=true` to enable coastlines and also `cmap, vmin, vmax, levels`.
"""
function earthsurface(A, projection = DEFPROJ;
    coast = true, kwargs...)
    fig = figure()
    fig.tight_layout()
    ax = subplot(111, projection=projection)
    ax, cb = earthsurface!(ax, A; kwargs...)
    if projection == LONLAT
        ax.set_xticks(-180:60:180, crs=LONLAT)
        ax.set_yticks(-90:30:90, crs=LONLAT)
    end
    if A.name ≠ Symbol("")
        ax.set_title(string(A.name))
    end
    return fig, ax, cb
end

function earthsurface!(ax, A::ClimArray; kwargs...)
    ax.gridlines(alpha = 0.25)
    lon, lat, s = _getwrappeddata(A)
    earthsurface!(ax, lon, lat, s; kwargs...)
end

function earthsurface!(ax, lon, lat, A;
    vmin = minimum(A), vmax = maximum(A), levels = 11,
    coast = true, specific_contour = nothing, kwargs...)

    lon, lat, A = _getwrappeddata(copy(Array(lon)), copy(Array(lat)), A)
    lvls = range(vmin, vmax, length = levels)

    # TODO: Add colorbar as a separate function
    lvls = range(vmin, vmax, length = levels)
    cmap = matplotlib.cm.get_cmap(get(kwargs, :cmap, :viridis), length(lvls)-1)
    # cmap.set_under("k")
    norm = matplotlib.colors.Normalize(vmin=lvls[1], vmax=lvls[end])
    cax, kw = matplotlib.colorbar.make_axes(ax,location="right",pad=0.02,shrink=0.8)
    cb = matplotlib.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm)
    coast && ax.coastlines(; color = "gray")

    ax.contourf(Array(lon), Array(lat), clamp.(A', vmin, vmax), lvls, cmap = cmap,
                transform=LONLAT, vmin=vmin, vmax = vmax)

    if specific_contour isa Array
        ax.contour(Array(lon), Array(lat), clamp.(A', vmin, vmax);
        levels = specific_contour, colors = "r", transform=LONLAT)
    end
    return ax, cb
end

function _getwrappeddata(A::ClimArray)
    lon = copy(Array(dims(A, Lon).val))
    lat = copy(Array(dims(A, Lat).val))
    s = copy(Array(A.data)) # assumes data have Lon x Lat dimensions
    _getwrappeddata(lon, lat, s)
end

function _getwrappeddata(lon, lat, s)
    if lon[end] == 359.5 && lon[1] == 0.5
        # add one more entry to have perfect overlap in plotting
        x = s[1, :]
        s = vcat(Array(s), x')
        push!(lon, 360.5)
    end
    return lon, lat, s
end

_locator = matplotlib.dates.AutoDateLocator(interval_multiples=true)
_locator.intervald
dateformater = matplotlib.dates.ConciseDateFormatter(
    matplotlib.dates.AutoDateLocator(interval_multiples=true)
)

PyPlot.plot(A::ClimArray{T, 1}; kwargs...) where {T} = plot!(gca(), A; kwargs...)
function plot!(ax, A::ClimArray{T, 1}; kwargs...) where {T}
    r = ax.plot(dims(A, 1).val, A.data; kwargs...)
    ax.set_xlabel(_xlabel(A))
    return r
end
function plot!(ax, A::ClimArray{T, 1, <: Tuple{<:Time}}; addlabels = true, kwargs...) where {T}
    t = dims(A, 1)
    x = 1:length(t.val)
    la = ["$(month(a))/$(year(a))" for a in t]
    r = ax.plot(x, A.data; kwargs...)
    ax.set_xlabel(_xlabel(A))
    i = [1, length(la)÷3, 2length(la)÷3, length(la)]
    ax.set_xticks(x[i])
    ax.set_xticklabels(la[i])
    return r
end

function plot!(ax, x, y; kwargs...) where {T}
    r = ax.plot(x, y; kwargs...)
    ax.set_xlabel(_xlabel(x))
end

function heatmap!(ax, A::ClimArray{T, 2}; kwargs...) where {T}
    x = dims(A, 2).val
    y = dims(A, 1).val
    a = A.data#[1:end-1, 1:end-1]
    m = ax.pcolormesh(x, y, a; kwargs...)
    ax.set_xlabel(_xlabel(dims(A)[2]))
    ax.set_ylabel(_xlabel(dims(A)[1]))
    return m
end

_xlabel(x) = ""
_xlabel(x::ClimArray) = _xlabel(dims(x)[1])
_xlabel(::Lon) = "longitude λ [ᵒ]"
_xlabel(::Lat) = "latitude φ [ᵒ]"
_xlabel(::Time) = "time [months]"

LANDMARKS_LON = Dict(
    350 => "West Africa",
    15 => "Europe",
    87 => "Himalaya",
    43 => "Saudi Ar.",
    76 => "India",
    135 => "Australia",
    282 => "Andes",
    300 => "Amazon",
    318 => "Greenland",
    334 => "Atlantic",
    209 => "Alaska",
    165 => "Pacific (W)",
    246 => "Pacific (E)",
    260 => "USA",
)

LANDMARKS_LAT = Dict(
    24 => "North\nAfrica",
    -25 => "Australia\nAndes",
    7 => "ITCZ",
    -5 => "Amazon\nCentral Africa",
    48 => "Europe\nMongolia",
    35 => "Himalaya",
    -2 => "Indian oc.",
    62 => "Canada\nScandinavia",
    -75 => "Antartica",
)

LANDMARKS_LAT_ABS = Dict(abs.(keys(LANDMARKS_LAT)) .=> values(LANDMARKS_LAT))

function landmarks!(ax, lmars)
    i = 0
    for (l, s) in lmars
        ax.axvline(l, ls = "dashed", lw = 1.0, alpha = 0.25, color = "k")
        ax.text(l, 0.98 - (i%3)*0.02, s, transform = ax.get_xaxis_transform(),
        size = 10, va="top", ha = "center")
        i+=1
    end
end

function plot_meanstd!(ax, y; color = "C0", kwargs...)
    m = mean(y); s = std(y)
    ax.axhline(m, lw = 1.5, color = color, ls = "dotted")
    ax.axhspan(m-s, m+s, alpha = 0.1, color = color,
        label = "mean±std: $(rdspl(m)) ± $(rdspl(s))"
    )
    for v in (m-s, m+s)
        ax.axhline(v, lw = 1.0, alpha = 0.5, ls = "solid",
        color = color)
    end
    ax.legend(loc = "lower left")
    return m, s
end

function plottrend!(ax, x, y; label="s", addlabel=true, kwargs...)
    if eltype(x) <: TimeType
        s, o, f, xf = linreg(1:length(x), y)
    else
        s, o, f, xf = linreg(x, y)
    end
    if addlabel
        ax.plot(x, f; label = label*"=$(round(s;sigdigits=3))", kwargs...)
        ax.legend()
    else
        ax.plot(x, f; kwargs...)
    end
    return s, o, f
end

function latitudinal_axis!(ax, llats = -70:20:70)
    ax.set_xticks(sind.(llats))
    ax.set_xticklabels(llats)
    ax.set_xlim(-1,1)
end

#########################################################################
# Color coding
#########################################################################
FIELDCOLORS = Dict(
    "toa_sw_all_mon" => ("C0", "viridis"),
    "toa_sw_clr_t_mon" => ("C3", "YlOrBr_r"),
    "toa_sw_clr_c_mon" => ("C3", "YlOrBr_r"),
    "solar_mon" => ("C5", "plasma"),
    "cldarea_total_daynight_mon" => ("C2", "Blues_r"),
    "sfc_sw_up_all_mon" => ("C4", "RdPu_r"),
    "sfc_sw_down_clr_c_mon" => ("C4", "PuRd_r"),
    "toa_cre_sw_mon" => ("C2", "bone"),
    "toa_sw_all_mon_per" => ("C2", "Blues_r"),
    "seasonal" => ("C2", "Blues_r"),
    "toa_sw_all_mon_res" => ("C4", "plasma"),
    "residual" => ("C4", "plasma"),
    "original" => ("C0", "viridis"),
)

fieldcolors(f) = get(FIELDCOLORS, string(f), ("C0", "binary_r"))

#########################################################################
# Color coding
#########################################################################
function crosscorrelation_plot(x, y, lags; ax = (figure(); gca()))
    cc = crosscor(x, y, lags)
    ax.plot(lags, cc, color = "C1", label = "CrossCor", ls = "dashed")
    ax.set_ylabel("crosscorrelation")
    mi, center, σ = mutualinformation(Array(x), Array(y), lags)
    ax.axhspan(center-3σ, center+3σ, label = "95% null", color = "C1", alpha = 0.25, zorder = -9)
    ax.axhline(center, lw = 1.0, ls = "dashed", color = "C1")
    ax.plot(lags, mi, label = "MI", color = c)
    ax.set_xlabel("delay [months]")
    ax.grid("on")
    ax.legend()
    return ax
end

function hemispheric_component_plot(C)
    Cn, Cs = hemispheric_means(zonalaverage(C))
    hemispheric_component_plot(Cn, Cs)
end

function hemispheric_component_plot(Cn, Cs)
    fig, ax = subplots()
    hemispheric_component_plot!(ax, Cn, Cs)
    return fig, ax
end

function hemispheric_component_plot!(ax, Cn, Cs)
    Cg = 0.5(Cn .+ Cs)
    plot!(ax, Cn, label = "NH, μ=$(rdspl(timeaverage(Cn)))")
    plot!(ax, Cs, label = "SH, μ=$(rdspl(timeaverage(Cs)))", ls = "--")
    plot!(ax, Cg, label = "GL, μ=$(rdspl(timeaverage(Cg)))", ls = ":")
    ax.legend()
    return ax
end

max_stars_repo_name: Datseris/EarthAlbedoSymmetry

max_issues_repo_name: Datseris/EarthAlbedoSymmetry

max_forks_repo_name: Datseris/EarthAlbedoSymmetry


text:
#=
# start of assignment 5 code
"""
    complement(base)

Get the DNA complement of the provided base:

    A <-> T
    G <-> C

Accepts uppercase or lowercase `String` or `Char`,
but always returns an uppercase `Char`.
If a valid base is not provided, the function throws an error.
"""
function complement(base)
    complements = Dict("A" => 'T',
                       "T" => 'A',
                       "G" => 'C',
                       "C" => 'G')

    base = uppercase(string(base))

    !(base in keys(complements)) && error("Invalid base $base")
    return complements[base]
end

"""
    isreverse(word1, word2)

    Checks to see if word1 is the reverse of word2
        Returns true if so, false otherwise
"""
function isreverse(word1, word2)
    if length(word1) != length(word2)
        return false
    end
    i = firstindex(word1)
    j = lastindex(word2)
    while j > 0

        #@show i j
        if word1[i] != word2[j]
            return false
        end
        i = nextind(word1, i)
        j = prevind(word2, j)
    end
    true
end
"""
    isreversecomplement(seq1, seq2)

Boolean function that checks whether two DNA seqences
are the reverse complement of one another, irrespective of capitalization.
Returns true if yes, false otherwise.

If any invalid bases are encountered,
or if sequences are different length, throws an error.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> isreversecomplement("aaatttcg", "cgaaattt")
    true

    julia> if isreversecomplement("C", "A")
               println("Yes!")
           else
               println("No!")
           end
    No!

    julia> isreversecomplement("TX", "AG")
    Error: Invalid base X

    julia> isreversecomplement("G", "CC")
    Error: Cannot compare sequuences of different length
"""
function isreversecomplement(seq1, seq2)
    ## need to check lengths, then reverse, as reverse if wrong base found then error
    seq2 = uppercase(seq2)
    if length(seq1) != length(seq2)
        throw(ErrorException("Cannot compare sequences of different length"))
    end

    #first find comp of 1, then find out if that is reverse of 2
    seq1comp = ""
    for letter in seq1
        seq1comp = seq1comp * complement(letter)
    end

    return isreverse(seq1comp, seq2)
end

"""
    reverse_complement(sequence)

Takes a DNA sequence and returns the reverse complement
of that sequence.

Takes lowercase or uppercase sequences,
but always returns uppercase.

Examples
≡≡≡≡≡≡≡≡≡≡
    julia> reverse_complement("AAATTT")
    "AAATTT"

    julia> reverse_complement("GCAT")
    "ATGC"

    julia> rc = reverse_complement("TTGGG");

    julia> println(rc)
    CCCAA
"""
function reverse_complement(sequence)
    ## your code here
    reverseseq = ""
    for letter in sequence
        reverseseq = reverseseq * complement(letter)
    end

    reverseseq = reverse(reverseseq)
end

"""
    isDNA(sequence)

checks the given sequence, and if one of its letters is not a base as stored in the bases array, then returns false.
otherwise, returns true
"""
function isDNA(sequence)
    bases = ['A', 'C', 'G', 'T']
    sequence = uppercase(sequence)
    for letter in sequence
        if !(letter ∈ bases)
            return false
        end
    end
    return true
end
# end of assignment 5 code

"""
   gc_content(sequence)

   Calculates the GC ratio of a DNA sequence.
   The GC ratio is the total number of G and C bases divided by the total length of the sequence.
   For more info about GC content, see here: https://en.wikipedia.org/wiki/GC-content

   Example
   ≡≡≡≡≡≡≡≡≡≡

       julia> question3("AATG")
       0.25

       julia> question3("CCCGG")
       1.0

       julia> question3("ATTA")
       0.0
"""
function gc_content(sequence)
    sequence = uppercase(sequence)

    seqlength = length(sequence)

    ## count the number of G's
    gs = count(==('G'), sequence)
    ## count the number of C's
    cs = count(==('C'), sequence)

    return (gs + cs) / seqlength
end
=#


text: if 1 -2 3 end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using HypothesisTests, Distributions, Base.Test

# sample drawn from uniform distribution
x = [0.3500, 0.1966, 0.2511, 0.6160, 0.4733,
     0.3517, 0.8308, 0.5853, 0.5497, 0.9172,
     0.2858, 0.7572, 0.7537, 0.3804, 0.5678,
     0.0759, 0.0540, 0.5308, 0.7792, 0.9340,
     0.1299, 0.5688, 0.4694, 0.0119, 0.3371
]
t = ApproximateOneSampleKSTest(x, Uniform())
@test t.δ ≈ 0.1440
@test t.δn ≈ 0.0571
@test t.δp ≈ 0.1440
@test pvalue(t) ≈ 0.6777349664784745
@test pvalue(t; tail=:left) ≈ 0.849573771973747
@test pvalue(t; tail=:right) ≈ 0.3545875485608989
show(IOBuffer(), t)

t = ApproximateTwoSampleKSTest(x, [(0:24)/25...])
@test t.δ ≈ 0.12
@test t.δn ≈ 0.08
@test t.δp ≈ 0.12
@test pvalue(t) ≈ 0.993764859699076
@test pvalue(t; tail=:left) ≈ 0.8521437889662113
@test pvalue(t; tail=:right) ≈ 0.697676326071031
show(IOBuffer(), t)

t = ExactOneSampleKSTest(x, Uniform())
@test t.δ ≈ 0.1440
@test t.δn ≈ 0.0571
@test t.δp ≈ 0.1440
@test pvalue(t) ≈ 0.6263437768244742
@test pvalue(t; tail=:left) ≈ 0.8195705417998183
@test pvalue(t; tail=:right) ≈ 0.32350648882777194
show(IOBuffer(), t)

## check fit to normal distribution
t = ApproximateOneSampleKSTest(x, Normal())
@test t.δ ≈ 0.5047473010922947
@test t.δn ≈ 0.5047473010922947
@test t.δp ≈ 0.17515194649718513
@test pvalue(t) ≈ 5.871827067532435e-6
@test pvalue(t; tail=:left) ≈ 2.9359135337662175e-6
@test pvalue(t; tail=:right) ≈ 0.21569061887162347

## check unequal sample size
t = ApproximateTwoSampleKSTest(x, [(0:5)/6...])
@test t.δ ≈ 0.22
@test t.δn ≈ 0.22
@test t.δp ≈ 0.09333333333333346
@test pvalue(t) ≈ 0.973300892518972
@test pvalue(t; tail=:left) ≈ 0.6260111498528065
@test pvalue(t; tail=:right) ≈ 0.9191544797498837

# http://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-fall-2006/lecture-notes/lecture14.pdf
x = [0.58, 0.42, 0.52, 0.33, 0.43, 0.23, 0.58, 0.76, 0.53, 0.64]
t = ApproximateOneSampleKSTest(x, Uniform())
@test t.δ ≈ 0.26
@test t.δn ≈ 0.23
@test t.δp ≈ 0.26
@test pvalue(t) ≈ 0.5084937988981307
@test pvalue(t; tail=:left) ≈ 0.3471494153245104
@test pvalue(t; tail=:right) ≈ 0.25872229825964005

t = ApproximateTwoSampleKSTest(x, [(0:9)/10...])
@test t.δ ≈ 0.3
@test t.δn ≈ 0.3
@test t.δp ≈ 0.2
@test pvalue(t) ≈ 0.7590978384203948
@test pvalue(t; tail=:left) ≈ 0.406569659740599
@test pvalue(t; tail=:right) ≈ 0.6703200460356393

t = ExactOneSampleKSTest(x, Uniform())
@test t.δ ≈ 0.26
@test t.δn ≈ 0.23
@test t.δp ≈ 0.26
@test pvalue(t) ≈ 0.4351284228580825
@test pvalue(t; tail=:left) ≈ 0.3013310572470338
@test pvalue(t; tail=:right) ≈ 0.2193143479950862


text:
# ---
# title: 9. Palindrome Number
# id: problem9
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/palindrome-number/description/>
# hidden: true
# ---
#
# Determine whether an integer is a palindrome. An integer is a palindrome when
# it reads the same backward as forward.
#
# **Follow up:** Could you solve it without converting the integer to a string?
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 121
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: x = -121
#     Output: false
#     Explanation: From left to right, it reads -121. From right to left, it becomes 121-. Therefore it is not a palindrome.
#
#
# **Example 3:**
#
#
#
#     Input: x = 10
#     Output: false
#     Explanation: Reads 01 from right to left. Therefore it is not a palindrome.
#
#
# **Example 4:**
#
#
#
#     Input: x = -101
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `-231 <= x <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function is_palindrome(x::Int)::Bool
    if x < 0
        false
    else
        y = x
        y_reverse = 0
        while y > 0
            y, n = divrem(y, 10)
            y_reverse = y_reverse * 10 + n
        end
        x == y_reverse
    end
end
## @lc code=end


text:
using Memoize, Base.Test

run = 0
@memoize function simple(a)
	global run += 1
	a
end
@test simple(5) == 5
@test run == 1
@test simple(5) == 5
@test run == 1
@test simple(6) == 6
@test run == 2
@test simple(6) == 6
@test run == 2

run = 0
@memoize function typed(a::Int)
	global run += 1
	a
end
@test typed(5) == 5
@test run == 1
@test typed(5) == 5
@test run == 1
@test typed(6) == 6
@test run == 2
@test typed(6) == 6
@test run == 2

run = 0
@memoize function default(a=2)
	global run += 1
	a
end
@test default() == 2
@test run == 1
@test default() == 2
@test run == 1
@test default(2) == 2
@test run == 1
@test default(6) == 6
@test run == 2
@test default(6) == 6
@test run == 2

run = 0
@memoize function default_typed(a::Int=2)
	global run += 1
	a
end
@test default_typed() == 2
@test run == 1
@test default_typed() == 2
@test run == 1
@test default_typed(2) == 2
@test run == 1
@test default_typed(6) == 6
@test run == 2
@test default_typed(6) == 6
@test run == 2

run = 0
@memoize function kw(; a=2)
	global run += 1
	a
end
@test kw() == 2
@test run == 1
@test kw() == 2
@test run == 1
@test kw(a=2) == 2
@test run == 1
@test kw(a=6) == 6
@test run == 2
@test kw(a=6) == 6
@test run == 2

run = 0
@memoize function kw_typed(; a::Int=2)
	global run += 1
	a
end
@test kw_typed() == 2
@test run == 1
@test kw_typed() == 2
@test run == 1
@test kw_typed(a=2) == 2
@test run == 1
@test kw_typed(a=6) == 6
@test run == 2
@test kw_typed(a=6) == 6
@test run == 2

run = 0
@memoize function default_kw(a=1; b=2)
	global run += 1
	(a, b)
end
@test default_kw() == (1, 2)
@test run == 1
@test default_kw() == (1, 2)
@test run == 1
@test default_kw(1, b=2) == (1, 2)
@test run == 1
@test default_kw(1, b=3) == (1, 3)
@test run == 2
@test default_kw(1, b=3) == (1, 3)
@test run == 2
@test default_kw(2, b=3) == (2, 3)
@test run == 3
@test default_kw(2, b=3) == (2, 3)
@test run == 3

run = 0
@memoize function default_kw_typed(a::Int=1; b::Int=2)
	global run += 1
	(a, b)
end
@test default_kw_typed() == (1, 2)
@test run == 1
@test default_kw_typed() == (1, 2)
@test run == 1
@test default_kw_typed(1, b=2) == (1, 2)
@test run == 1
@test default_kw_typed(1, b=3) == (1, 3)
@test run == 2
@test default_kw_typed(1, b=3) == (1, 3)
@test run == 2
@test default_kw_typed(2, b=3) == (2, 3)
@test run == 3
@test default_kw_typed(2, b=3) == (2, 3)
@test run == 3

run = 0
@memoize function ellipsis(a, b...)
	global run += 1
	(a, b)
end
@test ellipsis(1) == (1, ())
@test run == 1
@test ellipsis(1) == (1, ())
@test run == 1
@test ellipsis(1, 2, 3) == (1, (2, 3))
@test run == 2
@test ellipsis(1, 2, 3) == (1, (2, 3))
@test run == 2

run = 0
@memoize Dict function kw_ellipsis(;a...)
	global run += 1
	a
end
@test kw_ellipsis() == []
@test run == 1
@test kw_ellipsis() == []
@test run == 1
@test kw_ellipsis(a=1) == Any[(:a, 1)]
@test run == 2
@test kw_ellipsis(a=1) == Any[(:a, 1)]
@test run == 2
@test kw_ellipsis(a=1, b=2) == Any[(:a, 1), (:b, 2)]
@test run == 3
@test kw_ellipsis(a=1, b=2) == Any[(:a, 1), (:b, 2)]
@test run == 3

run = 0
@memoize function multiple_dispatch(a::Int)
	global run += 1
	1
end
@memoize function multiple_dispatch(a::Float64)
	global run += 1
	2
end
@test multiple_dispatch(1) == 1
@test run == 1
@test multiple_dispatch(1) == 1
@test run == 1
@test multiple_dispatch(1.0) == 2
@test run == 2
@test multiple_dispatch(1.0) == 2
@test run == 2

function outer()
	run = 0
	@memoize function inner(x)
		run += 1
		x
	end
	@memoize function inner(x, y)
		run += 1
		x+y
	end
	@test inner(5) == 5
	@test run == 1
	@test inner(5) == 5
	@test run == 1
	@test inner(6) == 6
	@test run == 2
	@test inner(6) == 6
	@test run == 2
	@test inner(5, 1) == 6
	@test run == 3
	@test inner(5, 1) == 6
	@test run == 3
end
outer()
@test !isdefined(:inner)


text:
@doc doc"""
    GraphManifoldType

This type represents the type of data on the graph that the [`GraphManifold`](@ref)
represents.
"""
abstract type GraphManifoldType end

@doc doc"""
    EdgeManifoldManifold <: GraphManifoldType

A type for a [`GraphManifold`](@ref) where the data is given on the edges.
"""
struct EdgeManifold <: GraphManifoldType end

@doc doc"""
    VectexGraphManifold <: GraphManifoldType

A type for a [`GraphManifold`](@ref) where the data is given on the vertices.
"""
struct VertexManifold <: GraphManifoldType end

@doc doc"""
    GraphManifold{G, M, T} <: AbstractPowerManifold{M,NestedPowerRepresentation}

Build a manifold, that is a [`PowerManifold`](@ref) of the [`Manifold`](@ref) `M` either on the edges or vertices
of a graph `G` depending on the [`GraphManifoldType`](@ref) `T`.

# Fields
* `G` is an `AbstractSimpleGraph`
* `M` is a [`Manifold`](@ref)
"""
struct GraphManifold{G<:AbstractGraph,TM,T<:GraphManifoldType} <: AbstractPowerManifold{TM,NestedPowerRepresentation}
    graph::G
    manifold::TM
end

function GraphManifold( g::G, M::TM, ::VertexManifold) where {G<:AbstractGraph,TM<:Manifold}
    return GraphManifold{G,TM,VertexManifold}(g, M)
end
function GraphManifold(g::G, M::TM, ::EdgeManifold) where {G<:AbstractGraph,TM<:Manifold}
    return GraphManifold{G,TM,EdgeManifold}(g, M)
end

const EdgeGraphManifold = GraphManifold{<:AbstractGraph,<:Manifold,EdgeManifold}
const VertexGraphManifold = GraphManifold{<:AbstractGraph,<:Manifold,VertexManifold}

@doc doc"""
    check_manifold_point(M::GraphManifold, x)

Check whether `x` is a valid point on the [`GraphManifold`](@ref), i.e. its
length equals the number of vertices (for [`VertexManifold`](@ref)s) or
the number of edges (for [`EdgeManifold`](@ref)s) and that each element of `x`
passes the [`check_manifold_point`](@ref) test for the base manifold `M.manifold`.
"""
check_manifold_point(::GraphManifold, ::Any...)
function check_manifold_point(M::VertexGraphManifold, x; kwargs...)
    if size(x) != (nv(M.graph),)
        return DomainError(
            length(x),
            "The number of points in `x` ($(length(x))) does not match the number of nodes in the graph ($(nv(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), nv(M.graph))
    return check_manifold_point(PM, x; kwargs...)
end
function check_manifold_point(M::EdgeGraphManifold, x; kwargs...)
    if size(x) != (ne(M.graph),)
        return DomainError(
            length(x),
            "The number of points in `x` ($(size(x))) does not match the number of edges in the graph ($(ne(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), ne(M.graph))
    return check_manifold_point(PM, x; kwargs...)
end

@doc doc"""
    check_tangent_vector(M::GraphManifold, x, v)

Check whether `x` is a valid point on the [`GraphManifold`](@ref), and
`v` it from its tangent space, i.e. its
length equals the number of vertices (for [`VertexManifold`](@ref)s) or
the number of edges (for [`EdgeManifold`](@ref)s) and that each element of `v`
together with its corresponding einty of `x` passes the
[`check_tangent_vector`](@ref) test for the base manifold `M.manifold`.
"""
check_tangent_vector(::GraphManifold, ::Any...)
function check_tangent_vector(M::VertexGraphManifold, x, v; kwargs...)
    if size(x) != (nv(M.graph),)
        return DomainError(
            length(x),
            "The number of points in `x` ($(size(x)) does not match the number of nodes in the graph ($(nv(M.graph))).",
        )
    end
    if size(v) != (nv(M.graph),)
        return DomainError(
            length(v),
            "The number of points in `v` ($(size(v)) does not match the number of nodes in the graph ($(nv(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), nv(M.graph))
    return check_tangent_vector(PM, x, v; kwargs...)
end
function check_tangent_vector(M::EdgeGraphManifold, x, v; kwargs...)
    if size(x) != (ne(M.graph),)
        return DomainError(
            length(x),
            "The number of elements in `x` ($(size(x)) does not match the number of edges in the graph ($(ne(M.graph))).",
        )
    end
    if size(v) != (ne(M.graph),)
        return DomainError(
            length(v),
            "The number of elements in `v` ($(size(v)) does not match the number of edges in the graph ($(ne(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), ne(M.graph))
    return check_tangent_vector(PM, x, v; kwargs...)
end

get_iterator(M::EdgeGraphManifold) = 1:ne(M.graph)
get_iterator(M::VertexGraphManifold) = 1:nv(M.graph)

@doc doc"""
    incident_log(M::GraphManifold, x)

Return the tangent vector on the (vertex) [`GraphManifold`](@ref), where at
each node the sum of the [`log`](@ref)s to incident nodes is computed.
For a `SimpleGraph`, an egde is interpreted as double edge in the corresponding
SimpleDiGraph

If the internal graph is a `SimpleWeightedGraph` the weighted sum of the
tangent vectors is computed.
"""
function incident_log(M::VertexGraphManifold, x)
    v = zero_tangent_vector(M, x)
    return incident_log!(M, v, x)
end

function incident_log!(M::VertexGraphManifold, v, x)
    rep_size = representation_size(M.manifold)
    for e in edges(M.graph)
        vw = _write(M, rep_size, v, src(e))
        v[src(e)] +=
            log(M.manifold, _read(M, rep_size, x, src(e)), _read(M, rep_size, x, dst(e)))
        if !is_directed(M.graph)
            v[dst(e)] +=
                log(
                    M.manifold,
                    _read(M, rep_size, x, dst(e)),
                    _read(M, rep_size, x, src(e))
                )
        end
    end
    return v
end
function incident_log!(
    M::GraphManifold{<:AbstractSimpleWeightedGraph,<:Manifold,VertexManifold},
    v,
    x,
)
    rep_size = representation_size(M.manifold)
    for e in edges(M.graph)
        v[src(e)] += (
            get_weight(M.graph, src(e), dst(e)) *
                log(
                    M.manifold,
                    _read(M, rep_size, x, src(e)),
                    _read(M, rep_size, x, dst(e))
                )
        )
        if !is_directed(M.graph)
            v[dst(e)] += (
                get_weight(M.graph, dst(e), src(e)) *
                    log(
                        M.manifold,
                        _read(M, rep_size, x, dst(e)),
                        _read(M, rep_size, x, src(e))
                    )
            )
        end
    end
    return v
end

@doc doc"""
    manifold_dimension(N::GraphManifold{G,M,VertexManifold})

returns the manifold dimension of the [`GraphManifold`](@ref) `N` on the vertices of
a graph $G=(V,E)$, i.e.
````math
d_{\mathcal N} = \lvert V \rVert d_{\mathcal M}.
````
"""
function manifold_dimension(M::VertexGraphManifold)
    return manifold_dimension(M.manifold) * nv(M.graph)
end
@doc doc"""
    manifold_dimension(N::GraphManifold{G,M,EdgeManifold})

returns the manifold dimension of the [`GraphManifold`](@ref) `N` on the edges of
a graph $G=(V,E)$, i.e.
````math
d_{\mathcal N} = \lvert E \rVert d_{\mathcal M}.
````
"""
function manifold_dimension(M::EdgeGraphManifold)
    return manifold_dimension(M.manifold) * ne(M.graph)
end


text:
#=
Visualisation of a Galton Board, also known as a Bean Machine or quincunx (https://en.wikipedia.org/wiki/Bean_machine).

Invented by Francis Galton. the Galton Board is used to demonstrate how the central limit theorem works.
Specifically, it is used to show that for_ a sufficiently large sample size, the binomial distribution
will approximate the normal distribution.
=#

using Javis
using Colors
using Random

"""
    ground(args...)

Define the background and the initial "foreground" color.
"""
function ground(args...)
    background("black")
    sethue("white")
end

"""
    draw_separators(n, xpos, ypos, offset, height; color = "white")

Draws `n` vertical lines of length `height`, each spaced `offset` pixels apart. `xpos` represents the
x-coordinate of the leftmost line, and `ypos` represents the y-coordinate of the upper point used to
construct each line.
"""
function draw_separators(n, xpos, ypos, offset, height; color = "white")
    for i in 1:n
        x = xpos + offset * (i - 1)
        Object(JLine(Point(x, ypos), Point(x, ypos + height); color = color))
    end
end

"""
    draw_pins(n, xpos, ypos, offset; color = "white")

Draws `n` circles of radius 3, each spaced `offset` pixels apart. `xpos` represents the x-coordinate of
the leftmost circle, and `ypos` represents the y-coordinate of the centre of each circle.
"""
function draw_pins(n, xpos, ypos, offset; color = "white")
    for i in 1:n
        Object(
            JCircle(Point(xpos + offset * (i - 1), ypos), 3; color = color, action = :fill),
        )
    end
end

"""
    flip_coin(p)

Simulates flipping a coin, returning 1 with probability `p` and 1 with probability `(1 - p)`. Helper function
used in `move_ball`.
"""
function flip_coin(p = 0.5)
    return rand() < p ? 1 : -1
end

"""
    move_ball(ball, first, last, offset, p, fno)

Moves `ball` from the `first` level of pins to the `last` level of pins (i.e. top to bottom). `offset` should
be the same offset used when creating the bins using `draw_separators`. The value `p` is the probability that
a ball will go right at each pin, and the value `fno` is the frame number from which the ball should start moving.

Returns (fno, pos), where `fno` represents the last frame of the ball when it has reached it's destination,
and pos represents the "sum" of the direction movements. For example, a ball that moved 5 times left, and
7 times right would have a "sum" of -5 + 7 = 2.
"""
function move_ball(ball, first, last, offset, p, fno)
    pos = 0
    act!(ball, Action((fno + 1):(fno + 5), anim_translate(0, offset / 2)))
    fno += 5
    for i in first:last
        direction = flip_coin(p)
        pos += direction
        act!(ball, Action((fno + 1):(fno + 5), anim_translate(direction * offset / 2, 0)))
        fno += 5
        act!(ball, Action((fno + 1):(fno + 5), anim_translate(0, offset)))
        fno += 5
    end
    return (fno, pos)
end

"""
    galton(seed)

Create the Galton Board animation and saves it in the file `galton.gif`. `seed` is used to determine the seed
value for the pseudorandom number generation, and is used for reproducability.
"""
function galton(seed, p)
    nframes = 1200
    video = Video(250, 600)
    Background(1:nframes, ground)

    logocolors = Colors.JULIA_LOGO_COLORS # Use the Julia logo colors
    Random.seed!(seed)

    # Parameters
    n = 15
    xpos = -100
    ypos = -75
    width = 200
    gap = 15
    height = 350
    offset = width / (n - 1)
    radius = offset / 4

    # Draws the line separators or "bins"
    draw_separators(n, xpos, ypos + gap, offset, height)

    # Create "base" of the board
    Object(
        JRect(
            Point(xpos - offset, ypos + gap + height),
            width + 2 * offset,
            7;
            color = logocolors.purple,
            action = :fill,
        ),
    )

    # Draws the pins starting from the bottom most level (closest to base)
    for level in n:-1:1
        draw_pins(level, xpos, ypos, offset; color = logocolors.blue)
        xpos += offset / 2
        ypos -= offset
    end

    nballs = 200
    colors = range(logocolors.green, logocolors.red, length = nballs)
    # Frequency table (i.e. bins), used to determine how much to translate each ball.
    # Total of n + 1 bins.
    slots = zeros(Int, n + 1)
    for i in 1:nballs
        # Set start frame of each ball, interval of 5 frames between each ball
        fno = 5 * i
        # Start each ball at the topmost level
        ball = Object(JCircle(Point(0, ypos), radius; color = colors[i], action = :fill))
        # Move the ball to the top of a bin
        fno_end, pos = move_ball(ball, 1, n, offset, p, fno)

        # Map from "sum" of directions to bin number
        index = (pos + n) ÷ 2 + 1
        # Move each ball into it's "bin"
        act!(
            ball,
            Action(
                (fno_end + 1):(fno_end + 10),
                anim_translate(0, height - 2 * radius * slots[index]),
            ),
        )
        # Update frequency table
        slots[index] += 1
    end

    render(video; pathname = "galton.gif", framerate = 30)
    return slots
end

# Generate the "galton.gif" file.
galton(94, 0.5)


text:
include("utils.jl");

# Source to Source Reverse Mode
# =============================
#
# [Forward mode](./forward.ipynb) works well because all of the symbolic
# operations happen at Julia's compile time; Julia can then optimise the
# resulting program (say, by applying SIMD instructions) and we get very fast
# derivative code. Although we can differentiate Julia code by [compiling it to
# a Wengert list](./tracing.ipynb), we'd be much better off if we could handle
# Julia code directly; then reverse mode can benefit from these optimisations
# too.
#
# However, Julia code is much more complex than a Wengert list, with constructs
# like control flow, data structures and function calls. To do this we'll have
# to handle each of these things in turn.
#
# The first thing to realise is that Julia code is much closer to a Wengert list
# than it looks. Despite its rich syntax, the compiler works with a Wengert-like
# format. The analyses and optimisations that compilers already carry out also
# benefit from this easily-work-with structure.

f(x) = x / (1 + x^2)

@code_typed f(1.0)

# Code with control flow is pnly a little different. We add `goto` statements
# and a construct called the "phi function"; the result is called [SSA
# form](https://en.wikipedia.org/wiki/Static_single_assignment_form).

function pow(x, n)
  r = 1
  while n > 0
    n -= 1
    r *= x
  end
  return r
end

pow(2, 3)

@code_typed pow(2, 3)

# The details of this format are not too important. SSA form is powerful but
# somewhat fiddly to work with in practice, so the aim of this notebook is
# to give a broad intuition for how we handle this.

max_stars_repo_name: backwardn/diff-zoo

max_issues_repo_name: backwardn/diff-zoo

max_forks_repo_name: backwardn/diff-zoo


text:
using TenPuzzle
using Test

@testset "Basic fetures" begin
  @testset "Typical problems" begin
      @test maketen(10) == "10"
      @test maketen(5, 5) == "5 + 5"

      @test maketen(1, 1, 9, 9) in (
          "(1 + 1 / 9) * 9", "(1 / 9 + 1) * 9",
          "9 * (1 + 1 / 9)", "9 * (1 / 9 + 1)"
      )
  end

end


text:
export compute_planestressplasticity_stress_and_stiffness_matrix
@doc raw"""
    planestressplasticity_get_stress_and_sensitivity(
        ε::Array{Float64}, ε0::Array{Float64}, σ0::Array{Float64},
        α::Array{Float64}, K::Float64, σY::Float64, H::Array{Float64}
    )

Returns
- The stress $\sigma$
- The sensitivity $\frac{d\sigma}{d\varepsilon}$
- New internal variable $\alpha_1$

Yield function is defined as
$$y = f(\sigma) - (\sigma_Y+K\alpha)$$
where $\alpha$ is the internal variable per Gauss points. `f` is the J2 function.

$\varepsilon$ is the proposed strain. $\varepsilon_0$ and $\sigma_0$ are strain and stress from the last step.
"""
function compute_planestressplasticity_stress_and_stiffness_matrix(
    u::Array{Float64}, ε0::Array{Float64}, σ0::Array{Float64},
    α::Array{Float64}, K::Float64, σY::Float64, H::Array{Float64}, m::Int64, n::Int64, h::Float64
)
    I = Int64[]; J = Int64[]; V = Float64[]
    function add(ii, jj, kk)
        for  i = 1:length(ii)
            for j = 1:length(jj)
                push!(I,ii[i])
                push!(J,jj[j])
                push!(V,kk[i,j])
            end
        end
    end
    ε = eval_strain_on_gauss_pts(u, m, n, h)
    out_rhs = zeros(2(m+1)*(n+1))
    out_α = zeros(4*m*n)
    out_σ = zeros(4*m*n, 3)

    B = zeros(4, 3, 8)
    for i = 1:2
        for j = 1:2
            ξ = pts[i]; η = pts[j]
            B[(i-1)*2+j,:,:] = [
                -1/h*(1-η) 1/h*(1-η) -1/h*η 1/h*η 0.0 0.0 0.0 0.0
                0.0 0.0 0.0 0.0 -1/h*(1-ξ) -1/h*ξ 1/h*(1-ξ) 1/h*ξ
                -1/h*(1-ξ) -1/h*ξ 1/h*(1-ξ) 1/h*ξ -1/h*(1-η) 1/h*(1-η) -1/h*η 1/h*η
            ]
        end
    end

    for i = 1:m
        for j = 1:n
            idx = [(j-1)*(m+1)+i;(j-1)*(m+1)+i+1;j*(m+1)+i;j*(m+1)+i+1]
            idx = [idx; idx .+ (m+1)*(n+1)]
            for p = 1:2
                for q = 1:2
                    k = (p-1)*2 + q
                    gauss_k = 4*((j-1)*m+i-1)+k
                    Bk = B[k, :, :]
                    σ, dΔσΔε, out_α[gauss_k] = _planestressplasticity_get_stress_and_sensitivity(
                        ε[gauss_k,:], ε0[gauss_k,:], σ0[gauss_k,:], α[gauss_k], K, σY, H
                    )
                    out_rhs[idx] += (σ'*Bk)[:]*0.25*h^2 # 1 x 8
                    out_σ[gauss_k,:] = σ
                    add(idx, idx, Bk'*dΔσΔε*Bk*0.25*h^2)
                end
            end
        end
    end
    return out_rhs, sparse(I, J, V, (m+1)*(n+1)*2, (m+1)*(n+1)*2), out_α, out_σ
end

function _f(σ, α, σY, K)
    return sqrt(σ[1]^2-σ[1]*σ[2]+σ[2]^2+3*σ[3]^2)-σY-K*α
end

function _fσ(σ)
    σ1, σ2, σ3 = σ[1], σ[2], σ[3]
    J2 = sqrt(σ1^2-σ1*σ2+σ2^2+3*σ3^2)
    z = [(σ1 - σ2/2)/J2;
        (-σ1/2 + σ2)/J2;
        3*σ3/J2]
end

function _fσσ(σ)
    σ1, σ2, σ3 = σ[1], σ[2], σ[3]
    J2 = sqrt(σ1^2-σ1*σ2+σ2^2+3*σ3^2)
    [     (-σ1 + σ2/2)*(σ1 - σ2/2)/J2^3 + 1/J2 (σ1/2 - σ2)*(σ1 - σ2/2)/J2^3 - 1/(2*J2)                                   -3*σ3*(σ1 - σ2/2)/J2^3;
    (-σ1 + σ2/2)*(-σ1/2 + σ2)/J2^3 - 1/(2*J2)    (-σ1/2 + σ2)*(σ1/2 - σ2)/J2^3 + 1/J2                                  -3*σ3*(-σ1/2 + σ2)/J2^3;
    3*σ3*(-σ1 + σ2/2)/J2^3                                                        3*σ3*(σ1/2 - σ2)/J2^3 -9*σ3^2/J2^3 + 3/J2]
end

function _planestressplasticity_get_stress_and_sensitivity(
    ε::Array{Float64}, ε0::Array{Float64}, σ0::Array{Float64},
    α0::Float64, K::Float64, σY::Float64, H::Array{Float64}
)
    local dΔσdΔε
    σ = σ0 + H*(ε-ε0)
    Δγ = 0.0
    r2 = _f(σ, α0, σY, K)
    if r2<0
        # @info "elasticity"
        dΔσdΔε = H
    else
        # @info "plasticity"
        function compute(σ, Δγ)
            α = α0 + Δγ
            r1 = σ + Δγ * H* _fσ(σ) - σ0 - H*ε + H*ε0
            r2 = _f(σ, α, σY, K)
            J = [UniformScaling(1.0)+Δγ*H*_fσσ(σ) H*_fσ(σ)
                reshape(_fσ(σ),1,3) -K]
            return [r1;r2], J
        end

        function compute_sensitivity(σ, Δγ)
            α = α0 + Δγ
            J = [UniformScaling(1.0)+Δγ*H*_fσσ(σ) H*_fσ(σ)
                reshape(_fσ(σ),1,3) -K]
            δ = J\[H;zeros(1,3)]
            return δ[1:3,:]
        end
        res0, _ = compute(σ, Δγ)
        for i = 1:100
            res, J = compute(σ, Δγ)
            δ = -J\res
            σ += δ[1:3]; Δγ += δ[4]
            if norm(res)/norm(res0) < 1e-6 || norm(res) < 1e-6
                break
            end
            if i==100
                @warn "Element not converging..."
            end
        end

        dΔσdΔε = compute_sensitivity(σ, Δγ)
    end
    α = α0+Δγ
    return σ[:], dΔσdΔε, α
end