text:
export pick_ins2Opt, pick_ins, pick_ins2, pick_ins2N, pick_ins1

# Tricky implementation to keep the predecessor iterate for inclusion in the secant formula
#
#

let (tpred, ϕpred, dϕpred) = (Inf, Inf, Inf)
    γ = 0.7
    γ2 = 0.999
    global function resetSec!(;t::Float64 = Inf, vϕ::Float64 = Inf, dϕ::Float64 = Inf)
        (tpred, ϕpred, dϕpred) = (t, vϕ, dϕ)
    end

    # test version that filters first the first secant step and then tries the second.
    # Should give the same iterates, slightly faster.

#    """
#        pick_ins2Opt is the scheme which picks the secant point whenever possible and improves it
#        using a pseudo further secant iteration using an Hermite interpolant of ϕ using ϕand dϕ at
#        the points defining the secant.
#    """
    global function pick_ins2Opt(a, b, ϕ, ϕa, ϕb, dϕa, dϕb, forward)
        L = b - a
        if forward
            s = tpred - a
            y = dϕpred - dϕa
            dS = - (dϕa * s / y)
            t2 = a + dS
            # first check secant improvement
            if (dS > 0) & (dS < γ*L)
                # prepare to improve the secant step
                α = coefHIN3([a;tpred],[ϕa dϕa; ϕpred dϕpred])
                p = Polynomial(α)
                dp = Polynomials.derivative(p)
                dϕ2 = dp(dS)
                y2 = dϕ2 - dϕa
                d2 = - (dϕ2 * dS / y2)
                dS2 = dS + d2
                if (dS2 > 0) & (dS2 < γ*L)
                    t2 = a + dS2
                else
                end
            elseif b<Inf
                t2 = (a + b)/2.0
            else
                t2 = pick_inInf(a)
            end
            tpred, ϕpred, dϕpred = a, ϕa, dϕa
            return t2
        else
            s = tpred - b
            y = dϕpred - dϕb
            dS =  - (dϕb * s / y)
            t2 = b + dS
            # first check secant improvement
            if (dS < 0) & (dS > -γ*L)
                # prepare to improve the secant step
                α = coefHIN3([b;tpred],[ϕb dϕb; ϕpred dϕpred])
                p = Polynomial(α)
                dp = Polynomials.derivative(p)
                dϕ2 = dp(dS)
                y2 = dϕ2 - dϕb
                d2 = - (dϕ2 * dS / y2)
                dS2 = dS + d2
                if (dS2 < 0) & (dS2 > -γ*L)
                    t2 = b + dS2
                else
                    t2 = b + dS
                end
            else
                t2 = (a + b)/2.0
            end
            tpred, ϕpred, dϕpred = b, ϕb, dϕb
            return t2

        end
    end

    """
    pick\\_ins is the picking scheme which picks the secant point whenever possible.
    """
    global function pick_ins(a, b, ϕ, ϕa, ϕb, dϕa, dϕb, forward)
        if b < Inf
            L = b - a
            if forward
                s = tpred - a
                y = dϕpred - dϕa
                dS = - (dϕa * s / y)
                tpred, dϕpred = a, dϕa
                anext = a + dS
                if (dS > 0) & (dS < γ*L)
                    return anext
                else
                    return (a + b)/2.0
                end
            else
                s = tpred - b
                y = dϕpred - dϕb
                dS =  - (dϕb * s / y)
                tpred, dϕpred = b, dϕb
                bnext = b + dS
                if (dS < 0) & (dS > -γ*L)
                    return bnext
                else
                    return (a + b)/2.0
                end
            end
        else
            return pick_inInf(a)
        end
    end

#    """
#            pick\\_ins2 is the picking scheme which picks the secant point whenever possible and improve#s it
#        using a pseudo further secant iteration using an Hermite interpolant of ϕ using ϕ and dϕ at
#        the points defining the secant.
#    """
    global function pick_ins2(a, b, ϕ, ϕa, ϕb, dϕa, dϕb, forward)
        L = b - a
        if forward
            s = tpred - a
            y = dϕpred - dϕa
            dS = - (dϕa * s / y)
            t2 = a + dS
            α = coefHIN3([a;tpred],[ϕa dϕa; ϕpred dϕpred])
            p = Polynomial(α)
            dp = Polynomials.derivative(p)
            dϕ2 = dp(dS)
            y2 = dϕ2 - dϕa
            d2 = - (dϕ2 * dS / y2)
            dS2 = dS + d2

            anext = a + dS2

            tpred, ϕpred, dϕpred = a, ϕa, dϕa
            if (dS2 > (1-γ2)*L) & (dS2 < γ*L)
                return anext
            elseif (dS > (1-γ2)*L) & (dS < γ*L)
                return a + dS
            elseif b<Inf
                return (a + b)/2.0
            else
                return pick_inInf(a)
            end
        else
            s = tpred - b
            y = dϕpred - dϕb
            dS =  - (dϕb * s / y)
            t2 = b + dS

            α = coefHIN3([b;tpred],[ϕb dϕb; ϕpred dϕpred])
            p = Polynomial(α)
            dp = Polynomials.derivative(p)
            dϕ2 = dp(dS)

            y2 = dϕ2 - dϕb
            d2 = - (dϕ2 * dS / y2)
            dS2 = dS + d2

            bnext = b + dS2

            tpred, ϕpred, dϕpred = b, ϕb, dϕb
            if (dS2 < (1-γ2)*L) & (dS2 > -γ*L)
                return bnext
            elseif (dS < (1-γ2)*L) & (dS > -γ*L)
                return b + dS
            else
                return (a + b)/2.0
            end
        end
    end

#    """
#        pick\\_ins2N is the picking scheme picking the secant point whenever possible and improves it
#        using a pseudo further Newton iteration using an Hermite interpolant of ϕ using ϕ and dϕ at
#        the points defining the secant.
#    """
    global function pick_ins2N(a, b, ϕ, ϕa, ϕb, dϕa, dϕb, forward)
        L = b - a
        if forward
            s = tpred - a
            y = dϕpred - dϕa
            dS = - (dϕa * s / y)
            t2 = a + dS
            α = coefHIN3([a;tpred],[ϕa dϕa; ϕpred dϕpred])
            p = Polynomial(α)
            dp = Polynomials.derivative(p)
            d2p = Polynomials.derivative(dp)
            dϕ2 = dp(dS)
            d2ϕ2 = d2p(dS)
            d2 = - (dϕ2 / d2ϕ2)
            dS2 = dS + d2

            anext = a + dS2

            tpred, ϕpred, dϕpred = a, ϕa, dϕa
            if (dS2 > 0) & (dS2 < γ*L)
                return anext
            elseif (dS > 0) & (dS < γ*L)
                return a + dS
            elseif b<Inf
                return (a + b)/2.0
            else
                return pick_inInf(a)
            end
        else
            s = tpred - b
            y = dϕpred - dϕb
            dS =  - (dϕb * s / y)
            t2 = b + dS

            α = coefHIN3([b;tpred],[ϕb dϕb; ϕpred dϕpred])
            p = Polynomial(α)
            dp = Polynomials.derivative(p)
            d2p = Polynomials.derivative(dp)
            dϕ2 = dp(dS)

            d2ϕ2 = d2p(dS)
            d2 = - (dϕ2 / d2ϕ2)
            dS2 = dS + d2

            bnext = b + dS2

            tpred, ϕpred, dϕpred = b, ϕb, dϕb
            if (dS2 < 0) & (dS2 > -γ*L)
                return bnext
            elseif (dS < 0) & (dS > -γ*L)
                return b + dS
            else
                return (a + b)/2.0
            end
        end
    end

#    """
#        pick\\_ins2 is the picking scheme which picks the secant point whenever possible and improves it
#        using a pseudo further secant iteration using an Hermite interpolant of ϕ using ϕ and dϕ at
#        the points defining the secant.
#    """
    global function pick_ins1(a, b, ϕ, ϕa, ϕb, dϕa, dϕb, forward)
        L = b - a
        if forward
            s = tpred - a
            y = dϕpred - dϕa
            dS = - (dϕa * s / y)
            t = a + dS

            tpred, ϕpred, dϕpred = a, ϕa, dϕa
            if (dS > (1-γ2)*L) & (dS < γ*L)
                return t
            elseif b<Inf
                return (a + b)/2.0
            else
                return pick_inInf(a)
            end
        else
            s = tpred - b
            y = dϕpred - dϕb
            dS =  - (dϕb * s / y)
            t = b + dS
            tpred, ϕpred, dϕpred = b, ϕb, dϕb
            if (dS < (1-γ2)*L) & (dS > -γ*L)
                return t
            else
                return (a + b)/2.0
            end
        end
    end

end  #  let

# Repeat the doc, the actual behavior of docstrings is not as documented:(

#    """
#    pick\\_ins is the picking scheme which picks the secant point whenever possible.
#    """
#pick_ins

    """
        pick\\_ins2N is the picking scheme picking the secant point whenever possible and improves it
        using a pseudo further Newton iteration using an Hermite interpolant of ϕ using ϕ and dϕ at
        the points defining the secant.
    """
pick_ins2N

    """
        pick_ins2Opt is the scheme which picks the secant point whenever possible and improves it
        using a pseudo further secant iteration using an Hermite interpolant of ϕ using ϕand dϕ at
        the points defining the secant.
    """
pick_ins2Opt

    """
        pick\\_ins1 is the picking scheme which picks the secant point whenever possible and improves it
        using a pseudo further secant iteration using an Hermite interpolant of ϕ using ϕ and dϕ at
        the points defining the secant.
    """
pick_ins1

"""
pick\\_ins2 is the picking scheme which picks the secant point whenever possible and improves it
        using a pseudo further secant iteration using an Hermite interpolant of ϕ using ϕ and dϕ at
        the points defining the secant.
"""
pick_ins2


text:
"""
Example of using HttpServer with HTTPS protocol:

1) Create a self-signed SSL Certificate
    - Generate a private key
        openssl genrsa -des3 -out server.key 1024

    - Generate a CSR (Certificate Signing Request)
        openssl req -new -key server.key -out server.csr -subj "/L=github.com/O=JuliaWeb/OU=HttpServer/CN=localhost"

    - Remove Passphrase from Key
        cp server.key server.key.org
        openssl rsa -in server.key.org -out server.key

    - Generating a Self-Signed Certificate
        openssl x509 -req -days 365 -in server.csr -signkey server.key -out server.crt

2) Add certificate/private key pair to a certificate store

3) Run server with the store which contains self-signed SSL certificate
"""

# Adding certificate/private key pair to certificate store
using GnuTLS
cert_store = GnuTLS.CertificateStore()
GnuTLS.load_certificate(cert_store, "server.crt", "server.key", true)

# Setting up simple server
using HttpServer
http = HttpHandler() do req::Request, res::Response
    Response("Hello Secure World!")
end
http.events["listen"] = (saddr) -> println("Running on https://$saddr (Press CTRL+C to quit)")

# Running server in HTTPS mode
server = Server(http)
run(server, port=8000, ssl=cert_store)

ext: jl


text:
"""
    PowerMethod way of finding the leading boundary mps
"""
@with_kw struct PowerMethod <: Algorithm
    tol_galerkin::Float64 = Defaults.tol
    tol_gauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    orthmaxiter::Int = Defaults.maxiter
    finalize::Function = (iter,state,ham,pars) -> (state,pars);
    verbose::Bool = Defaults.verbose
end

function leading_boundary(state::MPSMultiline, H,alg::PowerMethod,pars=params(state,H))
    galerkin  = 1+alg.tol_galerkin
    iter       = 1

    newAs = similar(state.AL)

    while true

        @threads for col in 1:size(state,2)

            vac = let state=state,pars=pars
                circshift([ac_prime(ac,row,col,state,pars) for (row,ac) in enumerate(state.AC[:,col])],1)
            end
            vc  = let state=state,pars=pars
                circshift([c_prime(c,row,col,state,pars) for (row,c) in enumerate(state.CR[:,col])],1)
            end

            for row in 1:size(state,1)
                QAc,_ = leftorth!(vac[row], alg=TensorKit.Polar())
                Qc,_  = leftorth!(vc[row], alg=TensorKit.Polar())
                newAs[row,col] = QAc*adjoint(Qc)
            end

        end

        state = MPSMultiline(newAs; leftgauged=true,tol = alg.tol_gauge, maxiter = alg.orthmaxiter)
        galerkin   = calc_galerkin(state, pars)
        alg.verbose && println("powermethod @iteration $(iter) galerkin = $(galerkin)")

        if galerkin <= alg.tol_galerkin || iter>=alg.maxiter
            iter>=alg.maxiter && println("powermethod didn't converge $(galerkin)")
            return state, pars, galerkin
        end

        (state,pars) = alg.finalize(iter,state,H,pars);

        iter += 1
    end
end


text:
"""
    AbstractRanking

Supertype for rankings.
"""
abstract type AbstractRanking end

"""
    Ranking{Ts<:Real} <: AbstractRanking

Ranking.
"""
struct Ranking{Ts<:Real} <: AbstractRanking
    e::Set{Int}
    q::Set{Int}
    S::Matrix{Ts}
    rl::Vector{Int}
    rs::Vector{Ts}

    function Ranking{Ts}(e::Set{Int}, q::Set{Int}, S::Matrix{Ts},
                        rl::Vector{Int}, rs::Vector{Ts}) where {Ts<:Real}
        return new(e, q, S, rl, rs)
    end
end
Ranking(e::Set{Int}, q::Set{Int}, S::Matrix{Ts}, rl::Vector{Int},
        rs::Vector{Ts}) where {Ts} =
    Ranking{Ts}(e, q, S, rl, rs)


text:
module Tools

using Plots
pyplot()

using Statistics
using StatsBase
using Bootstrap

using ..Sessions

export learning_curve, Uniform, curve_conf, curves, average_series

function moving_average(a::AbstractVecOrMat, n::Int)
    @assert 1 <= n <= size(a, 1)
    out = similar(a, size(a, 1) - n + 1, size(a, 2))
    out[1, :] = mean(a[1:n, :], dims=1)
    for i in axes(out, 1)[2:end]
        out[i, :] = (n * out[i - 1, :] .- a[i - 1, :] .+ a[i + n - 1, :]) / n
    end
    return out
end

# Plot a learning curve from the results of a experiment
function learning_curve(experiment_result, ma::Int=1)
    smoothed = moving_average(experiment_result, ma)

    mean_series = mean(smoothed; dims=2)
    sd = std(smoothed; dims=2, mean=mean_series)

    plot(mean_series, ribbon=sd, legend=false)
end

function average_series(returns::AbstractMatrix; bs=1000, ci=0.95, window=5)
    bs_mean = Matrix{Float64}(undef, size(returns, 1), 3)
    for (i, samples) in Iterators.enumerate(eachrow(returns))
        bs_mean[i, :] .= confint(
            bootstrap(
                mean,
                samples,
                BasicSampling(bs)
            ),
            BasicConfInt(ci)
        )[1]
    end

    bs_mean = moving_average(bs_mean, window)

    m = bs_mean[:, 1]
    lo = abs.(m - bs_mean[:, 2])
    hi = abs.(m - bs_mean[:, 3])

    return m, lo, hi
end

average_series(trial::Trial; bs=1000, ci=0.95, window=5) = average_series(trial.cummulative_rewards; bs, ci, window)

# function average_series(trial; bs=1000, ci=0.95, window=5)
#     bs_mean = Matrix{Float64}(undef, size(trial.cummulative_rewards, 1), 3)
#     for (i, samples) in Iterators.enumerate(eachrow(trial.cummulative_rewards))
#         bs_mean[i, :] .= confint(
#             bootstrap(
#                 mean,
#                 samples,
#                 BasicSampling(bs)
#             ),
#             BasicConfInt(ci)
#         )[1]
#     end

#     bs_mean = moving_average(bs_mean, window)

#     m = bs_mean[:, 1]
#     lo = abs.(m - bs_mean[:, 2])
#     hi = abs.(m - bs_mean[:, 2])

#     return m, lo, hi
# end

function curve_conf(trial; bs=1000, ci=0.95, window=5)
    m, lo, hi = average_series(trial; bs, ci, window)

    plot(m, ribbon=(lo, hi), legend=false)
end

function curves(trials...; window=5)
    bs = 1000
    ci = 0.95
    p = plot(ylims=(0, 500))

    for trial in trials
        bs_mean = Matrix{Float64}(undef, size(trial.cummulative_rewards, 1), 3)
        for (i, samples) in Iterators.enumerate(eachrow(trial.cummulative_rewards))
            bs_mean[i, :] .= confint(
                bootstrap(
                    mean,
                    samples,
                    BasicSampling(bs)
                ),
                BasicConfInt(ci)
            )[1]
        end

        bs_mean = moving_average(bs_mean, window)

        m = bs_mean[:, 1]
        lo = abs.(m - bs_mean[:, 2])
        hi = abs.(m - bs_mean[:, 3])

        plot!(m, ribbon=(lo, hi), label=trial.name)
    end

    return p
end

learning_curve(s::Session, ma::Int=1) = learning_curve(s.cummulative_rewards, ma)

learning_curve(t::Trial, ma::Int=1) = learning_curve(t.cummulative_rewards, ma)

end


text:
module SparseADRules

using ChainRulesCore: include
using Base: promote_eltype
using LinearAlgebra
using SparseArrays
using NiLang
using NiLang.AD
using ChainRulesCore

include("compat.jl")
include("linalg.jl")
include("lowranksvd.jl")
include("srcutils.jl")
include("chainrules.jl")
end


text:
module TestMain_VariableTimestep

using Test
using Mimi

import Mimi:
    reset_compdefs, reset_variables, @defmodel,
    variable, variable_names, external_param, build,
    compdef, compdefs, dimensions, dimension, compinstance

reset_compdefs()

@defcomp foo1 begin
    index1 = Index()

    par1 = Parameter()
    par2::Bool = Parameter(index=[time,index1], description="description par 1")
    par3 = Parameter(index=[time])

    var1 = Variable()
    var2 = Variable(index=[time])
    var3 = Variable(index=[time,index1])

    idx3 = Index()
    idx4 = Index()
    var4::Bool = Variable(index=[idx3])
    var5 = Variable(index=[index1, idx4])
end

@defmodel x1 begin
    index[index1] = [:r1, :r2, :r3]
    index[time] = [2010, 2015, 2030]
    index[idx3] = 1:3
    index[idx4] = 1:4
    component(foo1)

    foo1.par1 = 5.0
    # foo1.par2 = []
end

@test length(compdefs()) == 4   # adder, 2 connectors, and foo1

# x1 = foo1(Float64, Dict{Symbol, Int}(:time=>10, :index1=>3))
# x1 = foo1(Float64, Val{1}, Val{1}, Val{10}, Val{1}, Val{1}, Val{1}, Val{1}, Dict{Symbol, Int}(:time=>10, :index1=>3))

@test length(dimension(x1.md, :index1)) == 3

# @test_throws MethodError x1.Parameters.par1 = Array{Float64}(1,2)

par1 = external_param(x1, :par1)
@test par1.value == 5.0

set_param!(x1, :foo1, :par1, 6.0)
par1 = external_param(x1, :par1)
@test par1.value == 6.0

set_param!(x1, :foo1, :par2, [true true false; true false false; true true true])

set_param!(x1, :foo1, :par3, [1.0, 2.0, 3.0])

build(x1)

ci = compinstance(x1, :foo1)
reset_variables(ci)

# Check all variables are defaulted
@test isnan(get_var_value(ci, :var1))

m = Model()
set_dimension!(m, :time, 20)
set_dimension!(m, :index1, 5)
add_comp!(m, foo1)

@test :var1 in variable_names(x1, :foo1)

end #module


text:
# Loss function, gradients and hessian vector products for SVMs with quadratic loss

println("Loading SVM module")

function loss(data::Matrix, theta::Vector, lambda::Float64)
  result = 0.0
  n = size(data, 2)
  for i = 1:n
    x = data[1:end-1,i]
    y = data[end,i]
    result += 0.5 * max(0.0, 1 - y * dot(vec(x), theta))^2/n
  end
  result += lambda/2 * dot(theta, theta)
  return result
end

function grad!(datapoint::Vector, theta::Vector, lambda::Float64, grad::Vector)
  x = datapoint[1:end-1]
  y = datapoint[end]
  alpha = dot(x, theta)
  if y * alpha < 1.0
    BLAS.axpy!(1.0, lambda * theta -  (1 - y*alpha) * y * x, grad)
  else
    BLAS.axpy!(1.0, lambda * theta, grad)
  end
end

function hvp!(datapoint::Vector, theta::Vector, v::Vector, lambda::Float64, result::Vector)
  x = datapoint[1:end-1]
  y = datapoint[end]
  alpha = dot(x, theta)
  if y * alpha < 1.0
    BLAS.axpy!(1.0, x * dot(x, v) + lambda * v, result)
  else
    BLAS.axpy!(1.0, lambda * v, result)
  end
end

require("checks.jl")
ad_hoc_check("svm")

max_stars_repo_name: pcmoritz/slbfgs

max_issues_repo_name: pcmoritz/slbfgs

max_forks_repo_name: pcmoritz/slbfgs


text:
using GridapP4est
using Documenter

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

makedocs(;
    modules=[GridapP4est],
    authors="Alberto F. Martin <alberto.martin@monash.edu>",
    repo="https://github.com/gridap/GridapP4est.jl/blob/{commit}{path}#{line}",
    sitename="GridapP4est.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gridap.github.io/GridapP4est.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/gridap/GridapP4est.jl",
    devbranch="main",
)


text:
## Wrap to standard API spec

# Leaving as comment for future plans, see https://github.com/NavAbility/nva-sdk/issues/21
# TODO What to use:
# - NavAbilityPlatform
# - NavAbilityConnection
# - NavAbilityClientContext
# -
# mutable struct NavAbilityPlatform{T}
#     client
#     context::T
# end

# function NavAbilityPlatform(;
#             apiUrl::AbstractString="https://api.navability.io",
#             UserId::AbstractString="guest@navability.io",
#             RobotId::AbstractString=ENV["USER"],
#             SessionId::AbstractString="Session_" * string(uuid4())[1:8])
#     #
#     client = NavAbilityHttpsClient(apiUrl)
#     context = Client(UserId, RobotId, SessionId)
#     return NavAbilityPlatform(client, context)
# end

"""
  addVariable
Add a variable to the NavAbility Platform service
Example
```julia
addVariable(client, context, "x0", NVA.Pose2)
```
"""
function addVariable(client,
                     context,
                     label::Union{<:AbstractString,Symbol},
                     varType::Union{<:AbstractString,Symbol};
                     tags::Vector{String}=String[],
                     timestamp::String = string(now(Dates.UTC))*"Z")
                    # TODO
                    # solvable::Int=1
                    # nanosecondtime,
                    # smalldata,
    #
    union!(tags, ["VARIABLE"])
    v = Variable(string(label), Symbol(varType), tags, timestamp)
    return addVariable(client, context, v)
end

function assembleFactorName(xisyms::Union{Vector{String},Vector{Symbol}})
    return string(xisyms...,"f_",string(uuid4())[1:4])
end

function getFncTypeName(fnc::InferenceType)
    return split(string(typeof(fnc)),".")[end]
end

#TODO solverParams

function addFactor(client,
                   context,
                   xisyms::Union{Vector{String},Vector{Symbol}},
                   fnc::InferenceType;
                   multihypo::Vector{Float64}=Float64[],
                   nullhypo::Float64=0.0,
                   solvable::Int=1,
                   tags::Vector{String}=String[],
                   # timestamp::Union{DateTime,ZonedDateTime}=now(localzone()), #TODO why timestamp difference from IIF
                   timestamp::String = string(now(Dates.UTC))*"Z",
                   inflation::Real=3.0,
                   namestring::String = assembleFactorName(xisyms),
                   nstime::String="0"
                  )
                   # TODO maybe
                   #  graphinit::Bool=true,
                   #  threadmodel=SingleThreaded,
                   #  suppressChecks::Bool=false,
                   #  _blockRecursion::Bool=!getSolverParams(dfg).attemptGradients
    #
    # create factor data
    factordata = FactorData(;fnc, multihypo, nullhypo, inflation)

    fncType = getFncTypeName(fnc)

    union!(tags, ["FACTOR"])
    # create factor
    factor = Factor(
        namestring,
        nstime,
        fncType,
        string.(xisyms),
        factordata,
        solvable,
        tags,
        timestamp,
        DFG_VERSION
    )
    # add factor
    resultId = addFactor(client, context, factor)

    return resultId
end


text:
@testset "Initialize from rest" begin

    # original speedy T30 with 96x48 grid
    P = Parameters(NF=Float64)
    G = GeoSpectral(P)
    B = Boundaries(P,G)

    Prog = SpeedyWeather.initialize_from_rest(P,B,G)

    @test all(Prog.vor .== 0)
    @test all(Prog.div .== 0)

    # test surface layer only at the moment
    temp_grid = gridded(Prog.temp[:,:,end],G)
    pres_surf_grid = gridded(Prog.pres_surf,G)
    humid_grid = gridded(Prog.humid[:,:,end],G)

    # temperature between 200K and 350K everywhere
    @test all(temp_grid .> 200)
    @test all(temp_grid .< 350)

    # surface pressure between log(300hPa) and log(1200hPa) everywhere
    @test all(pres_surf_grid .> log(300))
    @test all(pres_surf_grid .< log(1200))

    # humidity non-negative everywhere
    # humidity has currently values of O(1e9)...
    @test_skip all(humid_grid .>= 0)
end


text:
Pkg.add("FastGaussQuadrature")
Pkg.add("DifferentialEquations")
Pkg.add("Parameters")


text:
#fltk.h
@getf libfltk mgl_create_graph_fltk(draw::DrawCallback, title::Ptr{Cchar}, par::Ptr{Void}, load::LoadCallback)::HMGL
@getf libfltk mgl_fltk_run()::Cint
@getf libfltk mgl_fltk_thr()::Cint


text:
"
  Name     : c4_19_scatter_plot_PyPlot.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 1/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
"
using PyPlot
n=50
x = 100*rand(n)
y = 100*rand(n)
areas = 800*rand(n)
fig = figure("pyplot_scatterplot",figsize=(10,10))
ax = axes()
scatter(x,y,s=areas,alpha=0.5)
title("Scatter Plot")
xlabel("X")
ylabel("Y")
grid("on")


text:
#=
day8:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-08
=#

using AoC, Test
using StructEquality # https://docs.juliahub.com/StructEquality/TwsrV/1.0.0/
using Base.Iterators: flatten
using Memoize

DIGITS = Dict('0'=>"abcefg",'1'=>"cf",'2'=>"acdeg",'3'=>"acdfg",'4'=>"bcdf",'5'=>"abdfg",'6'=>"abdefg",'7'=>"acf",'8'=>"abcdefg",'9'=>"abcdfg")
DIGIT_DECODER = Dict((v=>k) for (k,v) in DIGITS)
LENGTHS = Dict((k=>length(v)) for (k,v) in DIGITS)
UNIQUE_LENGTHS = Dict((v=>k) for (k,v) in filter((pair)->count(l->l==pair[2],values(LENGTHS))==1, LENGTHS))
SEGMENT_NAMES = ['a','b','c','d','e','f','g']

@memoize Dict allpermutations(n::Number) = allpermutations([1:n...])
@memoize Dict allpermutations(n::Vector) = allpermutationssorted(sort(n))
@memoize Dict function allpermutationssorted(n::Vector)
    if length(n) == 1
        return [n]
    end
    result = zeros(Int, length(n),0)
    for i in n
        remainder = setdiff(n, [i])

        for p in eachcol(allpermutations(remainder))
            p = [i;p...]
            result = hcat(result, p)
        end
    end
    return result
end

@test allpermutations(2) == [[1,2] [2,1]]
@test allpermutations(3) == [[1,2,3] [1,3,2] [2,1,3] [2,3,1] [3,1,2] [3,2,1]]
@test size(allpermutations(3),2) == 3*2*1
@test size(allpermutations(4),2) == 4*3*2*1
@test size(allpermutations(5),2) == 5*4*3*2*1

struct Display
    combos :: Vector{String}
    display :: Vector{String}
end
@def_structequal Display

function parseDisplay(line)
    (left,right) = split(line, " | ")
    return Display(split(left," "), split(right, " "))
end

@test parseDisplay("be cfbegad cbdgef fgaecd cgeb fdcge agebfd fecdb fabcd edb | fdgacbe cefdb cefbgd gcbe") ==
    Display(["be","cfbegad","cbdgef","fgaecd","cgeb","fdcge","agebfd","fecdb","fabcd","edb"],["fdgacbe","cefdb","cefbgd","gcbe"])

function evaluateCipher(combos)
    # keys are the displayed segments, values are the actual values
    segmentMappings = Dict(c => Set(SEGMENT_NAMES) for c in SEGMENT_NAMES)
    return Dict(k=>[v...][1] for (k,v) in evaluateCipher(combos, segmentMappings))
end

function evaluateCipherCandidate(displayedCombos, displayedCombo, segmentMappings, expectedCombos)
    if isempty(displayedCombos)
        return segmentMappings
    end

    actualDisplay = [c for c in displayedCombo]
    for expectedDisplay in filter(d -> length(d) == length(displayedCombo), expectedCombos)
        candidateMappingTo = [c for c in expectedDisplay]

        for candidateMappingFrom in eachcol(allpermutations(actualDisplay))
            candidateMappings = copy(segmentMappings)
            okMapping = true
            for c in eachindex(candidateMappingFrom)
                candidateMappings[candidateMappingFrom[c]] = intersect(candidateMappings[candidateMappingFrom[c]], candidateMappingTo[c])
                if isempty(candidateMappings[candidateMappingFrom[c]])
                    okMapping = false
                    break
                end
            end
            if okMapping
                result = evaluateCipher(setdiff(displayedCombos,[displayedCombo]), candidateMappings, setdiff(expectedCombos, [expectedDisplay]))
                if !isnothing(result)
                    return result
                end
            end
        end
    end
    return nothing
end

function evaluateCipher(displayedCombos, segmentMappings, expectedCombos=[values(DIGITS)...])

    if isempty(displayedCombos)
        return segmentMappings
    end

    fixedLengthCombos = filter(c -> length(c) ∈ keys(UNIQUE_LENGTHS), displayedCombos); lt=(x,y)->length(x)<length(y)

    if !isempty(fixedLengthCombos)
        displayedCombo = argmin(c->length(c), fixedLengthCombos)
        return evaluateCipherCandidate(displayedCombos, displayedCombo, segmentMappings, expectedCombos)
    end

    for displayedCombo in displayedCombos
        result = evaluateCipherCandidate(displayedCombos, displayedCombo, segmentMappings, expectedCombos)
        if !isnothing(result)
            return result
        end
    end
    return nothing
end
@test evaluateCipher(["acedgfb","cdfbe","gcdfa","fbcad","dab","cefabd","cdfgeb","eafb","cagedb","ab"]) == Dict('d'=>'a','e'=>'b','a'=>'c','f'=>'d','g'=>'e','b'=>'f','c'=>'g')

decipher(displayedCombos,cipher) = [decipher(d,cipher) for d in displayedCombos]
function decipher(displayedCombo::AbstractString, cipher)
    return String(sort([cipher[c] for c in displayedCombo]))
end
@test decipher("acedgfb", Dict('d'=>'a','e'=>'b','a'=>'c','f'=>'d','g'=>'e','b'=>'f','c'=>'g') ) == "abcdefg"

decode(display) = String([DIGIT_DECODER[d] for d in display])
@test decode(["abdfg"]) == "5"
@test decode(["abdfg","acdfg"]) == "53"

function evaluateDisplay(display::Display)
    cipher = evaluateCipher(display.combos)
    correctedDisplay = decipher(display.display, cipher)
    decode(correctedDisplay)
end

@test evaluateDisplay(Display(["acedgfb","cdfbe","gcdfa","fbcad","dab","cefabd","cdfgeb","eafb","cagedb","ab"],["cdfeb","fcadb","cdfeb","cdbaf"])) == "5353"

part1(lines) = part1([parseDisplay(line) for line in lines])
function part1(displays :: Array{Display})
    allDigits = [flatten(displays .|> d->d.display)...]
    uniqueLengthDigits = filter(allDigits) do digit
        length(digit) ∈ keys(UNIQUE_LENGTHS)
    end
    return length(uniqueLengthDigits)
end

@test part1(exampleLines(8,2)) == 26

part2(lines) = part2([parseDisplay(line) for line in lines])
function part2(displays :: Array{Display})
    displays = [evaluateDisplay(d) for d in displays]
    return sum(parse.(Int, displays))
end

@test part2(exampleLines(8,2)) == 61229

lines(8) |> ll -> @time part1(ll) |> show
exampleLines(8,2) |> ll -> @time part2(ll) |> show
lines(8) |> ll -> @time part2(ll) |> show


text:
using GeoStatsBase
using InverseDistanceWeighting
using Plots, VisualRegressionTests
using Distances
using Test, Pkg

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

@testset "Basic problem" begin
  geodata = georef((variable=[1.,0.,1.],), [25. 50. 75.;  25. 75. 50.])
  domain  = RegularGrid(100,100)
  problem = EstimationProblem(geodata, domain, :variable)

  solver = InvDistWeight(:variable => (neighbors=3,))

  solution = solve(problem, solver)

  if visualtests
    @plottest contourf(solution,size=(800,400)) joinpath(datadir,"solution.png") !istravis
  end
end

@testset "Haversine" begin
  geodata = georef((variable=[4.0,-1.0,3.0],), [50.0 100.0 200.0; -30.0 30.0 10.0])
  domain  = RegularGrid((1.0, -89.0), (359.0, 89.0), dims=(200, 100))
  problem = EstimationProblem(geodata, domain, :variable)

  solver = InvDistWeight(:variable => (distance=Haversine(1.0),))

  solution = solve(problem, solver)

  if visualtests
    @plottest contourf(solution,size=(800,200)) joinpath(datadir,"haversine.png") !istravis
  end
end


text:
using DelimitedFiles
using Plots

tensor = readdlm("C:path")
tensor = reshape(tensor, (420,420,31))
tensor = permutedims(tensor, [2,1,3])
tensor = tensor
heatmap(tensor[:,:,30])
contour(tensor[:,:,30], levels = 0.1:0.5)


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 19cc5ec6-f7cb-4736-98c3-498d91eb2414
using BarycentricInterpolation

# ╔═╡ 53285b2c-3704-4adc-9bdd-9dea619b5140
md" ## BarycentricInterpolation.jl"

# ╔═╡ 7ba7dd56-4ee9-4c5c-97a8-5f99a9ba507e
p = Chebyshev2{20}()           # create a Chebyshev type 2 polynomial of order 20

# ╔═╡ ddacc062-c4cb-40e1-8645-a1f8619eabea
x = nodes(p)                   # get the nodes

# ╔═╡ 7fdcd309-5e3e-4c65-9d8c-96d19fdc923c
y = sinpi.(x)                  # generate y values at the nodes

# ╔═╡ 9c9fdfad-44eb-4f5b-a5e0-82eaade02ed1
x_new = rand()*2 -1            # a random number in [-1, +1]

# ╔═╡ f3f7c1ea-f4aa-4dc2-a630-0fb37e4bdadb
println(interpolate(p, y, x_new) ≈ sinpi(x_new))       # hopefully true!

# ╔═╡ d0ca3867-a266-4299-ac49-9745154b3ad4
D = differentiation_matrix(p)  # get the differentiation matrix

# ╔═╡ 880d1ae8-5a5a-463a-8e35-a73bdc40208a
println(interpolate(p, D*y, x_new) ≈ pi*cospi(x_new))  # hopefully true!

# ╔═╡ c2fcf437-bd86-4b24-8c77-29702cf5373d
md" ## End of BarycentricInterpolation.jl"

# ╔═╡ Cell order:
# ╟─53285b2c-3704-4adc-9bdd-9dea619b5140
# ╠═19cc5ec6-f7cb-4736-98c3-498d91eb2414
# ╠═7ba7dd56-4ee9-4c5c-97a8-5f99a9ba507e
# ╠═ddacc062-c4cb-40e1-8645-a1f8619eabea
# ╠═7fdcd309-5e3e-4c65-9d8c-96d19fdc923c
# ╠═9c9fdfad-44eb-4f5b-a5e0-82eaade02ed1
# ╠═f3f7c1ea-f4aa-4dc2-a630-0fb37e4bdadb
# ╠═d0ca3867-a266-4299-ac49-9745154b3ad4
# ╠═880d1ae8-5a5a-463a-8e35-a73bdc40208a
# ╟─c2fcf437-bd86-4b24-8c77-29702cf5373d


text:
points(S::Jacobi,n)=fromcanonical(S,gaussjacobi(n,S.a,S.b)[1])

immutable JacobiTransformPlan{DD,T,TT}
    space::Jacobi{TT,DD}
    points::Vector{T}
    weights::Vector{T}
end

plan_transform(S::Jacobi,v::Vector) = JacobiTransformPlan(S,gaussjacobi(length(v),S.a,S.b)...)
function transform(S::Jacobi,vals,plan::JacobiTransformPlan)
#    @assert S==plan.space
    x,w = plan.points, plan.weights
    V=jacobip(0:length(vals)-1,S.a,S.b,x)'
    w2=w.*(1-x).^(S.a-plan.space.a).*(1+x).^(S.b-plan.space.b)   # need to weight if plan is different
    nrm=(V.^2)*w2

    V*(w2.*vals)./nrm
end

immutable JacobiITransformPlan{DD,T}
    space::Jacobi{Float64,DD}
    points::Vector{T}
end

plan_itransform(S::Jacobi,cfs::Vector) = JacobiITransformPlan(S,points(S,length(cfs)))
function itransform(S::Jacobi,cfs,plan::JacobiITransformPlan)
#    @assert S==plan.space
    jacobip(0:length(cfs)-1,S.a,S.b,tocanonical(S,plan.points))*cfs
end

function coefficients(f::Vector,a::Jacobi,b::Chebyshev)
    if domain(a) == domain(b) && (!isapproxinteger(a.a-0.5) || !isapproxinteger(a.b-0.5))
        cjt(f,a.a,a.b)
    else
        defaultcoefficients(f,a,b)
    end
end
function coefficients(f::Vector,a::Chebyshev,b::Jacobi)
    if domain(a) == domain(b) && (!isapproxinteger(b.a-0.5) || !isapproxinteger(b.b-0.5))
        icjt(f,b.a,b.b)
    else
        defaultcoefficients(f,a,b)
    end
end

function coefficients(f::Vector,a::Jacobi,b::Jacobi)
    if domain(a) == domain(b) && (!isapproxinteger(a.a-b.a) || !isapproxinteger(a.b-b.b))
        jjt(f,a.a,a.b,b.a,b.b)
    else
        defaultcoefficients(f,a,b)
    end
end


text:
# Test calculations of properties of tensors

using Test, CIJ
using LinearAlgebra: Symmetric

@testset "Properties" begin
    @testset "Is iso" begin
        let ciso = CIJ.iso(vp=8000, vs=4500), caniso = CIJ.thom(8000, 4500, 0.1, 0.1, 0.1)
            @test is_iso(ciso)
            ciso[1,1] = nextfloat(ciso[1,1])
            @test !isapprox(ciso[1,1], ciso[2,2], atol=eps(eltype(ciso)))
            @test !is_iso(ciso)
            @test is_iso(ciso, atol=√eps(eltype(ciso)))
            @test !is_iso(caniso)
        end
    end

    @testset "Is stable" begin
        let c = CIJ.iso(vp=8000, vs=4500), c′ = CIJ.iso(vp=8000, vs=0)
            @test is_stable(c)
            @test is_stable(c.data)
            @test is_stable(Array(c))
            @test is_stable(Array(c.data))
            @test !is_stable(c′)
            @test !is_stable(c′.data)
            @test !is_stable(Array(c′))
            @test !is_stable(Array(c′.data))

            c[1,1] = -c[1,1]
            @test !is_stable(c)
            c[1,1] = c[2,2]
            c[1,3] = c[1,2] = c[2,3] = c[1,1]
            @test !is_stable(c)
        end
        @testset "Catch wrong exception" begin
            let c = CIJ.ol()[1]
                @test_logs (:warn, "matrix not symmetrical: using upper half only"
                    ) is_stable(rand(6,6))
                @test_throws ArgumentError is_stable(rand(5,5))
            end
        end
    end

    @testset "Is 6×6" begin
        @test CIJ.is_6x6(CIJ.ol()[1])
        @test CIJ.is_6x6(rand(6, 6))
        @test !CIJ.is_6x6(rand(5, 5))
    end

    @testset "Is symm" begin
        @test CIJ.is_symm(CIJ.ol()[1])
        @test CIJ.is_symm(Symmetric(rand(6, 6)))
        @test !CIJ.is_symm((x = Array(CIJ.ol()[1]); x[1,2] = -x[1,2]; x))
    end
end


text:
using Turing, Distributions, Test
using MacroTools

# unit test model macro
expr = Turing.generate_observe(:x, :y)
@test expr.head == :escape
@test expr.args[1].head == :block
@test :(vi.logp += Turing.observe(sampler, y, x, vi)) in expr.args[1].args

@model testmodel_comp(x, y) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x ~ Normal(m, sqrt(s))
    y ~ Normal(m, sqrt(s))

    return x, y
end

testmodel_comp(1.0, 1.2)
c = deepcopy(Turing._compiler_)

alias1 = Dict(
            :name => :testmodel_comp_model,
            :args => [:(vi::Turing.VarInfo)],
            :kwargs => [],
            :body => :(return testmodel_comp_model(vi, Turing.SampleFromPrior()))
           )
@test c[:alias1] == MacroTools.combinedef(alias1)

alias2 = Dict(
            :name => :testmodel_comp_model,
            :args => [:(sampler::Turing.AnySampler)],
            :kwargs => [],
            :body => :(return testmodel_comp_model(Turing.VarInfo(), Turing.SampleFromPrior()))
           )
@test c[:alias2] == MacroTools.combinedef(alias2)

alias3 = Dict(
            :name => :testmodel_comp_model,
            :args => [],
            :kwargs => [],
            :body => :(return testmodel_comp_model(Turing.VarInfo(), Turing.SampleFromPrior()))
           )
@test c[:alias3] == MacroTools.combinedef(alias3)
@test length(c[:closure].args[2].args[2].args) == 6
@test mapreduce(line -> line.head == :macrocall, +, c[:closure].args[2].args[2].args) == 4

# check if drawing from the prior works
@model testmodel0(x) = begin
    x ~ Normal()
    return x
end
f0_mm = testmodel0()
@test mean(f0_mm() for _ in 1:1000) ≈ 0. atol=0.1

@model testmodel01(x) = begin
    x ~ Bernoulli(0.5)
    return x
end
f01_mm = testmodel01()
@test mean(f01_mm() for _ in 1:1000) ≈ 0.5 atol=0.1

# test if we get the correct return values
@model testmodel1(x1, x2) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x1 ~ Normal(m, sqrt(s))
    x2 ~ Normal(m, sqrt(s))

    return x1, x2
end
f1_mm = testmodel1(1., 10.)
@test f1_mm() == (1, 10)


text:
using HDF5
using DocOpt
using HPAT
#HPAT.DistributedPass.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#using ParallelAccelerator
#ParallelAccelerator.CGen.set_debug_level(3)
#ParallelAccelerator.DomainIR.set_debug_level(3)

@acc hpat function generate_file(hdf5_file, txt_file, N)
    D = 10
    A = randn(D,N)
    w = rand(1,D)
    print(w)
    S = (w*A .> .0)+0.0
    B = ParallelAccelerator.API.reshape(S, (N,))
    T = rand(N) .> .95
    Y = map((t,b)-> t?1.0-b:b, T, B)
    #Y = ones(N).+10.0
    #@par for i in 1:N
    #    Y[i] = T[i] ? B[i] : 1-B[i]
    #end
    DataSink(A, HDF5,"/points", hdf5_file)
    DataSink(Y, HDF5,"/responses", hdf5_file)
    points = [A; ParallelAccelerator.API.reshape(Y,1,N)]
    DataSink(points, TXT, txt_file)
end

function main()
    doc = """generate_logistic_regression.jl

generate random input for logistic regression example.

Usage:
  generate_logistic_regression.jl -h | --help
  generate_logistic_regression.jl [--instances=<instances>] [--path=<path>]

Options:
  -h --help                  Show this screen.
  --instances=<instances>    Specify number of instances; defaults to 2000000.
  --path=<path>              Specify output path for generated files; defaults to HPAT's default data path.
"""

    arguments = docopt(doc)

    if (arguments["--instances"] != nothing)
        instances = parse(Int, arguments["--instances"])
    else
        instances = 2*10^3
    end

    if (arguments["--path"] != nothing)
        path = arguments["--path"]
    else
        path = HPAT.getDefaultDataPath()
    end
    hdf5_file = path*"logistic_regression.hdf5"
    if isfile(hdf5_file)
        rm(hdf5_file)
    end
    txt_file = path*"logistic_regression.csv"
    if isfile(txt_file)
        rm(txt_file)
    end

    generate_file(hdf5_file, txt_file, instances)
end

main()

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
@testset "piecewisepolynomial" begin
    D = 2
    v1 = rand(D)
    v2 = rand(D)
    maha = Matrix{Float64}(I, D, D)
    v = 3
    k = PiecewisePolynomialKernel{v}(maha)

    k2 = PiecewisePolynomialKernel(v=v, maha=maha)

    @test k2(v1, v2) ≈ k(v1, v2) atol=1e-5

    @test_throws ErrorException PiecewisePolynomialKernel{4}(maha)

    @test repr(k) == "Piecewise Polynomial Kernel (v = $(v), size(maha) = $(size(maha)))"

    # Standardised tests.
    TestUtils.test_interface(k, ColVecs{Float64}; dim_in=2)
    TestUtils.test_interface(k, RowVecs{Float64}; dim_in=2)
    # test_ADs(maha-> PiecewisePolynomialKernel(v=2, maha = maha), maha)
    @test_broken "Nothing passes (problem with Mahalanobis distance in Distances)"

    test_params(k, (maha,))
end


text: square(x) = x^2


text:
# const Scalar = Union{Number, StaticArrays.Scalar, SVector{1}}

jacobian(f, w::Number, n) = ForwardDiff.derivative(v -> f(v, n), w)
jacobian(f, w::AbstractVector, n) = ForwardDiff.jacobian(v -> f(v, n), w)

eigvals(f, w, n) = real.(LinearAlgebra.eigen(jacobian(f, w, n)).values)
function eigen(f, w, n)
    λc, Rc = LinearAlgebra.eigen(jacobian(f, w, n))
    return real.(λc), real.(Rc)
end

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

struct LinearAdvectionFlux{V}
    velocity::V
end
LinearAdvectionFlux(v::AbstractVector) = LinearAdvectionFlux{SVector{length(v), eltype(v)}}(SVector(v...))

(f::LinearAdvectionFlux)(w, n) = w * (f.velocity' * n)
jacobian(f::LinearAdvectionFlux, w, n) = f.velocity' * n
eigvals(f::LinearAdvectionFlux, w, n) = f.velocity' * n

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

struct AdvectionFlux{D, VF}
    velocity_at_face::VF
end

AdvectionFlux{D}(vf) where D = AdvectionFlux{D, typeof(vf)}(vf)

function numerical_flux(flux::AdvectionFlux, mesh, w, scheme, i_face, dt)
    numerical_flux(LinearAdvectionFlux(flux.velocity_at_face(mesh, i_face)), mesh, w, scheme, i_face, dt)
end

function RotationFlux(center)
    function rotation_velocity(mesh, i_face)
        x = face_center(mesh, i_face) .- center
        return SVector{2, Float64}(-x[2], x[1])
    end
    return AdvectionFlux{2}(rotation_velocity)
end

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

struct ShallowWater{T}
    g::T
end
ShallowWater(; g=9.81) = ShallowWater{typeof(g)}(g)

function (f::ShallowWater)(v::SVector{2}, n::Union{Number, Scalar, SVector{1}})
    h, hu = v[1], v[2]
    SVector(hu * n, hu^2/h + h^2 * f.g/2) * n
end

function (f::ShallowWater)(v::SVector{3}, n::SVector{2})
    h, ux, uy = v[1], v[2]/v[1], v[3]/v[1]
    un = ux*n[1] + uy*n[2]
    return SVector(h*un, h*ux*un + h^2*n[1]*f.g/2, h*uy*un + h^2*n[2]*f.g/2)
end

function jacobian(f::ShallowWater, v::AbstractVector, n::SVector{2})
    h, ux, uy = v[1], v[2]/v[1], v[3]/v[1]
    un = ux * n[1] + uy * n[2]
    g = f.g
    @SMatrix [0                 n[1]        n[2];
              -ux*un+g*h*n[1]   un+ux*n[1]  ux*n[2];
              -uy*un+g*h*n[2]   uy*n[1]     un+uy*n[2]]
end

function eigvals(f::ShallowWater, v, n::SVector{2})
    h, ux, uy = v[1], v[2]/v[1], v[3]/v[1]
    un = ux * n[1] + uy * n[2]
    g = f.g
    c = sqrt(g*h)
    return SVector(un - c, un, un + c)
end

function eigen(f::ShallowWater, v, n::SVector{2})
    h, ux, uy = v[1], v[2]/v[1], v[3]/v[1]
    un = ux * n[1] + uy * n[2]
    ut = -ux * n[2] + uy * n[1]
    g = f.g
    c = sqrt(g*h)
    vals = SVector(un - c, un, un + c)
    vects = 0.5 .* @SMatrix [1.0                    0.0   1.0;
                             n[1]*(un-c) + n[2]*ut  n[2]  n[1]*(c+un) + n[2]*ut;
                             n[1]*ut + n[2]*(un-c)  n[1]  n[1]*ut + n[2]*(c+un)]
    return vals, vects
end


text:
import ..update!, ..select_action, ..reset!

using StaticArrays

using Ramnet.Models.AltDiscriminators:RegularizedFunctionalDiscriminator
using Ramnet.Optimizers
import Ramnet.Loss

struct RegularizedBinaryPolicyLoss <: Loss.AbstractLoss end

# Policy approximator with only two admissible actions: {-1, 1}.
struct RegularizedBinaryPolicy{OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    opt::FunctionalOptimizer{RegularizedBinaryPolicyLoss}
f::RegularizedFunctionalDiscriminator{1,T,E}
    rng::MersenneTwister
end

function Loss.grad(::RegularizedBinaryPolicyLoss, p::RegularizedBinaryPolicy{OS,T,O,E}, observation::O, action::Int, G) where {OS,T,O,E}
    sign(action) * (1 - probability(p, action, observation)) * G
end

function RegularizedBinaryPolicy(::Type{O}, n::Int, encoder::E; α::Float64=1.0, η::Float64=0.1, epochs::Int=1, partitioner::Symbol=:uniform_random, seed::Union{Nothing,Int}=nothing) where {OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    RegularizedBinaryPolicy{OS,T,O,E}(
        FunctionalOptimizer(RegularizedBinaryPolicyLoss(); learning_rate=η, epochs),
        RegularizedFunctionalDiscriminator{1}(OS, n, α, encoder, partitioner; seed),
        isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)
    )
end

# Compute π(action | observation)
function probability(p::RegularizedBinaryPolicy{OS,T,O,E}, action::Int, observation::O) where {OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    prob = 1 / (1 + exp(predict(p.f, observation)[1]))

    return action == 1 ? prob : 1 - prob
end

# TODO: Validate action (must be ether -1 or +1)
function update!(p::RegularizedBinaryPolicy{OS,T,O,E}, observation::O, action::Int, G) where {OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    gradient = Optimizers.grad(p.opt, p, observation, action, G)

    train!(p.f, observation, learning_rate(p.opt) * gradient)

    return nothing
end

function select_action(p::RegularizedBinaryPolicy{OS,T,O,E}, observation::O) where {OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    return rand(p.rng) < probability(p, 1, observation) ? 1 : -1
    end

function reset!(p::RegularizedBinaryPolicy{OS,T,O,E}) where {OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    Ramnet.reset!(p.f)

    return nothing
end


text:
using ComradeSoss
using Test

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


text:
# computes strongly connected components of a control flow graph `cfg`
# NOTE adapted from https://github.com/JuliaGraphs/LightGraphs.jl/blob/6134e21484f2312b33ba753d47dc2fc4d6b35fac/src/connectivity.jl#L221-L311
# since to load an entire LightGraphs.jl is a bit cost-ineffective in terms of a trade-off of latency vs. maintainability
function strongly_connected_components(g::CFG)
    T = Int
    zero_t = zero(T)
    one_t = one(T)
    nvg = nv(g)
    count = one_t

    index = zeros(T, nvg)         # first time in which vertex is discovered
    stack = Vector{T}()           # stores vertices which have been discovered and not yet assigned to any component
    onstack = zeros(Bool, nvg)    # false if a vertex is waiting in the stack to receive a component assignment
    lowlink = zeros(T, nvg)       # lowest index vertex that it can reach through back edge (index array not vertex id number)
    parents = zeros(T, nvg)       # parent of every vertex in dfs
    components = Vector{Vector{T}}()    # maintains a list of scc (order is not guaranteed in API)

    dfs_stack = Vector{T}()

    @inbounds for s in vertices(g)
        if index[s] == zero_t
            index[s] = count
            lowlink[s] = count
            onstack[s] = true
            parents[s] = s
            push!(stack, s)
            count = count + one_t

            # start dfs from 's'
            push!(dfs_stack, s)

            while !isempty(dfs_stack)
                v = dfs_stack[end] #end is the most recently added item
                u = zero_t
                @inbounds for v_neighbor in outneighbors(g, v)
                    if index[v_neighbor] == zero_t
                        # unvisited neighbor found
                        u = v_neighbor
                        break
                        #GOTO A push u onto DFS stack and continue DFS
                    elseif onstack[v_neighbor]
                        # we have already seen n, but can update the lowlink of v
                        # which has the effect of possibly keeping v on the stack until n is ready to pop.
                        # update lowest index 'v' can reach through out neighbors
                        lowlink[v] = min(lowlink[v], index[v_neighbor])
                    end
                end
                if u == zero_t
                    # All out neighbors already visited or no out neighbors
                    # we have fully explored the DFS tree from v.
                    # time to start popping.
                    popped = pop!(dfs_stack)
                    lowlink[parents[popped]] = min(lowlink[parents[popped]], lowlink[popped])

                    if index[v] == lowlink[v]
                        # found a cycle in a completed dfs tree.
                        component = Vector{T}()

                        while !isempty(stack) #break when popped == v
                            # drain stack until we see v.
                            # everything on the stack until we see v is in the SCC rooted at v.
                            popped = pop!(stack)
                            push!(component, popped)
                            onstack[popped] = false
                            # popped has been assigned a component, so we will never see it again.
                            if popped == v
                                # we have drained the stack of an entire component.
                                break
                            end
                        end

                        reverse!(component)
                        push!(components, component)
                    end

                else #LABEL A
                    # add unvisited neighbor to dfs
                    index[u] = count
                    lowlink[u] = count
                    onstack[u] = true
                    parents[u] = v
                    count = count + one_t

                    push!(stack, u)
                    push!(dfs_stack, u)
                    # next iteration of while loop will expand the DFS tree from u.
                end
            end
        end
    end

    # # assert with the original implementation
    # components′ = LightGraphs.strongly_connected_components(cfg_to_sdg(g))
    # @assert Set(Set.(components)) == Set(Set.(components′))

    return components
end

# compatibility with LightGraphs.jl interfaces
@inline nv(cfg::CFG) = length(cfg.blocks)
@inline vertices(cfg::CFG) = 1:nv(cfg)
@inline outneighbors(cfg::CFG, v) = cfg.blocks[v].succs

# import LightGraphs: SimpleDiGraph, add_edge!, LightGraphs
# function cfg_to_sdg(cfg::CFG)
#     g = SimpleDiGraph(length(cfg.blocks))
#     for (v, block) in enumerate(cfg.blocks)
#         for succ in block.succs
#             add_edge!(g, v, succ)
#         end
#     end
#     return g
# end


text:
export ASTTools

"""
    ASTTools

Tools for manipulating Julia ASTs. This module provides similar functionality
with MacroTools, but based on MLStyle.
"""
module ASTTools

using MLStyle
export @capture

# capture
function capturing_analysis(expr, out, is_literal)

    @match expr begin
        ::QuoteNode => capturing_analysis(expr.value, out, true)

        if is_literal
        end && Expr(:$, args...) => foreach(x -> capturing_analysis(x, out, false), args)

        if is_literal
        end && Expr(_, args...) => foreach(x -> capturing_analysis(x, out, true), args)

        if is_literal
        end && _ => nothing

        # not literal
        ::Symbol => (push!(out, expr); nothing)

        Expr(:quote, args...) => foreach(x -> capturing_analysis(x, out, true), args)

        :(Do($(args...))) => foreach(args) do arg
            @match arg begin
                Expr(:kw, key::Symbol, value) => begin
                    push!(out, key)
                end
                _ => nothing
            end
        end
        :($a || $b) => let out1 = Set(Symbol[]), out2 = Set(Symbol[])
            capturing_analysis(a, out1, false)
            capturing_analysis(b, out2, false)

            union!(out, intersect(out1, out2))
            nothing
        end
        # type pattern
        Expr(:(::), a, _) => capturing_analysis(a, out, false)
        # dict pattern
        :(Dict($(args...))) => foreach(args) do arg
            @match arg begin
                :($_ => $v) => capturing_analysis(v, out, false)
                _ => nothing
            end
        end
        # app pattern
        :($_($(args...))) => foreach(x -> capturing_analysis(x, out, false), args)
        # other expr
        Expr(_, args...) => foreach(x -> capturing_analysis(x, out, false), args)
        # ref pattern
        Expr(:&, _) ||
            # predicate
            Expr(:function, _...) ||
            Expr(:if, _...) ||
            x => nothing
    end
end

function capture(template)
    syms = Set(Symbol[])
    capturing_analysis(template, syms, true)
    out_expr =
        Expr(:call, Dict, (Expr(:call, =>, QuoteNode(each), each) for each in syms)...)
    arg_sym = gensym()
    let template = Expr(:quote, template)
        quote
            $arg_sym -> $MLStyle.@match $arg_sym begin
                $(template) => $(out_expr)
            end
        end
    end
end

"""
    @capture template
    @capture template expr

Template matching for expressions.

```julia
julia> @capture f(\$x)  :(f(2))
Dict{Symbol,Int64} with 1 entry:
  :x => 2
```
"""
:(@capture)

macro capture(template)
    capture(template) |> esc
end

macro capture(template, ex)
    Expr(:call, capture(template), ex) |> esc
end

end


text:
const TRY_BUT_ALLOW_FAILURES_URL_LIST = String[
    "https://gitlab.com/UnofficialJuliaMirror/GeoEfficiency.jl-b1f22607-7830-5ceb-9a8e-05ab1ac77008",
    "https://gitlab.com/UnofficialJuliaMirror/JuliaLang-julia",
    "https://gitlab.com/UnofficialJuliaMirror/julia.jl-1222c4b2-2114-5bfd-aeef-88e4692bbb3e",
    "https://gitlab.com/UnofficialJuliaMirror/libuv-libuv",
    ]


text:
struct AggregatedResults <: Results
    variables::Dict{Symbol, DataFrames.DataFrame}
    total_cost::Dict
    optimizer_log::Dict
    time_stamp::DataFrames.DataFrame
    duals::Dict{Symbol, Any}
end

get_res_variables(result::AggregatedResults) = result.variables
get_cost(result::AggregatedResults) = result.total_cost
get_log(result::AggregatedResults) = result.optimizer_log
get_time_stamp(result::AggregatedResults) = result.time_stamp
get_duals(result::AggregatedResults) = result.duals
get_variables(result::AggregatedResults) = result.dual_variables

# This function returns the length of time_stamp for each step
# that is unique to that step and not overlapping with the next time step.
function _count_time_overlap(stage::String,
    step::Array,
    variable::Array,
    references::Dict{Any,Any})
    date_df = references[stage][variable[1]]
    step_df = DataFrames.DataFrame(Date = Dates.DateTime[],
            Step = String[],
            File_Path = String[])
    for n in 1:length(step)
        step_df = vcat(step_df,date_df[date_df.Step .== step[n], :])
    end
    ref = DataFrames.DataFrame()
    for (ix,time) in enumerate(step_df.Date)
        file_path =  step_df[ix, :File_Path]
        time_file_path = joinpath(dirname(file_path), "time_stamp.feather")
        time_stamp = DataFrames.DataFrame(Feather.read("$time_file_path"))
        time_stamp = shorten_time_stamp(time_stamp)
        append!(ref,time_stamp)
    end
    if size(unique(ref),2) == size(ref,2)
        extra_time_length = 0
    else
        extra_time_length = size(unique(ref),1)./(length(step)+1)
    end
    return extra_time_length
end
# This is the count_time_overlap for if all results for a stage are desired

function _count_time_overlap(stage::String, references::Dict{Any,Any})
    variable = collect(keys(references[stage]))
    date_df = references[stage][variable[1]]
    ref = DataFrames.DataFrame()
    for (ix,time) in enumerate(date_df.Date)
        file_path =  date_df[ix, :File_Path]
        time_file_path = joinpath(dirname(file_path), "time_stamp.feather")
        time_stamp = DataFrames.DataFrame(Feather.read("$time_file_path"))
        time_stamp = shorten_time_stamp(time_stamp)
        append!(ref,time_stamp)
    end
    if size(unique(ref),2) == size(ref,2)
    extra_time_length = 0
    else
    extra_time_length = size(unique(ref),1)./(length(step)+1)
    end
    return extra_time_length
end

function _reading_references(results::Dict, duals::Array, stage::String, step::Array,
                             references::Dict, extra_time::Int64)

    for name in (dual)
        date_df = references[stage][name]
        step_df = DataFrames.DataFrame(Date = Dates.DateTime[], Step = String[], File_Path = String[])
        for n in 1:length(step)
            step_df = vcat(step_df,date_df[date_df.Step .== step[n], :])
        end
        results[name] = DataFrames.DataFrame()
        for (ix,time) in enumerate(step_df.Date)
            file_path = step_df[ix, :File_Path]
            var = Feather.read("$file_path")
            correct_var_length = size(1:(size(var,1) - extra_time),1)
            results[name] = vcat(results[name],var[1:correct_var_length,:])
        end
    end
    return results
end

function _reading_references(results::Dict, dual::Array, stage::String,
                             references::Dict, extra_time::Int64)
    for name in dual
        date_df = references[stage][name]
        results[name] = DataFrames.DataFrame()
        for (ix,time) in enumerate(date_df.Date)
            file_path = date_df[ix, :File_Path]
            var = Feather.read(file_path)
            correct_var_length = size(1:(size(var,1) - extra_time_length),1)
            results[name] = vcat(results[name],var[1:correct_var_length,:])
        end
    end
    return results
end

function _removing_extra_time(file_path::String, extra_time_length::Number)
    time_stamp = []
    time_file_path = joinpath(dirname(file_path), "time_stamp.feather")
    temp_time_stamp = DataFrames.DataFrame(Feather.read("$time_file_path"))
    non_overlap = size((1:size(temp_time_stamp, 1) - extra_time_length - 1),1)
    time_stamp = vcat(time_stamp,temp_time_stamp[(1:non_overlap),:])
    return time_stamp
end

""" This sums all of the rows in a result dataframe """
function rowsum(variable::DataFrames.DataFrame, name::String)
    variable = DataFrames.DataFrame(Symbol(name) => sum.(eachcol(variable)))
    return variable
end
""" This sums each column in a result dataframe """
function columnsum(variable::DataFrames.DataFrame)
    shortvar = DataFrames.DataFrame()
    varnames = collect(names(variable))
    eachsum = (sum.(eachrow(variable)))
    for i in 1: size(variable,1)
        df = DataFrames.DataFrame(Symbol(varnames[i]) => eachsum[i])
        shortvar = hcat(shortvar,df)
    end
    return shortvar
end

function _find_duals(variables::Array)
    duals = []
    for i in 1:length(variables)
        if occursin("dual", String.(variables[i]))
            duals = vcat(duals, variables[i])
        end
    end
    return duals
end

"""
    load_simulation_results(stage, step,variable,references)

This function goes through the reference table of file paths and
aggregates the results over time into a struct of type OperationsProblemResults
for the desired step range and variables

# Arguments
-`stage::String = "stage-1"``: The stage of the results getting parsed, stage-1 or stage-2
-`step::Array{String} = ["step-1", "step-2", "step-3"]`: the steps of the results getting parsed
-`variable::Array{Symbol} = [:P_ThermalStandard, :P_RenewableDispatch]`: the variables to be parsed
-`references::Dict = ref`: the reference dictionary created in execute!
or with make_references

# Example
```julia
stage = "stage-1"
step = ["step-1","step-2", "step-3"] # has to match the date range
variable = [:P_ThermalStandard, :P_RenewableDispatch]
results = load_simulation_results(stage,step, variable, references)
```
# Accepted Key Words
-`write::Bool`: if true, the aggregated results get written back to the results
file in the folder structure
"""
function load_simulation_results(stage::String,
         step::Array,
         variable::Array,
         references::Dict{Any,Any}; kwargs...)
    variables = Dict()
    duals = Dict()
    dual = _find_duals(collect(keys(references[stage])))
    time_stamp = DataFrames.DataFrame(Range = Dates.DateTime[])
    extra_time_length = _count_time_overlap(stage, step, variable, references)
    for l in 1:length(variable)
        date_df = references[stage][variable[l]]
        step_df = DataFrames.DataFrame(Date = Dates.DateTime[], Step = String[], File_Path = String[])
        for n in 1:length(step)
            step_df = vcat(step_df,date_df[date_df.Step .== step[n], :])
        end
        variables[(variable[l])] = DataFrames.DataFrame()
        for (ix,time) in enumerate(step_df.Date)
            file_path = step_df[ix, :File_Path]
            var = Feather.read("$file_path")
            correct_var_length = size(1:(size(var,1) - extra_time_length),1)
            variables[(variable[l])] = vcat(variables[(variable[l])],var[1:correct_var_length,:])
            if l == 1
                time_stamp = DataFrames.DataFrame(Range = _removing_extra_time(file_path, extra_time_length))
            end
        end
    end
    file_path = references[stage][variable[1]][1,:File_Path]
    optimizer = optimizer = JSON.parse(open(joinpath(dirname(file_path), "optimizer_log.json")))
    obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => optimizer["obj_value"])

    if !isempty(dual)
        duals = _reading_references(duals, dual, stage, step, references, extra_time_length)
        results = make_results(variables, obj_value, optimizer, time_stamp, duals)
    else
        results = make_results(variables, obj_value, optimizer, time_stamp)
    end
    if (:write in keys(kwargs))
        write_model_results(results, normpath("$file_path/../../../../"),"results")
    end
    return results
end

"""
    load_simulation_results(stage, references)

This function goes through the reference table of file paths and
aggregates the results over time into a struct of type OperationsProblemResults

# Arguments
-`stage::String = "stage-1"`: The stage of the results getting parsed, stage-1 or stage-2
-`references::Dict = ref`: the reference dictionary created in execute!
or with make_references

# Example
```julia
stage = "stage-1"
step = ["step-1","step-2", "step-3"] # has to match the date range
variable = [:P_ThermalStandard, :P_RenewableDispatch]
results = load_simulation_results(stage,step, variable, references)
```
# Accepted Key Words
-`write::Bool`: if true, the aggregated results get written back to the results
file in the folder structure
"""

function load_simulation_results(stage::String, references::Dict{Any,Any}; kwargs...)

    variables = Dict()
    duals = Dict()
    variable = (collect(keys(references[stage])))
    dual = _find_duals(variable)
    variable = setdiff(variable,duals)

    time_stamp = DataFrames.DataFrame(Range = Dates.DateTime[])
    extra_time_length = _count_time_overlap(stage,references)
    for l in 1:length(variable)
        date_df = references[stage][variable[l]]
        variables[(variable[l])] = DataFrames.DataFrame()
        for (ix,time) in enumerate(date_df.Date)
            file_path = date_df[ix, :File_Path]
            var = Feather.read(file_path)
            correct_var_length = size(1:(size(var,1) - extra_time_length),1)
            variables[(variable[l])] = vcat(variables[(variable[l])],var[1:correct_var_length,:])
            if l == 1
                time_stamp = DataFrames.DataFrame(Range = _removing_extra_time(file_path, extra_time_length))
            end
        end
    end

    file_path = references[stage][variable[1]][1,:File_Path]
    optimizer = optimizer = JSON.parse(open(joinpath(dirname(file_path), "optimizer_log.json")))
    obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => optimizer["obj_value"])

    if !isnothing(dual)
        duals = _reading_references(duals, dual, stage, references, extra_time_length)
        results = make_results(variables, obj_value, optimizer, time_stamp, duals)
    else
        results = make_results(variables, obj_value, optimizer, time_stamp)
    end
    if (:write in keys(kwargs)) == true
        file_type = get(kwargs, :file_type, Feather)
        write_model_results(results, normpath("$file_path/../../../../"),"results"; file_type = file_type)
    end
    return results
end

function _concat_string(duals::Vector{Symbol})
    duals = (String.(duals)).*"_dual"
    return duals
end
"""
    make_references(sim::Simulation, date_run::String)

Creates a dictionary of variables with a dictionary of stages
that contains dataframes of date/step/and desired file path
so that the results can be parsed sequentially by variable
and stage type.

**Note:** make_references can only be run after run_sim_model
or else, the folder structure will not yet be populated with results

# Arguments
-`sim::Simulation = sim`: simulation object created by Simulation()
-`date_run::String = "2019-10-03T09-18-00"``: the name of the file created
that contains the specific simulation run of the date run and "-test"

# Example
```julia
sim = Simulation("test", 7, stages, "/Users/yourusername/Desktop/";
verbose = true, system_to_file = false)
execute!(sim::Simulation; verbose::Bool = false, kwargs...)
references = make_references(sim, "2019-10-03T09-18-00-test")
```
"""
function make_references(sim::Simulation, date_run::String; kwargs...)
    sim.ref.date_ref[1] = sim.daterange[1]
    sim.ref.date_ref[2] = sim.daterange[1]
    references = Dict()
    for (ix, stage) in enumerate(sim.stages)
        variables = Dict{Symbol, Any}()
        interval = PSY.get_forecasts_interval(stage.sys)
        variable_names = (collect(keys(sim.stages[ix].canonical.variables)))
        if :dual_constraints in keys(kwargs) && !isnothing(get_constraints(stage.canonical))
            dual_cons = _concat_string(kwargs[:dual_constraint])
            variable_names = vcat(variable_names, dual_cons)
        end
        for name in variable_names
            variables[name] = DataFrames.DataFrame(Date = Dates.DateTime[],
                                           Step = String[], File_Path = String[])
        end
        for s in 1:(sim.steps)
            for run in 1:stage.executions
                sim.ref.current_time = sim.ref.date_ref[ix]
                for name in variable_names
                    full_path = joinpath(sim.ref.raw, "step-$(s)-stage-$(ix)",
                                replace_chars("$(sim.ref.current_time)",":","-"), "$(name).feather")
                    if isfile(full_path)
                        date_df = DataFrames.DataFrame(Date = sim.ref.current_time,
                                                       Step = "step-$(s)", File_Path = full_path)
                        variables[name] = vcat(variables[name], date_df)
                    else
                        println("$full_path, no such file path")
                     end
                end
                sim.ref.run_count[s][ix] += 1
                sim.ref.date_ref[ix] = sim.ref.date_ref[ix] + interval
            end
        end
        references["stage-$ix"] = variables
        stage.execution_count = 0
    end
    return references
end


text:
using IterativeClosestPoint
using Documenter

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

makedocs(;
    modules=[IterativeClosestPoint],
    authors="Tom McGrath <tmcgrath325@gmail.com> and contributors",
    repo="https://github.com/tmcgrath325/IterativeClosestPoint.jl/blob/{commit}{path}#{line}",
    sitename="IterativeClosestPoint.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tmcgrath325.github.io/IterativeClosestPoint.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/tmcgrath325/IterativeClosestPoint.jl",
    devbranch="main",
)


text:
# see correct result with slanted (narrow) L1 modes at y=+-5.  L1 should not be round blobs
# https://github.com/JuliaRobotics/RoME.jl/pull/434#issuecomment-817246038

using RoME
using Test
using TensorCast

##

@testset "Test correlation induced by partials" begin

##

# start with an empty factor graph object
N = 200
fg = initfg()
getSolverParams(fg).N = N

getSolverParams(fg).useMsgLikelihoods = true

addVariable!(fg, :x0, Pose2)
addVariable!(fg, :l1, Point2)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0; 0; 0], diagm([0.3;0.3;0.3].^2)) ))
ppr = Pose2Point2Range(MvNormal([7.3], diagm([0.3].^2)))
addFactor!(fg, [:x0; :l1], ppr )

addVariable!(fg, :x1, Pose2)
pp = Pose2Pose2(MvNormal([9.8;0;0.8], diagm([0.3;0.3;0.05].^2)))
ppr = Pose2Point2Range(MvNormal([6.78], diagm([0.3].^2)))
addFactor!(fg, [:x0; :x1], pp )
addFactor!(fg, [:x1; :l1], ppr )

##

tree = solveTree!(fg)

## check that stuff is where it should be

L1 = getPoints(getBelief(fg, :l1))
@cast L1_[j,i] := L1[i][j]

@test 0.3*N < sum(L1_[2,:] .< 0)
@test 0.3*N < sum(0 .< L1_[2,:])

L1_n = L1_[:, L1_[2,:] .< 0]
L1_p = L1_[:, 0 .< L1_[2,:]]

mvn = fit(MvNormal, L1_n)
mvp = fit(MvNormal, L1_p)

# check diagonal structure for correlation
@test isapprox(mvn.Σ.mat[1,1], 2.5, atol=2.2)
@test isapprox(mvn.Σ.mat[2,2], 2.5, atol=2.2)
@test isapprox(mvn.Σ.mat[1,2], 2.0, atol=1.6)

@test isapprox(mvp.Σ.mat[1,1], 2.5, atol=2.2)
@test isapprox(mvp.Σ.mat[2,2], 2.5, atol=2.2)
@test isapprox(mvp.Σ.mat[1,2], -2.0, atol=1.6)

# sanity check for symmetry
@test mvn.Σ.mat - mvn.Σ.mat' |> norm < 0.01

# test means in the right location
@test isapprox(mvn.μ[1], 5.4, atol=2.0)
@test isapprox(mvn.μ[2], -4.2, atol=2.0)

@test isapprox(mvp.μ[1], 5.4, atol=2.0)
@test isapprox(mvp.μ[2], 4.2, atol=2.0)

##

end

##

# using Cairo, RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)

##

# plotSLAM2D(fg, drawPoints=false)

# plotKDE(fg, :l1, levels=10)

##


text:
module MLFS
using MLStyle
using Setfield
using DataStructures: list, cons, nil, LinkedList, Cons, Nil
include("TypedIO.jl")
import .TypedIO: fromVec, toVec, @typedIO

"""
forward declaration of IR.ExprImpl
"""
abstract type ExprImpl end

struct CFunc{A, B}
    f :: Function
end

_targs = [Symbol(:a_, i) for i = 1:3]
_args = [Symbol(:t_, i) for i = 1:3]
function (f::CFunc{R, Tuple{}})()::R where R
    f.f()
end

for i in 1:3
    let args = _args[1:i],
        targs = _targs[1:i],
        params = [:($a :: $t) for (a, t) in zip(args, targs)]

        Base.eval(MLFS, quote
            function (f::CFunc{R, Tuple{$(targs...)}})($(params...))::R where {R, $(targs...)}
                f.f($(args...))
            end
        end)
    end
end

include("HM.jl")
using .HM
export T
@active T(x) begin
    @match x begin
        App(Nom(:Type), x) => Some(x)
        _ => nothing
    end
end
(::Type{T})(x::HMT) = App(Nom(:Type), x)

include("Store.jl")
include("FS.jl")
include("Core.jl")
include("Signals.jl")
include("AST/Surf.jl")
include("AST/IR.jl")
include("AST/SurfBuild.jl")
include("TypeClass.jl")
include("Infer.jl")

include("TypeErasure.jl")
include("Compiler/ToJulia.jl")
include("Modular.jl")
include("IR2Jl.jl")
# Write your package code here.

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AdminAppsRestrictedApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminAppsRestrictedList(_api::AdminAppsRestrictedApi; limit=nothing, cursor=nothing, team_id=nothing, enterprise_id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.apps.restricted.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.query, "enterprise_id", enterprise_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

List restricted apps for an org or workspace.
Param: limit::Int32
Param: cursor::String
Param: team_id::String
Param: enterprise_id::String
Return: Dict{String, Any}
"""
function adminAppsRestrictedList(_api::AdminAppsRestrictedApi; limit=nothing, cursor=nothing, team_id=nothing, enterprise_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminAppsRestrictedList(_api; limit=limit, cursor=cursor, team_id=team_id, enterprise_id=enterprise_id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminAppsRestrictedList(_api::AdminAppsRestrictedApi, response_stream::Channel; limit=nothing, cursor=nothing, team_id=nothing, enterprise_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminAppsRestrictedList(_api; limit=limit, cursor=cursor, team_id=team_id, enterprise_id=enterprise_id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminAppsRestrictedList


text:
@testset "101.symmetric-tree.jl" begin
    @test is_symmetric_tree(TreeNode{Int}([1, 2, 2, 3, 4, 4, 3])) == true
    @test is_symmetric_tree(TreeNode{Int}([1, 2, 2, nothing, 3, nothing, 3])) == false
end


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"General Utilities."

module Utils

export arrange, usage_and_exit, bin2dec, ∑, array2num, find_neighbors, find_all_neighbors

arrange(vs :: Vector, dtype = String) = [parse(dtype, v) for v in vs]

"If `is_exit`, print the usage of the program and, then, quit."
function usage_and_exit(is_exit = true)
    if is_exit
        println("Usage: ./dayN [INPUT]")

        exit(1)
    end
end

"Convert integer binary digit array `B` to a decimal integer."
bin2dec(B :: AbstractArray) = sum((*).(B, (^).(2, (length(B) - 1):-1:0)))

"Summation of all integer numbers between [0, `n`]."
∑(n :: Int64) = (n^2 + n)/2

"Transform an array `N` of digits in a number."
array2num(N :: AbstractArray) = N .* (^).(10, (length(N) - 1):-1:0) |> sum

"Find the neighbors (non-diagonal) of `(i, j)` in the matrix `A`."
function find_neighbors(A :: AbstractMatrix, (i, j) :: Tuple{Int, Int})
    neighbors = Array{Tuple{Int, Int}}([])
    idxs = CartesianIndices(A)

    axial = [(i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)]

    for idx in axial
        if CartesianIndex(idx) in idxs
            push!(neighbors, idx)
        end
    end

    return neighbors
end

"Find all the neighbors of `(i, j)` in the matrix `A`."
function find_all_neighbors(A :: AbstractMatrix, (i, j) :: Tuple{Int, Int})
    neighbors = Array{Tuple{Int, Int}}([])
    idxs = CartesianIndices(A)

    axial = [(i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)]
    diagonal = [(i - 1, j - 1), (i - 1, j + 1), (i + 1, j - 1), (i + 1, j + 1)]

    for idx in [axial; diagonal]
        if CartesianIndex(idx) in idxs
            push!(neighbors, idx)
        end
    end

    return neighbors
end

end


text:
@testset "Limited Length Cycles" begin
    completedg = CompleteDiGraph(4)
    pathdg = PathDiGraph(5)
    cycledg = CycleDiGraph(5)

    @test length(simplecycles_limited_length(completedg, 0)) == 0
    @test length(simplecycles_limited_length(completedg, 1)) == 0
    @test length(simplecycles_limited_length(completedg, 2)) == 6
    @test length(simplecycles_limited_length(completedg, 3)) == 14
    @test length(simplecycles_limited_length(completedg, 4)) == 20
    @test length(simplecycles_limited_length(completedg, 4, 10)) == 10
    @test length(simplecycles_limited_length(completedg, 4, typemax(Int))) == 20

    @test length(simplecycles_limited_length(pathdg, 1)) == 0
    @test length(simplecycles_limited_length(pathdg, 2)) == 0
    @test length(simplecycles_limited_length(pathdg, 3)) == 0
    @test length(simplecycles_limited_length(pathdg, 4)) == 0
    @test length(simplecycles_limited_length(pathdg, 5)) == 0

    @test length(simplecycles_limited_length(cycledg, 1)) == 0
    @test length(simplecycles_limited_length(cycledg, 2)) == 0
    @test length(simplecycles_limited_length(cycledg, 3)) == 0
    @test length(simplecycles_limited_length(cycledg, 4)) == 0
    @test length(simplecycles_limited_length(cycledg, 5)) == 1

    selfloopg = DiGraph([
        0 1 0 0;
        0 0 1 0;
        1 0 1 0;
        0 0 0 1;
    ])
    cycles = simplecycles_limited_length(selfloopg, nv(selfloopg))
    @test [3] in cycles
    @test [4] in cycles
    @test [1, 2, 3] in cycles || [2, 3, 1] in cycles || [3, 1, 2] in cycles
    @test length(cycles) == 3

    octag = smallgraph(:octahedral)
    octadg = DiGraph(octag)
    octalengths, _ = simplecycleslength(octadg)
    for k = 1:6
        @test sum(octalengths[1:k]) == length(simplecycles_limited_length(octag, k))
        @test sum(octalengths[1:k]) == length(simplecycles_limited_length(octadg, k))
    end
end


text:
module MDTrajAnalyse

export AbstractTrajectory,
       Trajectory,
       natoms,
       getatoms,
       dihedral,
       distances,
       distances!,
       cellvolume,
       angletoframe,
       sphericalview,
       subtrajectory,
       compute_rdf,
       get_close_atoms,

       AbstractAtomNames,
       AtomNames,
       NoAtomName,

       read_pdb,
       read_trajectory,
       read_dcd,
       read_xyz,
       write_xyz,

       AbstractUnitCell,
       AbtractPeriodicCell,
       AbstractOrthorombicCell,
       CubicCell,
       NonPeriodic,
       OrthorombicCell,
       TriclinicCell,
       VariableCell,
       celldiag,
       cellmatrix,
       parsecell,
       volume

include("cell.jl")
include("trajectory.jl")
include("fileaccess.jl")
include("properties.jl")
include("SubModules/visualize.jl")

end # module


text:
"""
    calculateSWI(data, options)

Returns the calculated SWI using 'data' and 'options'.
"""
function calculateSWI(data, options=Options())
    if !options.writesteps
        options.writedir = nothing
    end
    getswimag(data, options) .* getswiphase(data, options)
end

function getswimag(data, options)
    combined_echoes = combine_echoes_swi(data.mag, data.TEs, options.combination_type)
    @debug savenii(combined_echoes, "combined_mag", options.writedir, data.header)

    sensitivity = options.sensitivity
    if sensitivity === nothing
        sensitivity = getsensitivity(data.mag[:,:,:,1], pixdim=getpixdim(data))
    end

    @debug savenii(sensitivity, "sensitivity", options.writedir, data.header)
    combined_echoes ./ sensitivity
end

function combine_echoes_swi(mag, TEs, type)
    T2s, factor = gettissue_easy(:B7T)

    if ndims(mag) == 3 # only one echo
        copy(mag)
    elseif type == :SNR
        RSS(mag)
    elseif typeof(type) <: Pair
        type, (w1, w2) = type
        if type == :CNR
            S(TE, tissue) = factor[tissue] * exp(-TE / T2s[tissue])
            w(TE, w1, w2) = S(TE, w1) - S(TE, w2)

            combine_weighted(mag, w.(TEs, w1, w2))
        end
    else
        throw("ERROR: $type not defined for combination of echoes!")
    end
end

function combine_weighted(mag, w)
    combined = zeros(size(mag)[1:3])
    for ieco in 1:size(mag, 4)
        combined .+= mag[:,:,:,ieco] .* w[ieco]
    end
    combined
end

function getswiphase(data, options)
    mask = getrobustmask(view(data.mag,:,:,:,1))
    @debug savenii(mask, "maskforphase", options.writedir, data.header)

    combined = getcombinedphase(data, options, mask)
    scaleandthreshold!(combined, mask, options.mode, options.level)
end

function getcombinedphase(data, options, mask)
    phase = data.phase
    mag = data.mag
    TEs = data.TEs
    σ = options.σ

    # reshape TEs to the highest dimension, eg. "size(TEs) = (1, 1, 1, 3)"
    TEs = reshape(TEs, 1, 1, 1, length(TEs))
    @show typeof(phase)
    unwrapped = similar(phase)
    if options.unwrapping == :laplacian
        for iEco in 1:size(phase, 4)
            unwrapped[:,:,:,iEco] .= laplacianunwrap(view(phase,:,:,:,iEco))
        end
        @debug savenii(unwrapped, "unwrappedphase", options.writedir, data.header)

        for iEco in 1:size(phase, 4)
            smoothed = gaussiansmooth3d(unwrapped[:,:,:,iEco], σ; mask=mask, dims=1:2)
            unwrapped[:,:,:,iEco] .-= smoothed
        end
        combined = combine_echoes(unwrapped, mag, TEs)
        @debug savenii(unwrapped, "filteredphase", options.writedir, data.header)
        @debug savenii(combined, "combinedphase", options.writedir, data.header)

    elseif unwrapping == :laplacianslice
        for iEco in 1:size(phase, 4), iSlc in 1:size(phase, 3)
            unwrapped[:,:,iSlc,iEco] .= laplacianunwrap(view(phase,:,:,iSlc,iEco))
            smoothed = gaussiansmooth3d(unwrapped[:,:,iSlc,iEco], σ; mask=mask[:,:,iSlc], dims=1:2)
            unwrapped[:,:,iSlc,iEco] .-= smoothed
        end
        combined = combine_echoes(unwrapped, mag, TEs)
        @debug savenii(unwrapped, "unwrapped", options.writedir, data.header)
        @debug savenii(combined, "combinedphase", options.writedir, data.header)

    elseif unwrapping == :romeo
        unwrapped = kunwrap(phase, mag; TEs = TEs)#, mask = mask)
        @debug savenii(unwrapped, "unwrappedphase", options.writedir, data.header)

        combined = combine_echoes(unwrapped, mag, TEs)
        @debug savenii(combined, "combinedphase", options.writedir, data.header)

        combined .-= gaussiansmooth3d(combined, σ; mask = mask, dims = 1:2)
        @debug savenii(combined, "filteredphase", options.writedir, data.header)

    else
        error("Unwrapping $unwrapping not defined!")
    end

    combined
end

combine_echoes(unwrapped::AbstractArray{T,3}, mag, TEs) where T = copy(unwrapped) # one echo
function combine_echoes(unwrapped::AbstractArray{T,4}, mag, TEs) where T
    dim = 4
    TEs = reshape(TEs, ones(Int, dim-1)..., length(TEs)) # size = (1,1,1,nEco)

    combined = sum(unwrapped .* mag; dims=dim)
    combined ./= sum(mag .* Float32.(TEs); dims=dim)
    dropdims(combined; dims=dim)
end

function scaleandthreshold!(swiphase, mask, mode, level)
    swiphase[.!mask] .= 0
    pos = swiphase .> 0

    if mode == :tanh
        m = median(swiphase[mask .& (swiphase .> 0)])
        m *= 10 / level
        f(x) = (1 + tanh(1 - x/m)) / 2
        swiphase .= f.(swiphase)

    elseif mode == :positive
        swiphase[.!pos] .= 1
        swiphase[pos] .= robustrescale(swiphase[pos], 1, 0.5) .^ level

    elseif mode == :negative
        swiphase[pos] .= 1
        swiphase[.!pos] .= robustrescale(swiphase[.!pos], 0, 1) .^ level

    elseif mode == :triangular
        swiphase[pos] .= robustrescale(swiphase[pos], 1, 0) .^ level
        swiphase[.!pos] .= robustrescale(swiphase[.!pos], 0, 1) .^ level

    else
        error("$mode not defined for SWI")
    end
    swiphase[swiphase .< 0] .= 0
    swiphase
end

function createMIP(S, d=7)
    [minimum(S[x,y,z:z+d-1]) for x in 1:size(S,1), y in 1:size(S,2), z in 1:size(S,3)-d+1]
end
getpixdim(data) = data.header.pixdim[2:(1+ndims(data.mag))]


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 5163d91e-7e01-11eb-3d85-1bde1b52bcc2
using Distributions

# ╔═╡ 4cac0536-7e01-11eb-0f19-29257e61fef3
using StatsPlots

# ╔═╡ 9bcb7778-7e01-11eb-38f5-0d934dfa6201
x = [0:0.01:1;]

# ╔═╡ 5605125a-7e01-11eb-13a6-596920930e6a
begin
	α = @bind a html"<input type='range' min='0.1' max='10' step='0.1' value='5'>"
	β = @bind b html"<input type='range' min='0.1' max='10' step='0.1' value='5'>"
	md"""**distribution parameters:**
	$$\alpha$$: $α
	$$\beta$$: $β
	"""
end

# ╔═╡ 0132d818-7e02-11eb-1a70-bdb0bbcefa7a
md"""
 $$\alpha$$ = $a

 $$\beta$$ = $b
"""

# ╔═╡ 8cb866a0-7e02-11eb-1fdb-eff565b95d5f
begin
	apb = a + b
	m = a / (apb)
	v = a * b / apb * apb * (apb + 1)
	md"""
	Mean: $$\frac{\alpha}{\alpha + \beta}$$ = $m

	Variance: $$\frac{\alpha\beta}{(\alpha + \beta)^2 * (\alpha + \beta +1)}$$ = $v
	"""
end

# ╔═╡ 8fa79d8c-7e01-11eb-19a7-d33fcde6ff4c
B = Beta(a, b)

# ╔═╡ 995a9474-7e01-11eb-348b-cb8a8233b489
plot(x, pdf.(B, x),
	label = "Beta($a, $b)", linewidth = 2)

# ╔═╡ Cell order:
# ╟─9bcb7778-7e01-11eb-38f5-0d934dfa6201
# ╟─5605125a-7e01-11eb-13a6-596920930e6a
# ╟─0132d818-7e02-11eb-1a70-bdb0bbcefa7a
# ╟─8cb866a0-7e02-11eb-1fdb-eff565b95d5f
# ╠═8fa79d8c-7e01-11eb-19a7-d33fcde6ff4c
# ╟─995a9474-7e01-11eb-348b-cb8a8233b489
# ╟─5163d91e-7e01-11eb-3d85-1bde1b52bcc2
# ╟─4cac0536-7e01-11eb-0f19-29257e61fef3

max_stars_repo_name: tim-au/pluto

max_issues_repo_name: tim-au/pluto

max_forks_repo_name: tim-au/pluto


text:
@doc raw"""
    ContextManager(func)

As ContextManager we denote a computation which has a pre-computation and possibly a cleaning up step.

The single argument is supposed to be a function which expects one argument, the continuation function.
Think of it like the following:
```julia
function contextmanagerready(cont)
  # ... do something before
  value = ... # create some value to work on later
  result = cont(value)  # pass the value to the continuation function (think like `yield`)
  # ... do something before exiting, e.g. cleaning up
  result # IMPORTANT: always return the result of the `cont` function
end
```
Now you can wrap it into `ContextManager(contextmanagerready)` and you can use all the context manager
functionalities right away.

There is a simple `@ContextManager` for writing less parentheses
```julia
mycontextmanager(value) = @ContextManager function(cont)
  println("got value = $value")
  result = cont(value)
  println("finished value = $value")
  result
end
```

----------------

You can run it in two ways, either by just passing `Base.identity` as the continuation function
```julia
julia> mycontextmanager(4)(x -> x)
got value = 4
finished value = 4
```
or for convenience we also overload `Base.run`
```julia
# without any extra arguments runs the contextmanager with Base.identity
run(mycontextmanager(4))
# also works with a given continuation, which makes for a nice do-syntax
run(x -> x, mycontextmanager(4))
run(mycontextmanager(4)) do x
  x
end
```
"""
struct ContextManager{F}
  f::F
end

@doc raw"""
    @ContextManager function(cont); ...; end

There is a simple `@ContextManager` for writing less parentheses

```julia
mycontextmanager(value) = @ContextManager function(cont)
  println("got value = $value")
  result = cont(value)
  println("finished value = $value")
  result
end
```
"""
macro ContextManager(func)
  quote
    ContextManager($(esc(func)))
  end
end

# ContextManager is just a wrapper
# pass through function call syntax
(c::ContextManager)(cont) = c.f(cont)
# in principle `cont` could be anything, as anything could be a function, however we have an ambiguity with Base AbstractCmd, and as this is mainly used for do-syntax
# it should be fine to restrict to `Function`
Base.run(cont::Function, c::ContextManager) = c(cont)
Base.run(c::ContextManager) = c(identity)

function Base.eltype(::Type{<:ContextManager{F}}) where F
  Base.promote_op((f, arg) -> f(arg), F, typeof(identity))
end
Base.eltype(::Type{<:ContextManager}) = Any

function Base.foreach(f, c::ContextManager)
  c(f)
  nothing
end

Base.map(f, c::ContextManager) = @ContextManager cont -> begin
  c(x -> cont(f(x)))
end

Iterators.flatten(contextmanager::ContextManager) = @ContextManager cont -> begin
  # execute both nested ContextManagers in the nested manner
  contextmanager() do inner_contextmanager
    convert(ContextManager, inner_contextmanager)(cont)
  end
end