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

@testset ExtendedTestSet "SingleKeypointRegression" begin
    method = BlockMethod(
        (Image{2}(), FastAI.Keypoints{2}(10)),
        (
            ProjectiveTransforms((16, 16), inferencefactor=8),
            ImagePreprocessing(),
            KeypointPreprocessing((16, 16)),
        )
    )
    DLPipelines.checkmethod_core(method)
    FastAI.checkmethod_plot(method)
end


text:
@testset "Function: nrbdegelev" begin
  line = nrbdegelev(nrbline(),2);
  @test line.coefs ≈ [0.0 1/3 2/3 1.0;
                      0.0 0.0 0.0 0.0;
                      0.0 0.0 0.0 0.0;
                      1.0 1.0 1.0 1.0];
  @test line.knots == vec([0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0]);
  @test line.number == 4;
  @test line.order == 4;
  @test typeof(line)<:NURBS1D;

  @test_throws ArgumentError nrbdegelev(line,[2;5]);
  @test_throws ArgumentError nrbdegelev(nrbtestsrf(),5);
end # testset

ext: jl


text:
module OneDimensionalCorrelation

using ArgCheck
using BlockBandedMatrices
using QuadGK
using LinearAlgebra
using PolynomialBases
using ProgressMeter
using Tullio

include("basis.jl")
include("hartree_fock.jl")
include("integrals.jl")
include("ci.jl")

export  Basis,
        BasisLobatto,
        CIHamilton,
        CIOverlap,
        Element1D,
        ElementGrid,
        ElementGridLobatto

export  block_ci,
        bracket,
        ci_matrix,
        ci_vector_product,
        coulomb_matrix,
        derivative_matrix,
        element_ci_orbitals,
        eri,
        erig,
        exchange_matrix,
        fock_matrix,
        fock_matrix!,
        full_ci,
        get_element,
        get_element_indexes,
        get_identity,
        get_length,
        get_weight,
        gram_schmit,
        g_tensor,
        initial_orbitals,
        kinetic_energy,
        metric_tensor,
        one_electron_operator,
        particle_in_box,
        reduced_ci_orbitals,
        solve_hartree_fock

# Write your package code here.

end

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
# To build the documentation:
#    - julia --project="." make.jl
#    - empty!(ARGS); include("make.jl")
# To build the documentation without running the tests:
#    - julia --project="." make.jl preview
#    - push!(ARGS,"preview"); include("make.jl")

# !!! note "An optional title"
#    4 spaces idented
# note, tip, warning, danger, compat

# Format notes:

# # A markdown H1 title
# A non-code markdown normal line

## A comment within the code chunk

#src: line exclusive to the source code and thus filtered out unconditionally

using Pkg
cd(@__DIR__)
Pkg.activate(".")

#Pkg.resolve()
Pkg.instantiate()
#Pkg.add(["Documenter", "Literate", "Glob", "DataFrames", "OdsIO"])

using Test, Documenter, Literate, DocumenterMarkdown, Glob, DataFrames, OdsIO

const LESSONS_ROOTDIR = joinpath(@__DIR__, "lessonsSources")
# Important: If some lesson is removed but the md file is left, this may still be used by Documenter

const LESSONS_ROOTDIR_TMP = joinpath(@__DIR__, "lessonsSources_tmp")
# Where to save the lessons before they are preprocessed

MAKE_PDF = false

LESSONS_SUBDIR = Dict(
  "INTRO - Introduction to the course, Julia and ML"  => "00_-_INTRO_-_Introduction_julia_ml",
  "JULIA1 - Basic Julia programming"              => "01_-_JULIA1_-_Basic_Julia_programming",
  "JULIA2 - Scientific programming with Julia"    => "02_-_JULIA2_-_Scientific_programming_with_Julia",
  "ML1 - Introduction to Machine Learning"        => "03_-_ML1_-_Introduction_to_Machine_Learning",
  "NN - Neural Networks"                          => "04_-_NN_-_Neural_Networks",
  "DT - Decision trees based algorithms [DRAFT]"  => "05_-_DT_-_Decision_trees_based_algorithms"
)

# Utility functions.....

function link_example(content)
    edit_url = match(r"EditURL = \"(.+?)\"", content)[1]
    footer = match(r"^(---\n\n\*This page was generated using)"m, content)[1]
    content = replace(
        content, footer => "[View this file on Github]($(edit_url)).\n\n" * footer
    )
    return content
end

"""
    include_sandbox(filename)
Include the `filename` in a temporary module that acts as a sandbox. (Ensuring
no constants or functions leak into other files.)
"""
function include_sandbox(filename)
    mod = @eval module $(gensym()) end
    return Base.include(mod, filename)
end

function makeList(rootDir,subDirList)
    outArray = []
    for l in sort(collect(subDirList), by=x->x[2])
      #println(l)
      lessonName = l[1]
      lessonDir  = l[2]
      lessonName  = replace(lessonName,"_"=>" ")
      dirArray =[]
      for file in filter(file -> endswith(file, ".md"), sort(readdir(joinpath(rootDir,lessonDir))))
        displayFilename = replace(file,".md"=>"")
        displayFilename = replace(displayFilename,"_"=>" ")
        push!(dirArray,displayFilename=>joinpath(lessonDir,file))
      end
      push!(outArray,lessonName=>dirArray)
    end
    return outArray
end

"""
    rdir(string,match)

Return a vector of all files (full paths) of a given directory recursivelly taht matches `match`, recursivelly.

# example
filenames = getindex.(splitdir.(rdir(LESSONS_ROOTDIR,"*.jl")),2) #get all md filenames
"""
function rdir(dir::AbstractString, pat::Glob.FilenameMatch)
    result = String[]
    for (root, dirs, files) in walkdir(dir)
        append!(result, filter!(f -> occursin(pat, f), joinpath.(root, files)))
    end
    return result
end
rdir(dir::AbstractString, pat::AbstractString) = rdir(dir, Glob.FilenameMatch(pat))

function literate_directory(dir)
    # Removing old compiled md files...
    #for filename in filter(file -> endswith(file, ".md"), readdir(dir))
    #    rm(joinpath(dir,filename))
    #end

    for filename in filter(file -> endswith(file, ".jl"), readdir(dir))
        filenameNoPath = filename
        filename = joinpath(dir,filename)
        # if the md file exist, let's delete it first...
        filenameMD = replace(filename,".jl" => ".md")
        if isfile(filenameMD)
            rm(filenameMD)
        end

        # `include` the file to test it before `#src` lines are removed. It is
        # in a testset to isolate local variables between files.
        if ! ("preview" in ARGS)
            @testset "$(filename)" begin
                println(filename)
                include_sandbox(filename)
             end
             Literate.markdown(
                 filename,
                 dir;
                 documenter = true,
                 postprocess = link_example,
                 # default is @example -> evaluated by documenter at the end of the block
                 codefence =  "```@repl $filenameNoPath" => "```"
             )
        else
            Literate.markdown(
                filename,
                dir;
                documenter = true,
                postprocess = link_example,
                codefence =  "```text" => "```"
            )
        end
    end
    return nothing
end

function preprocess(rootDir)
    cd(@__DIR__)
    Pkg.activate(".")
    #rootDir = LESSONS_ROOTDIR
    files  = rdir(rootDir,"*.md")
    videos = ods_read(joinpath(@__DIR__,"videosList.ods");sheetName="videos",retType="DataFrame")
    for file in files
        #file = files[4]
        origContent = read(file,String)
        outContent = ""
        filename = splitdir(file)[2]
        segmentVideos = videos[videos.host_filename .== filename,:]
        #openVideosFlag = filename == "0001_-_Course_presentation.md" ? "open" : ""

        if size(segmentVideos,1) > 0
            outContent *= """
                        ```@raw html
                        <div id="ytb-videos">
                        <span style=font-weight:bold;>Videos related to this segment (click the title to watch)</span>
                        """
            for (i, video) in enumerate(eachrow(segmentVideos))
                #video = segmentVideos[1,:]
                openVideosFlag = (i == 1) ? "open" : ""
                outContent *= """
                    <details $(openVideosFlag)><summary>$(video.lesson_short_name) - $(video.segment_id)$(video.part_id): $(video.part_name) ($(video.minutes):$(video.seconds))</summary>
                    <div class="container ytb-container">
                        <div class="embed-responsive embed-responsive-16by9">
                            <iframe class="embed-responsive-item" src="https://www.youtube.com/embed/$(video.vid)" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen="" frameborder="0"></iframe>
                        </div>
                    </div>
                    </details>
                    """
            end # end of each video
            outContent *= """
                </div>
                ```
                ------
                """
        end # end of if there are videos
        outContent *= origContent
        if (filename != "index.md")
            commentCode = """
                ```@raw html
                <script src="https://utteranc.es/client.js"
                        repo="sylvaticus/SPMLJ"
                        issue-term="title"
                        label="💬 website_comment"
                        theme="github-dark"
                        crossorigin="anonymous"
                        async>
                </script>
                ```
                """
            addThisCode1 = """
                ```@raw html
                <div class="addthis_inline_share_toolbox"></div>
                ```
                """
            addThisCode2 = """
                ```@raw html
                <!-- Go to www.addthis.com/dashboard to customize your tools -->
                <script type="text/javascript" src="//s7.addthis.com/js/300/addthis_widget.js#pubid=ra-6256c971c4f745bc"></script>
                ```
                """
            # https://crowdsignal.com/support/rating-widget/
            ratingCode1 = """
                ```@raw html
                <div id="pd_rating_holder_8962705"></div>
                <script type="text/javascript">
                const pageURL = window.location.href;
                PDRTJS_settings_8962705 = {
                "id" : "8962705",
                "unique_id" : "$(file)",
                "title" : "$(filename)",
                "permalink" : pageURL
                };
                </script>
                ```
                """
            ratingCode2 = """
                ```@raw html
                <script type="text/javascript" charset="utf-8" src="https://polldaddy.com/js/rating/rating.js"></script>
                ```
                """
            outContent *= "\n---------\n"
            outContent *= ratingCode1
            outContent *= addThisCode1
            outContent *= "\n---------\n"
            outContent *= commentCode
            outContent *= ratingCode2
            outContent *= addThisCode2
        end
        write(file,outContent)
    end # end for each file
end #end preprocess function

# ------------------------------------------------------------------------------
# Saving the unmodified source to a temp directory
cp(LESSONS_ROOTDIR, LESSONS_ROOTDIR_TMP; force=true)

println("Starting literating tutorials (.jl --> .md)...")
literate_directory.(map(lsubdir->joinpath(LESSONS_ROOTDIR ,lsubdir),values(LESSONS_SUBDIR)))

if MAKE_PDF
    println("Starting making PDF...")
    makedocs(sitename="SPMLJ - Introduction to Scientific Programming and Machine Learning with Julia",
            authors = "Antonello Lobianco",
            pages = [
                "Index" => "index.md",
                "Lessons" => makeList(LESSONS_ROOTDIR,LESSONS_SUBDIR),
            ],
            format = Markdown(),
            source  = "lessonsSources", # Attention here !!!!!!!!!!!
            build   = "buildedDoc_PDF",
    )
end

println("Starting preprocessing markdown pages...")
preprocess(LESSONS_ROOTDIR)

println("Starting making the documentation...")
makedocs(sitename="SPMLJ",
         authors = "Antonello Lobianco",
         pages = [
            "Index" => "index.md",
            "Lessons" => makeList(LESSONS_ROOTDIR,LESSONS_SUBDIR),
         ],
         format = Documenter.HTML(
             prettyurls = false,
             analytics = "G-Q39LHCRBB6",
             assets = ["assets/custom.css"],
             ),
         #strict = true,
         #doctest = false
         source  = "lessonsSources", # Attention here !!!!!!!!!!!
         build   = "buildedDoc",
         #preprocess = preprocess
)

# Copying back the unmodified source
cp(LESSONS_ROOTDIR_TMP, LESSONS_ROOTDIR; force=true)

println("Starting deploying the documentation...")
deploydocs(
    repo = "github.com/sylvaticus/SPMLJ.git",
    devbranch = "main",
    target = "buildedDoc"
)

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
const GenericNastranDeck = Vector{GenericCard}
const NastranDeck = Vector{NastranCard}

function next_token(line::IOBuffer,N::Int64)
    token = String(read(line,N))
    n = length(token)
    i = searchindex(token,'\t')
    if i > 0
        skip(line,i-n)
        strip(token[1:i-1])
    else
        strip(token)
    end
end

function next_token_comma(line::IOBuffer)
    token = readuntil(line,',')
    if length(token) > 0 && token[end] == ','
        strip(token[1:end-1])
    else
        strip(token)
    end
end

function parse_line{T<:AbstractString}(line::T)
    buf = IOBuffer(line)
    card = T[]
    if ',' in line
        tok = next_token_comma(buf)
        if length(tok) > 0 && tok[end] == '*'
            n = 5
            tok = tok[1:end-1]
        else
            n = 9
        end
        push!(card,tok)
        for i in 1:n
            push!(card,next_token_comma(buf))
        end
    else
        tok = next_token(buf,8)
        if length(tok) > 0 && tok[end] == '*'
            l = 16
            n = 4
            tok = tok[1:end-1]
        else
            l = 8
            n = 8
        end
        push!(card,tok)
        for _ in 1:n
            push!(card,next_token(buf,l))
        end
        push!(card,next_token(buf,8))
    end
    card
end

function remove_comments(line::AbstractString)
    i = searchindex(line,'$')
    if i == 0
        line
    else
        line[1:i-1]
    end
end

isnastranfloat = r"^([+-]?[0-9]*(?:[0-9]\.|\.[0-9])[0-9]*)[dD]?([+-][0-9]+)$"
isfloat = r"^([+-]?[0-9]*(?:[0-9]\.|\.[0-9])[0-9]*)(?:[eE]([+-][0-9]+))?$"
isint = r"^(?:[+-][0-9]{1,7}|[0-9]{1,8})$"
function parse_number(field::String)
    m = match(isnastranfloat,field)
    if m != nothing
        parse(Float64,join(m.captures,"e"))
    elseif ismatch(isfloat,field)
        parse(Float64,field)
    elseif ismatch(isint,field)
        parse(Int64,field)
    else
        field
    end
end

function parse_numbers(card::Vector{String})
    map(parse_number,card)
end

function read_cards(filename::AbstractString)
    map(eachline(open(filename))) do line
        line |> chomp |> remove_comments |> parse_line |> parse_numbers
    end
end

immutable NastranCardIterator
    cards
end

Base.start(I::NastranCardIterator) = (Dict{String,GenericCard}(),start(I.cards))

function Base.next(I::NastranCardIterator,state)
    cont_cards, cards_state = state
    while !done(I.cards,cards_state)
        card, cards_state = next(I.cards,cards_state)
        if all(f -> f == "",card)
            if haskey(cont_cards,"+") && haskey(cont_cards,"")
                c = pop!(cont_cards,"")
                cont_cards[""] = pop!(cont_cards,"+")
                return (c,(cont_cards,cards_state))
            elseif haskey(cont_cards,"+")
                cont_cards[""] = pop!(cont_cards,"+")
            end
        elseif card[1] == "+" || card[1] == ""
            if !haskey(cont_cards,"+")
                @show card
                @show cont_cards
                error()
            end
            append!(cont_cards["+"],card[2:end-1])
        else
            if haskey(cont_cards,"+") && haskey(cont_cards,"")
                c = pop!(cont_cards,"")
                cont_cards[""] = pop!(cont_cards,"+")
                cont_cards["+"] = card[1:end-1]
                return (c,(cont_cards,cards_state))
            elseif haskey(cont_cards,"+")
                cont_cards[""] = pop!(cont_cards,"+")
            end
            cont_cards["+"] = card[1:end-1]
        end
    end
    if haskey(cont_cards,"")
        c = pop!(cont_cards,"")
        return (c,(cont_cards,cards_state))
    end
    if haskey(cont_cards,"+")
        c = pop!(cont_cards,"+")
        return (c,(cont_cards,cards_state))
    end
end
Base.done(I::NastranCardIterator,state) = done(I.cards,state[2]) && length(state[1]) == 0

function GenericNastranDeck(filename::AbstractString)
    cards = read_cards(filename)
    GenericNastranDeck(collect(NastranCardIterator(cards)))
end

function NastranDeck(filename::AbstractString)
    cards = read_cards(filename)
    deck = []
    for card in NastranCardIterator(cards)
        push!(deck,card)
    end
    deck
end


text:
function getPolynomial(X, Y)
    n = length(Y)
    p = string("P = x -> ", Y[1])
    for i::Int in 2:n
        P = string(p, " + ", Y[i])
        for j in 1:i-1
            p = string(p, "*(x-", X[j], ")")
        end
    end
    return eval(Meta.parse(p))
end

function getCoefficients(X, Y)

    n = length(X)

    for i::Int in 2:n
        for j::Int in 2:i
            Y[i, j] = (Y[i, j-1] - Y[i-1, j-1])/(X[i] - X[i-j+1])
        end
    end

    F = []

    for i in 1:n
        push!(F, Y[i, i])
    end

    return F
end

"""
    divid(X, Y, cff)

Computes interpolating polynomial passing through the points whoose coordinates are provided in `X[]` and `Y[]`
using Newton's divided difference method. Returns `::Array{Float64}` if `cff = true` (default), or the interpolating polynomial P if `cff = false`.

    divid(X::Array{Float64}, Y::Array{Float64}, cff::Bool=true)

### OUTPUT:
  * `F` - coefficiants of the interpolating polynomial of the type `Array{Float64}`,
  * `P` - interpolating polynomial.
"""
function divid(X::Array{Float64}, Y::Array{Float64}, cff::Bool=true)
    if cff
        return getCoefficients(X, Y)
    else
        return getPolynomial(X, Y)
    end
end


text:
# Constraints

"""
    supports_constraint(model::ModelLike, ::Type{F}, ::Type{S})::Bool where {F<:AbstractFunction,S<:AbstractSet}

Return a `Bool` indicating whether `model` supports `F`-in-`S` constraints, that is,
`copy_to(model, src)` does not return `CopyUnsupportedConstraint` when `src` contains `F`-in-`S` constraints.
If `F`-in-`S` constraints are only not supported in specific circumstances, e.g. `F`-in-`S` constraints cannot be combined with another type of constraint, it should still return `true`.
"""
supports_constraint(model::ModelLike, ::Type{<:AbstractFunction}, ::Type{<:AbstractSet}) = false

"""
    struct UnsupportedConstraint{F<:AbstractFunction, S<:AbstractSet} <: UnsupportedError
        message::String # Human-friendly explanation why the attribute cannot be set
    end

An error indicating that constraints of type `F`-in-`S` are not supported by
the model, i.e. that [`supports_constraint`](@ref) returns `false`.
"""
struct UnsupportedConstraint{F<:AbstractFunction, S<:AbstractSet} <: UnsupportedError
    message::String # Human-friendly explanation why the attribute cannot be set
end
UnsupportedConstraint{F, S}() where {F, S} = UnsupportedConstraint{F, S}("")

element_name(::UnsupportedConstraint{F, S}) where {F, S} = "`$F`-in-`$S` constraint"

"""
    struct AddConstraintNotAllowed{F<:AbstractFunction, S<:AbstractSet} <: NotAllowedError
        message::String # Human-friendly explanation why the attribute cannot be set
    end

An error indicating that constraints of type `F`-in-`S` are supported (see
[`supports_constraint`](@ref)) but cannot be added.
"""
struct AddConstraintNotAllowed{F<:AbstractFunction, S<:AbstractSet} <: NotAllowedError
    message::String # Human-friendly explanation why the attribute cannot be set
end
AddConstraintNotAllowed{F, S}() where {F, S} = AddConstraintNotAllowed{F, S}("")

operation_name(::AddConstraintNotAllowed{F, S}) where {F, S} = "Adding `$F`-in-`$S` constraints"

"""
    struct ScalarFunctionConstantNotZero{T, F, S} <: Exception
        constant::T
    end

An error indicating that the constant part of the function in the constraint
`F`-in-`S` is nonzero.
"""
struct ScalarFunctionConstantNotZero{T, F, S} <: Exception
    constant::T
end

function Base.showerror(io::IO,
                    err::ScalarFunctionConstantNotZero{T, F, S}) where {T, F, S}
    print(io, "In `$F`-in-`$S` constraint: Constant $(err.constant) of the ",
          "function is not zero. The function constant should be moved to the ",
          "set. You can use `MOI.Utilities.add_scalar_constraint` which does ",
          "this automatically.")
end

"""
    add_constraint(model::ModelLike, func::F, set::S)::ConstraintIndex{F,S} where {F,S}

Add the constraint ``f(x) \\in \\mathcal{S}`` where ``f`` is defined by `func`, and ``\\mathcal{S}`` is defined by `set`.

    add_constraint(model::ModelLike, v::VariableIndex, set::S)::ConstraintIndex{SingleVariable,S} where {S}
    add_constraint(model::ModelLike, vec::Vector{VariableIndex}, set::S)::ConstraintIndex{VectorOfVariables,S} where {S}

Add the constraint ``v \\in \\mathcal{S}`` where ``v`` is the variable (or vector of variables) referenced by `v` and ``\\mathcal{S}`` is defined by `set`.

* An [`UnsupportedConstraint`](@ref) error is thrown if `model` does not support
  `F`-in-`S` constraints,
* a [`AddConstraintNotAllowed`](@ref) error is thrown if it supports `F`-in-`S`
  constraints but it cannot add the constraint(s) in its current state and
* a [`ScalarFunctionConstantNotZero`](@ref) error may be thrown if
  `func` is an `AbstractScalarFunction` with nonzero constant and `set`
  is [`EqualTo`](@ref), [`GreaterThan`](@ref), [`LessThan`](@ref) or
  [`Interval`](@ref).
* a [`LowerBoundAlreadySet`](@ref) error is thrown if `F` is a
  [`SingleVariable`](@ref) and a constraint was already added to this variable
  that sets a lower bound.
* a [`UpperBoundAlreadySet`](@ref) error is thrown if `F` is a
  [`SingleVariable`](@ref) and a constraint was already added to this variable
  that sets an upper bound.
"""
function add_constraint(model::ModelLike, func::AbstractFunction,
                        set::AbstractSet)
    throw_add_constraint_error_fallback(model, func, set)
end

# throw_add_constraint_error_fallback checks whether func and set are both
# scalar or both vector. If it is the case, it calls
# `correct_throw_add_constraint_error_fallback`
function throw_add_constraint_error_fallback(model::ModelLike,
                                             func::AbstractScalarFunction,
                                             set::AbstractScalarSet;
                                             kwargs...)
    correct_throw_add_constraint_error_fallback(model, func, set; kwargs...)
end
function throw_add_constraint_error_fallback(model::ModelLike,
                                             func::AbstractVectorFunction,
                                             set::AbstractVectorSet;
                                             kwargs...)
    correct_throw_add_constraint_error_fallback(model, func, set; kwargs...)
end
function throw_add_constraint_error_fallback(model::ModelLike,
                                             func::AbstractScalarFunction,
                                             set::AbstractVectorSet;
                                             kwargs...)
    error("Cannot add a constraint of the form `ScalarFunction`-in-`VectorSet`")
end
function throw_add_constraint_error_fallback(model::ModelLike,
                                             func::AbstractVectorFunction,
                                             set::AbstractScalarSet;
                                             kwargs...)
    error("Cannot add a constraint of the form `VectorFunction`-in-`ScalarSet`")
end

# func and set are both scalar or both vector
function correct_throw_add_constraint_error_fallback(model::ModelLike,
                                                     func::AbstractFunction,
                                                     set::AbstractSet;
                                                     error_if_supported=AddConstraintNotAllowed{typeof(func), typeof(set)}())
    if supports_constraint(model, typeof(func), typeof(set))
        throw(error_if_supported)
    else
        throw(UnsupportedConstraint{typeof(func), typeof(set)}())
    end
end

# convenient shorthands TODO: document
add_constraint(model::ModelLike, v::VariableIndex, set::AbstractScalarSet) = add_constraint(model, SingleVariable(v), set)
add_constraint(model::ModelLike, v::Vector{VariableIndex}, set::AbstractVectorSet) = add_constraint(model, VectorOfVariables(v), set)

"""
    add_constraints(model::ModelLike, funcs::Vector{F}, sets::Vector{S})::Vector{ConstraintIndex{F,S}} where {F,S}

Add the set of constraints specified by each function-set pair in `funcs` and `sets`. `F` and `S` should be concrete types.
This call is equivalent to `add_constraint.(model, funcs, sets)` but may be more efficient.
"""
function add_constraints end

# default fallback
add_constraints(model::ModelLike, funcs, sets) = add_constraint.(model, funcs, sets)

"""
    LowerBoundAlreadySet{S1, S2}

Error thrown when setting a `SingleVariable`-in-`S2` when a
`SingleVariable`-in-`S1` has already been added and the sets `S1`, `S2` both
set a lower bound, i.e. they are [`EqualTo`](@ref), [`GreaterThan`](@ref),
[`Interval`](@ref), [`Semicontinuous`](@ref) or [`Semiinteger`](@ref).
"""
struct LowerBoundAlreadySet{S1, S2}
    vi::VariableIndex
end

function Base.showerror(io::IO, err::LowerBoundAlreadySet{S1, S2}) where {S1, S2}
    print(io, typeof(err), ": Cannot add `SingleVariable`-in`", S2,
          "` constraint for variable ", err.vi, " as a `SingleVariable`-in`",
          S1, "` constraint was already set for this variable and both",
          " constraints set a lower bound.")
end

"""
    UpperBoundAlreadySet{S1, S2}

Error thrown when setting a `SingleVariable`-in-`S2` when a
`SingleVariable`-in-`S1` has already been added and the sets `S1`, `S2` both
set an upper bound, i.e. they are [`EqualTo`](@ref), [`LessThan`](@ref),
[`Interval`](@ref), [`Semicontinuous`](@ref) or [`Semiinteger`](@ref).
"""
struct UpperBoundAlreadySet{S1, S2}
    vi::VariableIndex
end

function Base.showerror(io::IO, err::UpperBoundAlreadySet{S1, S2}) where {S1, S2}
    print(io, typeof(err), ": Cannot add `SingleVariable`-in`", S2,
          "` constraint for variable ", err.vi, " as a `SingleVariable`-in`",
          S1, "` constraint was already set for this variable and both",
          " constraints set an upper bound.")
end

"""
## Transform Constraint Set

    transform(model::ModelLike, c::ConstraintIndex{F,S1}, newset::S2)::ConstraintIndex{F,S2}

Replace the set in constraint `c` with `newset`. The constraint index `c`
will no longer be valid, and the function returns a new constraint index with
the correct type.

Solvers may only support a subset of constraint transforms that they perform
efficiently (for example, changing from a `LessThan` to `GreaterThan` set). In
addition, set modification (where `S1 = S2`) should be performed via the
`modify` function.

Typically, the user should delete the constraint and add a new one.

### Examples

If `c` is a `ConstraintIndex{ScalarAffineFunction{Float64},LessThan{Float64}}`,

```julia
c2 = transform(model, c, GreaterThan(0.0))
transform(model, c, LessThan(0.0)) # errors
```
"""
function transform end

# default fallback
function transform(model::ModelLike, c::ConstraintIndex, newset)
    f = get(model, ConstraintFunction(), c)
    delete(model, c)
    add_constraint(model, f, newset)
end


text:
DeepCFR.vectorized(::MatrixGame, I) = SA[Float32(I)]

function DeepCFR.vectorized(game::Kuhn, I)
    p, pc, hist = I
    h = convert(SVector{3,Float32}, hist)
    SA[Float32(p), Float32(pc), h...]
end

@testset "Solver" begin

    ## Basic operation + reasonable outputs
    game = Kuhn()
    sol = DeepCFRSolver(game; buffer_size=50, traversals=10)
    train!(sol, 100)
    I0 = infokey(game, initialhist(game))
    σ = @inferred sol(DeepCFR.vectorized(game, I0))

    @test σ isa Vector{<:AbstractFloat}
    @test sum(σ) ≈ 1.0
    @test DeepCFR.on_gpu(sol) == false
    @test DeepCFR.infotype(sol) == SVector{5,Float32}
    @test DeepCFR.regrettype(sol) == Vector{Float32}
    @test DeepCFR.strattype(sol) == Vector{Float32}

    game = MatrixGame()
    sol = DeepCFRSolver(game)
    train!(sol, 10)
    I0 = infokey(game, initialhist(game))
    σ = @inferred sol(DeepCFR.vectorized(game, I0))
    @test σ isa Vector{<:AbstractFloat}
    @test sum(σ) ≈ 1.0
    @test DeepCFR.on_gpu(sol) == false
    @test DeepCFR.infotype(sol) == SVector{1,Float32}
    @test DeepCFR.regrettype(sol) == Vector{Float32}
    @test DeepCFR.strattype(sol) == Vector{Float32}

    sol = DeepCFRSolver(game; traverser=DeepCFR.OutcomeSample())
    train!(sol, 10)
    I0 = infokey(game, initialhist(game))
    σ = @inferred sol(DeepCFR.vectorized(game, I0))
    @test σ isa Vector{<:AbstractFloat}
    @test sum(σ) ≈ 1.0
    @test DeepCFR.on_gpu(sol) == false
    @test DeepCFR.infotype(sol) == SVector{1,Float32}
    @test DeepCFR.regrettype(sol) == Vector{Float32}
    @test DeepCFR.strattype(sol) == Vector{Float32}

    game = MatrixGame()
    sol = DeepCFRSolver(game, on_gpu=true)
    train!(sol, 10)
    I0 = infokey(game, initialhist(game))
    σ = @inferred sol(DeepCFR.vectorized(game, I0))
    @test σ isa Vector{<:AbstractFloat}
    @test sum(σ) ≈ 1.0
    @test DeepCFR.on_gpu(sol) == true

    ## Printing

    io = IOBuffer()
    m = MIME"text/plain"()
    Base.show(io, m, sol)
    str = String(take!(io))
    @test occursin("GPU", str)
    @test occursin("Advantage Net", str)
    @test occursin("Advantage Memory", str)
    @test occursin("Strategy Net", str)
    @test occursin("Strategy Memory", str)
    @test occursin("Batch Size", str)
    @test occursin("Traversals", str)
    @test occursin("Advantage Optimizer", str)
    @test occursin("Strategy Optimizer", str)
    @test occursin("Game", str)
end

@testset "Regret Matching" begin
    v1 = Float32.([1,2,3])
    σ1 = DeepCFR.regret_match_strategy(v1)
    @test all(σ1 .≈ Float32.([1/6,2/6,3/6]))

    v2 = Float32.([-1,2,3])
    σ2 = DeepCFR.regret_match_strategy(v2)
    @test all(σ2 .≈ Float32.([0,2/5,3/5]))

    v3 = Float32.([1,-2,3])
    σ3 = DeepCFR.regret_match_strategy(v3)
    @test all(σ3 .≈ Float32.([1/4,0,3/4]))

    v4 = Float32.([1,2,-3])
    σ4 = DeepCFR.regret_match_strategy(v4)
    @test all(σ4 .≈ Float32.([1/3,2/3,0]))

    v5 = Float32.([-1,-2,-3])
    σ5 = DeepCFR.regret_match_strategy(v5)
    @test all(σ5 .≈ Float32.([1,0,0]))

    v6 = Float32.([-3,-2,-1])
    σ6 = DeepCFR.regret_match_strategy(v6)
    @test all(σ6 .≈ Float32.([0,0,1]))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 819e0a72-cc97-4910-8f05-91ebbb2273d8
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ bb2a66f4-1c5e-4a4a-a0e6-965e0996eae7
md"""
# Figure 4.18 Continuation diagram

**NOTE**

[Bifurcations.jl](https://github.com/tkf/Bifurcations.jl) does not work on Julia v1.6 (As I am writing this).
And [BifurcationKit.jl](https://github.com/rveltz/BifurcationKit.jl) might be too complex for this example.

"""

# ╔═╡ 0549e1c5-16b5-4b3e-aca2-75d41ee93451
p = (K1=20.0, K2=5.0, K3=5.0, K4=5.0, K5=2.0, N=4)

# ╔═╡ a80a7285-da0b-463e-8d69-165d2b797c3b
k1Range = LinRange(0.0, 1000.0, 100)

# ╔═╡ a32c3e1e-eac7-47fa-b9cb-2ac919cd4c56
u0 = LVector(A=0.0, B=0.0)

# ╔═╡ c4725c60-bded-11eb-2b5a-c57f48c64c55
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ a05c0b69-f761-4920-ac17-ebf437f2862e
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ 2c15803d-104a-4b87-8556-0238ce4042a2
function model!(du, u, p, t)
    @unpack K1, K2, K3, K4, K5, N = p
    @unpack A, B = u

    v1 = K1 * hill(1, B, N)
    v5 = K5 * A

    du.A = v1 - v5 - K3 * A
    du.B = K2 + v5 - K4 * B
end

# ╔═╡ 4577a459-7151-4e53-a235-1f35ffa49eff
# See also :
# Ensemble analysis: https://diffeq.sciml.ai/stable/features/ensemble/
a = map(k1Range) do k1
	p1 = @set p.K1 = k1
	prob = SteadyStateProblem(model!, u0, p1)
	sol = solve(prob)
	sol[1]
end

# ╔═╡ a9766d07-4ed0-40b5-a00d-a1694413e289
plot(k1Range, a, title = "Fig 4.18 Continuation diagram", xlabel = "K1" , ylabel= "Steady state [A]", leg=nothing, ylim=(0.0, 4.0))

# ╔═╡ Cell order:
# ╠═bb2a66f4-1c5e-4a4a-a0e6-965e0996eae7
# ╠═2c15803d-104a-4b87-8556-0238ce4042a2
# ╠═0549e1c5-16b5-4b3e-aca2-75d41ee93451
# ╠═a80a7285-da0b-463e-8d69-165d2b797c3b
# ╠═a32c3e1e-eac7-47fa-b9cb-2ac919cd4c56
# ╠═4577a459-7151-4e53-a235-1f35ffa49eff
# ╠═a9766d07-4ed0-40b5-a00d-a1694413e289
# ╠═c4725c60-bded-11eb-2b5a-c57f48c64c55
# ╠═a05c0b69-f761-4920-ac17-ebf437f2862e
# ╠═819e0a72-cc97-4910-8f05-91ebbb2273d8

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
export
    Conditional,
    extend_tuple_type

using Folds
using StatsFuns

@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)
tuplejoin(t::Tuple) = t

"""
    extend_tuple_type(T1, T2)

Given two types `T1` and `T2`, concatenate the types into a single
tuple type.

# Arguments
- `T1`: Any type
- `T2`: A tuple type

# Returns
- If `T1` is a tuple type, a tuple with the concatenation of the types in `T1` and `T2`
- If `T1` is not a tuple type, a tuple with `T1` prepended to the types in `T2`

# Examples
```
julia> extend_tuple_type(Int64, Tuple{Float64})
Tuple{Int64, Float64}

julia> extend_tuple_type(Tuple{Int64}, Tuple{Float64})
Tuple{Int64, Float64}

julia> extend_tuple_type(Tuple{Vector{Float64}}, Tuple{Symbol,Symbol})
Tuple{Vector{Float64}, Symbol, Symbol}
```
"""
function extend_tuple_type(T1, T2)
    if T1 <: Tuple
        Tuple{tuplejoin(fieldtypes(T1), fieldtypes(T2))...}
    else
        Tuple{T1, fieldtypes(T2)...}
    end
end

"""
    abstract type Conditional{I <: Tuple, J <: Tuple, K <: Tuple, O, S <: SFunc{J, O}} <: SFunc{K, O}

`Conditional` *sfuncs* represent the generation of an sfunc depending on the values of parents.  An
subtype of `Conditional` must provide a `gensf` method that takes an `I` and returns an
`SFunc{J,O}` (**important** the generated SFunc must not appear outside the Conditional.
It should not be a parent).

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `logcpdf`
- `make_factors`
- `compute_pi`
- `send_lambda`

# Type parameters
- `I`: the type of data used to generate the `Conditional`'s *sfunc*
- `J`: a tuple that represents the input types (the `I`) of the `Conditional`'s generated *sfunc*
- `K`: the input types of the `Conditional`; this is a tuple of types constructed from `I`,
and `J` using `extend_tuple_types`
- `O`: the output type(s) of both the `Conditional` and the `Conditional`'s generated *sfunc*
- `S`: the type of the `Conditional`'s generated *sfunc*

"""
abstract type Conditional{I, J <: Tuple, K <: Tuple, O, S <: SFunc{J, O}} <: SFunc{K, O}
end

function split_pars(::Conditional{I}, pars) where {I}
    if I <: Tuple
        n1 = length(fieldnames(I))
        (Tuple(pars[1:n1]), Tuple(pars[n1+1:length(pars)]))
    else
        (pars[1], Tuple(pars[2:length(pars)]))
    end
end

@impl begin
    struct ConditionalSupport end

    function support(sf::Conditional{I,J,K,O},
                    parranges::NTuple{N,Vector},
                    size::Integer,
                    curr::Vector{<:O}) where {I,J,K,O,N}
        (iranges, jranges) = split_pars(sf, parranges)
        irgs::Array{Array{Any, 1}, 1} = [r for r in iranges]
        isrange = cartesian_product(irgs)
        rng::Vector = isempty(curr) ? Vector{O}() : copy(curr)
        # inc = Int(ceil(size / length(isrange)))
        allcombranges= Dict{Any, Vector{O}}()

        # create complete range regardless of the size
        for is in isrange
            gsf = gensf(sf,tuple(is...))
            newrng = copy(support(gsf, jranges, size, curr))
            allcombranges[is] = unique(newrng)
        end

        # create range by mixing values from ranges created by each combination of parents values up to size
        allrng = collect(Iterators.flatten(values(allcombranges)))
        max_size = length(unique(allrng))
        while length(rng) < min(size, max_size)
            for is in isrange
                if(length(allcombranges[is]) > 0)
                    val = popfirst!(allcombranges[is])
                    push!(rng, val)
                end
            end
            rng = unique(rng)
        end
        return isempty(rng) ? O[] : rng
    end
end

@impl begin
    struct ConditionalSample end

    function sample(sf::Conditional{I,J,K,O}, i::K)::O where {I,J,K,O}
        (ivals, jvals) = split_pars(sf, i)
        sfg = gensf(sf,ivals)
        return sample(sfg, jvals)
    end
end

@impl begin
    struct ConditionalLogcpdf end

    function logcpdf(sf::Conditional{I,J,K,O}, i::K, o::O)::AbstractFloat where {I,J,K,O}
        (ivals, jvals) = split_pars(sf, i)
        sfgen = gensf(sf,ivals)
        return logcpdf(sfgen, jvals, o)
    end
end

@impl begin
    struct ConditionalSupportQuality end

    function support_quality(sf::Conditional{I,J,K,O,S}, parranges) where {I,J,K,O,S}
        q = support_quality_rank(:CompleteSupport)
        (iranges, jranges) = split_pars(sf, parranges)
        isrange = cartesian_product(iranges)
        for is in isrange
            gsf = gensf(sf,tuple(is...))
            imp = get_imp(MultiInterface.get_policy(), Support, gsf, jranges, 0, O[])
            q = min(q, support_quality_rank(support_quality(imp, gsf, jranges)))
        end
        return support_quality_from_rank(q)
    end
end

@impl begin
    mutable struct ConditionalMakeFactors
        numpartitions::Int64 = 1
    end

    function make_factors(sf::Conditional{I,J,K,O},
                        range::__OptVec{<:O},
                        parranges::NTuple{N,Vector},
                        id,
                        parids::Tuple)::Tuple{Vector{<:Scruff.Utils.Factor}, Vector{<:Scruff.Utils.Factor}} where {I,J,K,O,N}

        if any(isempty, parranges) # This is possible in lazy inference
            # Just return an empty factor so as not to break
            keys = (id,)
            dims = (0,)
            entries = Float64[]
            fs = [Factor(dims, keys, entries)]
            return (fs, fs)
        end

        # Here is a very general method for constructing factors for conditional sfuncs.
        # We introduce a switch variable.
        # For each of the generated sfuncs, we make its factors and extend them by saying
        # they are only relevant when the switch value corresponds to that choice.
        # We also add a factor saying that the switch variable takes on the appropriate value
        # for each value of the i inputs.
        # For this switch factor, we use another decomposition, relating the switch variable
        # to each of the i inputs separately. This avoids a quadratic blowup and brings the
        # complexity analogous to that of less general methods.
        # This can be improved by making a special case where there is only one i parent,
        # and returning a single factor.
        (iranges, jranges) = split_pars(sf, parranges)
        (iids, jids) = split_pars(sf, parids)
        icombos = cartesian_product(iranges)
        lfs = Factor[]
        ufs = Factor[]
        switchkey = nextkey()
        switchsize = length(icombos)
        for (switchval, ivals) in enumerate(icombos)
            function extend(factor)
                keys = [k for k in factor.keys]
                push!(keys, switchkey)
                keys = Tuple(keys)
                dims = [d for d in factor.dims]
                push!(dims, switchsize)
                dims = Tuple(dims)
                entries = Float64[]
                for e in factor.entries
                    for k in 1:switchsize
                        # 1.0 is the irrelevant value when factors are multiplied
                        push!(entries, k == switchval ? e : 1.0)
                    end
                end
                result = Factor(dims, keys, entries)
                return result
            end
            subsf = gensf(sf,tuple(ivals...))
            (sublfs, subufs) = make_factors(subsf, range, jranges, id, jids)
            append!(lfs, [extend(f) for f in sublfs])
            append!(ufs, [extend(f) for f in subufs])
        end

        for i in 1:length(iranges)
            switchfactorkeys = (parids[i], switchkey)
            parsize = length(parranges[i])
            switchfactordims = (parsize, switchsize)
            switchfactorentries = Float64[]
            for j in 1:parsize
                for ivals in icombos
                    # The effect of this is that the product will only be 1 for all the corresponding
                    # i values.
                    push!(switchfactorentries, j == indexin([ivals[i]], parranges[i])[1] ? 1.0 : 0.0)
                end
            end
            switchfact = Factor(switchfactordims, switchfactorkeys, switchfactorentries)
            push!(lfs, switchfact)
            push!(ufs, switchfact)
        end

        return (lfs, ufs)
    end
end

@impl begin
    struct ConditionalComputePi end

    function compute_pi(sf::Conditional{I,J,K,O},
                     range::__OptVec{<:O},
                     parranges::NTuple{N,Vector},
                     incoming_pis::Tuple)::Dist{<:O} where {N,I,J,K,O}

        (iranges, jranges) = split_pars(sf, parranges)
        # We need to correctly handle the case with duplicate parent values
        # logcpdf below will get the full incoming pi for all duplicate values
        irgs::Vector{Set{Any}} = [Set(r) for r in iranges]
        icombos = cartesian_product(irgs)
        (ipis, jpis) = split_pars(sf, incoming_pis)
        result_pieces = [Vector{Float64}() for x in range]

        for is in icombos
            gsf = gensf(sf,tuple(is...))
            ps = compute_pi(gsf, range, Tuple(jranges), jpis)
            if (!isempty(iranges))
                # ipi = sum([logcpdf(ipis[i], (), irgs[i][ind[i]]) for i in 1:length(iranges)])
                ipi = sum([logcpdf(ipis[j], (), is[j]) for j in 1:length(iranges)])
            else
                ipi = 1
            end
            for i in 1:length(range)
                push!(result_pieces[i], ipi + logcpdf(ps, (), range[i]))
            end
        end
        result = normalize([exp(StatsFuns.logsumexp(rps)) for rps in result_pieces])
        return Cat(range, result)
    end
end

@impl begin
    struct ConditionalSendLambda end

    function send_lambda(sf::Conditional{I,J,K,O},
                       lambda::Score{<:O},
                       range::__OptVec{<:O},
                       parranges::NTuple{N,Vector},
                       incoming_pis::Tuple,
                       parent_ix::Integer)::Score where {N,I,J,K,O}
        (iranges, jranges) = split_pars(sf, parranges)
        (ipis, jpis) = split_pars(sf, incoming_pis)
        iar :: Array{Array{Int, 1}} = collect(map(r -> collect(1:length(r)), iranges))
        iinds = cartesian_product(iar)
        jar :: Array{Array{Int, 1}} = collect(map(r -> collect(1:length(r)), jranges))
        jinds = cartesian_product(jar)
        result = zeros(Float64, length(parranges[parent_ix]))
        # Need to make sure the target parent range is a Vector{T} rather than a Vector{Any}
        T = typeof(parranges[parent_ix][1])
        target_parrange :: Vector{T} = parranges[parent_ix]

        if parent_ix <= length(iranges)
            rs::Array{Array{Int, 1}, 1} = [collect(1:length(rg)) for rg in iranges]
            deleteat!(rs, parent_ix)
            restranges = cartesian_product(rs)
            # We want to send a lambda message to the specific I input.
            # We can ignore the pi message from this input.
            # For each value of this input, we want to compute \sum_{i',j,o} \pi_J(j) P_ps.sf(o | j ; gen_params(ps, (i, i'))) \lambda(o).
            for i = 1:length(iranges[parent_ix])
                ival = iranges[parent_ix][i]
                ls = Vector{Float64}()
                for rest in restranges
                    ipi = 0
                    fullvals = []
                    for r = 1:parent_ix-1
                        ipi += logcpdf(ipis[r], (), parranges[r][rest[r]])
                        push!(fullvals, iranges[r][rest[r]])
                    end
                    push!(fullvals, ival)
                    for r = parent_ix+1:length(iranges)
                        ipi += logcpdf(ipis[r], (), parranges[r][rest[r-1]])
                        push!(fullvals, iranges[r][rest[r-1]])
                    end
                    gsf = gensf(sf,tuple(fullvals...))
                    for jind in jinds
                        jpi = 0
                        jval = []
                        for r = 1:length(jpis)
                            jpi += logcpdf(jpis[r], (), jind[r])
                            push!(jval, jranges[r][jind[r]])
                        end
                        jvaltypes = [typeof(j) for j in jval]
                        tjval = convert(Vector{typejoin(jvaltypes...)}, jval)
                        pi = ipi + jpi
                        for o in range
                            push!(ls, pi + logcpdf(gsf, Tuple(tjval), o) + get_log_score(lambda, o))
                        end
                    end
                end
                result[i] = StatsFuns.logsumexp(ls)
            end
        elseif parent_ix <= length(parranges)
            n = length(iranges)
            jix = parent_ix - n
                # We weight the lambda messages sent to the specific J parent by the pi message for the I parents
            result_pieces = [Vector{Float64}() for x in target_parrange]
            for iind in iinds
                ival = [iranges[r][iind[r]] for r in 1:n]
                gsf = gensf(sf,tuple(ival...))
                ipi = sum(logcpdf(ipis[r], (), iind[r]) for r in 1:n)
                lam = send_lambda(gsf, lambda, range, jranges, jpis, jix)
                for i in 1:length(target_parrange)
                    push!(result_pieces[i], ipi + get_log_score(lam, parranges[parent_ix][i]))
                end
            end
            result = [StatsFuns.logsumexp(rp) for rp in result_pieces]
        end

        return LogScore(target_parrange, result)
    end
end


text:
using AlphaZero

include("../main.jl")
using .ConnectFour

computer = ConnectFour.Solver.Player()
human = Human()

interactive!(GameSpec(), computer, human)


text:
using Documenter
using SpinMonteCarlo

makedocs( format=:html,
          sitename="SpinMonteCarlo.jl",
          modules=[SpinMonteCarlo],
          pages = [
                   "Home" => "index.md",
                   "Manual" => Any[
                                   "Run Monte Carlo" => "runmc.md",
                                   "Develop Monte Carlo" => "develop.md",
                                  ],
                   "Library" => Any[
                                    "Public" => "lib/public.md",
                                    "Internals" => "lib/internals.md",
                                   ]
                  ]
        )

deploydocs( repo = "github.com/yomichi/SpinMonteCarlo.jl.git",
            julia="0.7",
            target="build",
            deps = nothing,
            make = nothing,
          )


text:
# Tools for downloading test videos

module TestVideos

using VideoIO: VideoIO
import Base: show
import Downloads: download
export testvideo

videodir = joinpath(dirname(@__FILE__), "..", "videos")

mutable struct VideoFile{compression}
    name::AbstractString
    description::AbstractString
    license::AbstractString
    credit::AbstractString
    source::AbstractString
    download_url::AbstractString
    numframes::Int
    testframe::Int
    summarysize::Int
end

show(io::IO, v::VideoFile) = print(
    io,
    """
VideoFile:
   name:         $(v.name)  $(isfile(joinpath(videodir, v.name)) ? "(Downloaded)" : "")
   description:  $(v.description)
   license:      $(v.license)
   credit:       $(v.credit)
   source:       $(v.source)
   download_url: $(v.download_url)
   numframes:    $(v.numframes)
   summarysize:  $(v.summarysize)
 """,
)

VideoFile(name, description, license, credit, source, download_url, numframes, testframe, summarysize) =
    VideoFile{:raw}(name, description, license, credit, source, download_url, numframes, testframe, summarysize)

# Standard test videos
const videofiles = Dict(
    "ladybird.mp4" => VideoFile(
        "ladybird.mp4",
        "Ladybird opening wings (slow motion)",
        "Creative Commons: By Attribution 3.0 Unported (http://creativecommons.org/licenses/by/3.0/deed)",
        "CC-BY NatureClip (http://youtube.com/natureclip)",
        "https://downloadnatureclip.blogspot.com/p/download-links.html",
        "https://archive.org/download/LadybirdOpeningWingsCCBYNatureClip/Ladybird%20opening%20wings%20CC-BY%20NatureClip.mp4",
        397,
        13,
        3216,
    ),
    "annie_oakley.ogg" => VideoFile(
        "annie_oakley.ogg",
        "The \"Little Sure Shot\" of the \"Wild West,\" exhibition of rifle shooting at glass balls.",
        "pubic_domain (US)",
        "",
        "https://commons.wikimedia.org/wiki/File:Annie_Oakley_shooting_glass_balls,_1894.ogg",
        "https://upload.wikimedia.org/wikipedia/commons/8/87/Annie_Oakley_shooting_glass_balls%2C_1894.ogv",
        726,
        2,
        167311096,
    ),
    "crescent-moon.ogv" => VideoFile(
        "crescent-moon.ogv",
        "Moonset (time-lapse).",
        "Creative Commons Attribution 2.0 Generic (http://creativecommons.org/licenses/by/2.0/deed)",
        "Photo : Thomas Bresson",
        "https://commons.wikimedia.org/wiki/File:2010-10-10-Lune.ogv",
        "https://upload.wikimedia.org/wikipedia/commons/e/ef/2010-10-10-Lune.ogv",
        1213,
        1,
        9744,
    ),
    "black_hole.webm" => VideoFile(
        "black_hole.webm",
        "Artist’s impression of the black hole inside NGC 300 X-1 (ESO 1004c)",
        "Creative Commons Attribution 3.0 Unported (http://creativecommons.org/licenses/by/3.0/deed)",
        "Credit: ESO/L. Calçada",
        "https://www.eso.org/public/videos/eso1004a/",
        "https://upload.wikimedia.org/wikipedia/commons/1/13/Artist%E2%80%99s_impression_of_the_black_hole_inside_NGC_300_X-1_%28ESO_1004c%29.webm",
        597,
        1,
        4816,
    ),
)

"""
    testvideo(name, ops...)

Returns an `AVInput` object for the given video name. The video will be downloaded if it isn't available.
"""
function testvideo(name, ops...)
    videofile = joinpath(videodir, name)
    if !isfile(videofile)
        files = collect(keys(videofiles))
        # Bool is needed here for Julia v0.3
        ind = Bool[startswith(x, name) for x in files]
        count = sum(ind)

        if count == 1
            name = files[ind][1]
            videofile = joinpath(videodir, name)
            !isfile(videofile) && download(videofiles[name])
            !isfile(videofile) && throw(ErrorException("Error downloading test videofile \"$name\""))
        elseif count == 0
            throw(ErrorException("$name is not a known test videofile"))
        else
            throw(ErrorException("$name matches more than one test videofile"))
        end
    end
    return VideoIO.open(videofile, ops...) # Returns an AVInput, which in turn must be opened with openvideo, openaudio, etc.
end

function write_info(v::VideoFile)
    credit_file = joinpath(videodir, "$(v.name).txt")
    Base.open(credit_file, "w") do f
        return println(f, v)
    end
end

function download(v::VideoFile)
    write_info(v)
    println(stderr, "Downloading $(v.name) to $videodir")
    return download(v.download_url, joinpath(videodir, v.name))
end
"""
    available()

Print a list of all available test videos.
"""
function available()
    for v in values(videofiles)
        println(v)
    end
end

names() = collect(keys(videofiles))

"""
    download_all()

Downloads all test videos.
"""
function download_all()
    for (filename, v) in videofiles
        videofile = joinpath(videodir, filename)

        !isfile(videofile) && download(v)

        infofile = joinpath(videodir, "$(v.name).txt")
        !isfile(infofile) && write_info(v)
    end
end

"""
    remove_all()

Remove all test videos.
"""
function remove_all()
    for filename in keys(videofiles)
        videofile = joinpath(videodir, filename)
        isfile(videofile) && rm(videofile)
    end
end

end


text:
module _Plots
using Requires: @require

abstract type AbstractShimFunction{name} <: Function end
(::AbstractShimFunction)(_args...; _kwargs...) = error("please import Plots")

struct ShimFunction{name} <: AbstractShimFunction{name} end

macro shim(name::Symbol)
    esc(:(const $name = ShimFunction{$(QuoteNode(name))}()))
end

@shim grid

function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        (::ShimFunction{name})(args...; kwargs...) where {name} =
            getfield(Plots, name)(args...; kwargs...)
    end
end

end  # module _Plots

struct IntervalsPlotter
    stks::IntervalStacks
end

struct PointsPlotter
    stks::IntervalStacks
end

@recipe function _(stks::IntervalStacks)
    link := :x
    layout := _Plots.grid(2, 1; heights = [0.9, 0.1])

    @series begin
        subplot := 1
        IntervalsPlotter(stks)
    end

    @series begin
        subplot := 2
        yguide --> "Pt."
        PointsPlotter(stks)
    end
end

function time_converter(stks::IntervalStacks)
    t0 = time_ns()
    if !isempty(stks.trees)
        t0 = stks.events[stks.trees[1].eventid].start
    end
    astime(t) = (t - t0) / 1e9
    return astime
end

@recipe function _(pltr::IntervalsPlotter)
    stks = pltr.stks

    data = plot_data(stks)
    for (x, y, c) in zip(data.xs, data.ys, data.cs)
        # y = y .+ (rand() .- 0.5)
        @series begin
            seriestype := :line
            x := x
            y := y
            fillrange := y[1] + 1
            fillcolor := c
            linecolor := :white
            # linestyle := :dot
            # linewidth := 4
            label := nothing
            ()
        end
    end
    @series begin
        seriestype := :scatter
        x := data.root.xs
        y := data.root.ys
        # markersize := 2
        # markershape := :x
        markercolor := :white
        markerstrokecolor := :black
        label := nothing
        ()
    end

    @series begin
        yguide --> "Intervals"
        label := nothing
        ()
    end
end

@recipe function _(pltr::PointsPlotter)
    stks = pltr.stks
    astime = time_converter(stks)
    times = [astime(stks.events[i].start) for i in stks.points]

    seriestype := :vline
    y := times
    xguide --> "Time [Sec]"
    yguide --> "Points"
    yticks := nothing
    label --> nothing
    ()
end

function plot_data(stks::IntervalStacks)
    astime = time_converter(stks)

    # TODO: use matrix?
    xs = Vector{Float64}[]
    ys = Vector{Float64}[]
    cs = Int[]
    data = (; xs = xs, ys = ys, cs = cs, root = (; xs = Float64[], ys = Float64[]))

    colormap = Dict{UInt64,Int}()
    function getcolor(event)
        # k = hash(event.approx_taskid, event.locationid)
        # k = event.locationid
        k = event.approx_taskid
        get!(colormap, k) do
            length(colormap) + 1
        end
    end

    function add(tree, y0)
        event = stks.events[tree.eventid]
        color = getcolor(event)
        start = astime(event.start)
        stop = astime(event.stop)
        push!(xs, [start, stop])
        push!(ys, [y0, y0])
        push!(cs, color)
        tree.branches === nothing && return
        for b in tree.branches
            add(b, y0 + 1)
        end
    end

    # Edges of the "staircase":
    rights = TimeNS[typemin(TimeNS)]
    tops = Int[0]
    # h = 0
    for tree in stks.trees
        # add(tree, h)
        # h += tree.height + 2
        # continue

        # @show rights, tops
        event = stks.events[tree.eventid]
        i = searchsortedlast(rights, event.start)
        if tree.height + 2 < tops[i]
            y0 = 0
            j = searchsortedfirst(rights, event.stop)
            if j > length(rights)
                tops[j-1] = tree.height + 2
            end
            rights[j-1] = event.stop
            deleteat!(rights, 2:j-2)
            deleteat!(tops, 2:j-2)
            # @show isnothing(tree.branches), astime(event.start), astime(event.stop)
        else
            y0 = tops[end]
            if event.stop < rights[end]
                tops[end] += tree.height + 2
            else
                push!(rights, event.stop)
                push!(tops, tops[end] + tree.height + 2)
            end
        end
        add(tree, y0)
        push!(data.root.xs, astime(event.start))
        push!(data.root.ys, y0)
    end

    return data
end


text:
#
# Common data types and associated functions
#

# A custom `Unitful` dimension & unit for mass parameters.
Unitful.@derived_dimension MassParameter Unitful.𝐋^3/Unitful.𝐓^2

"""
Absract type for orbital coordinate frames.
"""
abstract type AbstractFrame end

# TODO // NOTE a future Julia patch will clarify how types and associated values are parameterized
"""
Abstract type for structures parameterized (in part) by `Unitful.Unit` types.
"""
abstract type AbstractUnitfulStructure{F <: AbstractFloat, LU <: Unitful.LengthFreeUnits, TU <: Unitful.TimeFreeUnits} end

"""
Abstract type for all structures describing orbital states (typically Cartesian or Keplerian).
"""
abstract type AbstractState{F, LU, TU, FR<:AbstractFrame} <: AbstractUnitfulStructure{F, LU, TU} end

"""
Abstract type for all structures describing orbital systems (R2BP systems, CR3BP systems, NBP systems, etc.)
"""
abstract type AbstractSystem{F, LU, TU} <: AbstractUnitfulStructure{F, LU, TU} end

"""
Abstract type for all structures which hold __both__ states and systems (aka the whole `Orbit`!)
"""
abstract type AbstractOrbit{F, LU, TU} <: AbstractUnitfulStructure{F, LU, TU} end

"""
An inertial coordinate frame, centered at _some_ location (typically the center of a body, or a barycenter).
"""
struct Inertial <: AbstractFrame end

"""
A Synodic coordinate frame. Located at the barycenter of a CR3BP system, and rotates with the two Celestial Bodies.
"""
struct Synodic <: AbstractFrame end

"""
A Perifocal coordinate frame. An `Inertial` frame which is 2D (orbit falls completely within the XY axis).
"""
struct Perifocal <: AbstractFrame end

"""
Cartesian state which describes a spacecraft or body's position and velocity with respect to _something_,
as described the the frame `FR`.
"""
mutable struct CartesianState{F, LU, TU, FR} <: AbstractState{F, LU, TU, FR}
    t::F
    r::SubArray{F, 1, MVector{6, F}, Tuple{UnitRange{Int64}}, true}
    v::SubArray{F, 1, MVector{6, F}, Tuple{UnitRange{Int64}}, true}
    rv::MVector{6,F}

    function CartesianState(r::AbstractVector{R}, v::AbstractVector{V}, epoch::E=0, frame=Bodycentric;
                            lengthunit=u"km", timeunit=u"s") where {R<:Real, V<:Real, E<:Real}
        F  = promote_type(R, V, E)
        if !(F <: AbstractFloat)
            @warn "Type provided ($(string(F))) is not a float: defaulting to Float64."
            F = Float64
        end
        @assert length(r) == length(v) == 3 "Both arguments must have length equal to 3!"
        rv = MVector{6, F}(r[1], r[2], r[3], v[1], v[2], v[3])
        pos = @views rv[1:3]
        vel = @views rv[4:6]
        return new{F, typeof(lengthunit), typeof(timeunit), frame}(F(epoch), pos, vel, rv)
    end

    function CartesianState(r::AbstractVector{R}, v::AbstractVector{V},
                            epoch::T = 0u"s", frame=Bodycentric) where {R<:Unitful.Length, V<:Unitful.Velocity, T<:Unitful.Time}

        # Get ready to commit a crime... we need the Time unit provided in velocity quantity V

        # V is a quantity with units Length / Time.  Let's reverse that so we have Time instead of Time^(-1)
        TL = inv(unit(V))

        # Now we need to find which index (1 or 2) contains the Time unit
        timeaxis = findfirst(T -> T isa Unitful.Dimension{:Time}, collect(typeof(typeof(TL).parameters[2]).parameters[1]))

        # Now that we have the proper index, let's select the time unit
        timeunit = Unitful.FreeUnits{(typeof(TL).parameters[1][timeaxis],), Unitful.𝐓, nothing}()

        # This is easy...
        lengthunit = unit(R)

        # Phew!
        return CartesianState(ustrip.(lengthunit, r), ustrip.(lengthunit/timeunit, v),
                              ustrip(timeunit, epoch), frame; lengthunit = lengthunit, timeunit = timeunit)
    end

    function CartesianState(cart::CartesianState{F,LU,TU,FR}) where {F,LU,TU,FR}
        return CartesianState(cart.r,  cart.v, cart.t, FR; lengthunit = LU(), timeunit = TU())
    end

    function CartesianState(arr::StaticVector{6})
        return CartesianState(arr[1:3], arr[4:6])
    end
end

"""
Convert `CartesianState` values between different `eltype` and `lengthunit` and `timeunit` types.
"""
function Base.convert(::Type{CartesianState{F,LU,TU}}, cart::CartesianState) where {F,LU,TU}
    r = F.(ustrip.(LU(), position_vector(cart)))
    v = F.(ustrip.(LU()/TU(), velocity_vector(cart)))
    t = F.(ustrip.(TU(), epoch(cart)))
    return CartesianState(r, v, t, coordinateframe(cart); lengthunit = LU(), timeunit = TU())
end

"""
Convert `CartesianState` values between the same `eltype` and `lengthunit` and `timeunit` types.
"""
function Base.convert(::Type{CartesianState{F,LU,TU}}, cart::CartesianState{F,LU,TU}) where {F,LU,TU}
    return cart
end

# Overrides `Unitful` unit conversions for `AbstractUnitfulStructure` instances.
(u::Unitful.LengthFreeUnits)(state::CartesianState{F,LU,TU}) where {F,LU,TU} = convert(CartesianState{F, typeof(u), TU}, state)

# Overrides `Unitful` unit conversions for `AbstractUnitfulStructure` instances.
(u::Unitful.TimeFreeUnits)(state::CartesianState{F,LU,TU}) where {F,LU,TU} = convert(CartesianState{F, LU, typeof(u)}, state)

"""
Print a `CartesianState` to `io`.
"""
function Base.show(io::IO, state::CartesianState{F,LU,TU,FR}) where {F,LU,TU,FR}
    println(io, "  ", string(FR), " Cartesian State:")
    println("")
    println(io, "    t =  ", round(state.t; digits=6), " ", string(TU()))
    println(io, "    r = ", [round(state.r[1]; digits=6) round(state.r[2]; digits=6) round(state.r[3]; digits=6)], " ", string(LU()))
    println(io, "    v = ", [round(state.v[1]; digits=6) round(state.v[2]; digits=6) round(state.v[3]; digits=6)], " ", string(LU()/TU()))
end

"""
Print a `CartesianState` to `io`.
"""
Base.show(io::IO, ::MIME"text/plain", state::CartesianState{F,LU,TU,FR}) where {F,LU,TU,FR} = show(io, state)

"""
A normalized length unit, typically used in normalized CR3B problems.
"""
const NormalizedLengthUnit = Unitful.FreeUnits{(), Unitful.𝐋, nothing}

"""
A normalized time unit, typically used in normalized CR3B problems.
"""
const NormalizedTimeUnit = Unitful.FreeUnits{(), Unitful.𝐓, nothing}

"""
Returns the timepoint associated with a `CartesianState`.
"""
epoch(cart::CartesianState) = cart.t * timeunit(cart)

"""
Returns the first paremeter (commonly labeled as `F`) of a `AbstractUnitfulStructure`.
"""
Base.eltype(::AbstractUnitfulStructure{F}) where F = F

"""
Returns the `Unitful.Length` unit associated with a structure.
"""
lengthunit(::AbstractUnitfulStructure{F, LU, TU}) where {F,LU,TU} = LU()

"""
Returns the `Unitful.Time` unit associated with a structure.
"""
timeunit(::AbstractUnitfulStructure{F, LU, TU}) where {F, LU, TU} = TU()

"""
Returns the `Unitful.Velocity` unit associated with a structure.
"""
velocityunit(state::AbstractUnitfulStructure) = lengthunit(state) / timeunit(state)

"""
Returns the `MassParameter` unit associated with a structure.
"""
massparameterunit(state::AbstractUnitfulStructure) = lengthunit(state)^3 / timeunit(state)^2

"""
Returns the coordinate frame associated with a state.
"""
coordinateframe(::AbstractState{F, LU, TU, FR}) where {F, LU, TU, FR} = FR

"""
Returns the `Unitful` position vector of the `Cartesianstate`.
"""
position_vector(state::CartesianState{F,LU,TU}) where {F, LU, TU} = state.r * lengthunit(state)

"""
Returns the `Unitful` velocity vector of the `CartesianState`.
"""
velocity_vector(state::CartesianState{F,LU,TU}) where {F, LU, TU} = state.v * velocityunit(state)

"""
Returns the `Unitful` scalar position of the `CartesianState`.
"""
scalar_position(state::CartesianState{F,LU,TU}) where {F, LU, TU} = norm(position_vector(state))

"""
Returns the `Unitful` scalar velocity of the `CartesianState`.
"""
scalar_velocity(state::CartesianState{F,LU,TU}) where {F, LU, TU} = norm(velocity_vector(state))

"""
A `Trajectory` is just a vector of `Orbit` instances!
"""
const Trajectory{T<:AbstractOrbit} = Vector{T}

"""
Prints a `Trajectory` instance to `io`.
"""
Base.show(io::IO, ::MIME"text/plain", traj::Trajectory) = show(io,traj)

"""
A `Manifold` is a vector of `Trajectory` instances.
"""
const Manifold{T} = Vector{Trajectory{T}} where T <: AbstractOrbit

"""
An alias for `CartesianState` instances with `Synodic` coordinate frames.
"""
const SynodicCartesianState{F, LU, TU} = CartesianState{F, LU, TU, Synodic}

"""
An alias for `CartesianState` instances with `Inertial` coordinate frames.
"""
const InertialCartesianState{F, LU, TU} = CartesianState{F, LU, TU, Inertial}

"""
A common reference frame for astrodynamics problems that humans solve!
The Earth-centered Inertial frame.
"""
struct ECI <: AbstractFrame end;

"""
A common reference frame for astrodynamics problems that humans solve!
The Sun-centered (heliocentric) Inertial frame.
"""
struct HCI <: AbstractFrame end;

"""
An inertial frame centered at the CR3BP system barycenter.
"""
struct BarycentricInertial  <: AbstractFrame end;

Base.string(::Type{Inertial}) = "Inertial"
Base.string(::Type{BarycentricInertial}) = "BarycentricInertial"
Base.string(::Type{ECI}) = "ECI"
Base.string(::Type{HCI}) = "HCI"
Base.string(::Type{Synodic}) = "Synodic"


text:
using Test
using WordTokenizers

@testset "Tweet Tokenize" begin
    @testset "Basic Tests" begin
        s0 = "This is a cooool #dummysmiley: :-) :-P <3 and some arrows < > -> <--"
        @test tweet_tokenize(s0) ==
            ["This", "is", "a", "cooool", "#dummysmiley", ":", ":-)", ":-P", "<3", "and", "some", "arrows", "<", ">", "->", "<--"]

        s1 = "@Joyster2012 @CathStaincliffe Good for you, girl!! Best wishes :-)"
        @test tweet_tokenize(s1) ==
            ["@Joyster2012", "@CathStaincliffe", "Good", "for", "you", ",", "girl", "!", "!", "Best", "wishes", ":-)"]

        s2 = "3Points for #DreamTeam Gooo BAILEY! :) #PBB737Gold @PBBabscbn"
        @test tweet_tokenize(s2) ==
            ["3Points", "for", "#DreamTeam", "Gooo", "BAILEY", "!", ":)", "#PBB737Gold", "@PBBabscbn"]

        s3 = "@Insanomania They do... Their mentality doesn't :("
        @test tweet_tokenize(s3) ==
            ["@Insanomania", "They", "do", "...", "Their", "mentality", "doesn't", ":("]

        s4 = "RT @facugambande: Ya por arrancar a grabar !!! #TirenTirenTiren vamoo !!"
        @test tweet_tokenize(s4) ==
            ["RT", "@facugambande", ":", "Ya", "por", "arrancar", "a", "grabar", "!", "!", "!", "#TirenTirenTiren", "vamoo", "!", "!"]

        s5 = "@crushinghes the summer holidays are great but I'm so bored already :("
        @test tweet_tokenize(s5, reduce_len=true) ==
            ["@crushinghes", "the", "summer", "holidays", "are", "great", "but", "I'm", "so", "bored", "already", ":("]

        s6 = "@jrmy: I'm REALLY HAPPYYY about that! NICEEEE :D :P"
        @test tweet_tokenize(s6) ==
            ["@jrmy", ":", "I'm", "REALLY", "HAPPYYY", "about", "that", "!", "NICEEEE", ":D", ":P"]
    end

    @testset "Remove Handles and Reduce Length" begin
        s7 = "@remy: This is waaaaayyyy too much for you!!!!!!"
        @test tweet_tokenize(s7, strip_handle=true, reduce_len=true) ==
            [":", "This", "is", "waaayyy", "too", "much", "for", "you", "!", "!", "!"]

        s8 = "@_willy65: No place for @chuck tonight. Sorry."
        @test tweet_tokenize(s8, strip_handle=true, reduce_len=true) ==
            [":", "No", "place", "for", "tonight", ".", "Sorry", "."]

        s9 = "@mar_tin is a great developer. Contact him at mar_tin@email.com."
        @test tweet_tokenize(s9, strip_handle=true, reduce_len=true) ==
            ["is", "a", "great", "developer", ".", "Contact", "him", "at", "mar_tin@email.com", "."]
       end

    @testset "Test long sentences" begin
        s10 = "Photo: Aujourd'hui sur http://t.co/0gebOFDUzn Projet... http://t.co/bKfIUbydz2.............................. http://fb.me/3b6uXpz0L"
        @test tweet_tokenize(s10) ==
            ["Photo", ":", "Aujourd'hui", "sur", "http://t.co/0gebOFDUzn", "Projet", "...", "http://t.co/bKfIUbydz2", "...", "http://fb.me/3b6uXpz0L"]
    end

    @testset "Replace HTML Entities" begin
        @test tweet_tokenize("An HTML Entity - Δ") ==
        ["An", "HTML", "Entity", "-", "Δ"]

        @test tweet_tokenize("Another HTML Entity - Δ") ==
        ["Another", "HTML", "Entity", "-", "Δ"]

        @test tweet_tokenize("Another HTML Entity - Δ") ==
        ["Another", "HTML", "Entity", "-", "Δ"]

        @test tweet_tokenize("Price: £100") ==
        [ "Price", ":", "£", "100"]

        @test tweet_tokenize("Check out this invalid symbol �") ==
        [ "Check", "out", "this", "invalid", "symbol", "\u81"]

        @test tweet_tokenize("A•B = B•A ") ==
        [ "A", "•", "B", "=", "B", "•", "A"]

        @test tweet_tokenize("Check out this symbol in Windows-1252 encoding €") ==
        [ "Check", "out", "this", "symbol", "in", "Windows", "-", "1252", "encoding", "€"]
    end

    @testset "Token Buffer functions" begin
        # Emoticons reverse and HTML_Entities
        @test tweet_tokenize("(-:> (-: <head> Hi there! </head>") ==
        ["(-:>", "(-:", "<head>", "Hi", "there", "!", "</head>"]

        # Phone numbers and email
        @test tweet_tokenize("+0 (123) 333-5553 and 1111222 are valid for ab@cd.e.") ==
        ["+0 (123) 333-5553", "and", "1111222", "are", "valid", "for", "ab", "@cd", ".", "e", "."]

        @test tweet_tokenize("+0      0------      333 333 2222 333 33 3 33333") ==
        ["+0", "0---   333 333 2222","333", "33", "3", "33333"]

        @test tweet_tokenize("11111112222 11112222") ==
        ["11111112222", "11112222"]

        # Hashtags and Twitter usernames
        @test tweet_tokenize("#a- @ ") ==
        ["#", "a-", "@"]

        # URLs
        @test tweet_tokenize("abc.com/(xyz)(xyz)a(a( and abc.co.com/ and http:  ") ==
        ["abc.com/(xyz)(xyz)a", "(", "a", "(", "and", "abc.co.com/", "and", "http", ":"]

    end
end


text:
using jInv.Mesh
using jInv.LinearSolvers
using Multigrid.ParallelJuliaSolver
using SparseArrays
using LinearAlgebra
using Test

speye = (n) -> SparseMatrixCSC(1.0I,n,n);

println("===  Test Julia wrapper with symmetric (real) matrix ====");
domain = [0.0, 1.0, 0.0, 1.0];
n      = [20,23];

Mr     = getRegularMesh(domain,n)
G      = getNodalGradientMatrix(Mr);
m      = sparse(Diagonal(exp.(randn(size(G,1)))));
Ar     = G'*m*G;
Ar     = Ar + 1e-1*norm(Ar,1)*speye(size(Ar,2));
N      = size(Ar,2);
b      = Ar*rand(N);
bs     = Ar*rand(N,5);
Bs     = (convert(Array{Float64,1}, b),convert(Array{Float64,2}, bs));

LU = getParallelJuliaSolver(Float64,Int64,numCores=1,backend=1);

for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	t = @elapsed local x,  = solveLinearSystem(Ar,Bs[j],LU);
	@test norm(Ar*x-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-8
	#println("V0: took ",t);
	println("LU solve time: ",LU.solveTime);
	LU.solveTime=0.0;
end

clear!(LU);
LU = getParallelJuliaSolver(Float64,Int64,numCores=1,backend=2);

for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	t = @elapsed local x,  = solveLinearSystem(Ar,Bs[j],LU);
	@test norm(Ar*x-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-8
	#println("V1: took ",t);
	println("LU solve time: ",LU.solveTime);
	LU.solveTime=0.0;
end

LU = getParallelJuliaSolver(Float64,Int64,numCores=2,backend=3);

for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	t = @elapsed local x,  = solveLinearSystem(Ar,Bs[j],LU);
	@test norm(Ar*x-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-8
	#println("V1: took ",t);
	println("LU solve time: ",LU.solveTime);
	LU.solveTime=0.0;
end

Bs     = (convert(Array{Float32,1}, b),convert(Array{Float32,2}, bs));
LU = getParallelJuliaSolver(Float32,UInt32,numCores=2,backend=3);

for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	t = @elapsed local x,  = solveLinearSystem(Ar,Bs[j],LU);
	@test norm(Ar*x-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-4
	#println("V1: took ",t);
	println("LU solve time: ",LU.solveTime);
	LU.solveTime=0.0;
end

println("===  Test Julia wrapper with shifted (complex) Laplacian ====");
Ar     = G'*m*G + (1+1im)*speye(size(Ar,2));
b      = Ar*rand(ComplexF64,N);
bs     = Ar*rand(ComplexF64,N,5);
Bs     = (b,bs);

LU = copySolver(LU);
clear!(LU);

LU = getParallelJuliaSolver(ComplexF64,Int64,numCores=2,backend=3);

for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	local x,  = solveLinearSystem(Ar,Bs[j],LU);
	@test norm(Ar*x-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-10
	print("\n")
end
println("")
Bs     = (convert(Array{ComplexF32,1}, b),convert(Array{ComplexF32,2}, bs));
LU = getParallelJuliaSolver(ComplexF32,UInt32,numCores=2,backend=3);
for j=1:length(Bs)
	print("nrhs=$(size(Bs[j],2)), issparse(rhs)=$(issparse(Bs[j])) : ")
	local x,  = solveLinearSystem(Ar,Bs[j],LU);
	@test norm(Ar*x-Bs[j],Inf)/norm(Bs[j],Inf) < 1e-4
	print("\n")
end
println("")

println("===  Test Parellel Julia Wrapper: nonsymmetric matrices ====");
n = 50
sNonSym  = getParallelJuliaSolver(Float64,Int64,numCores=1,backend=1);
sNonSym  = copySolver(sNonSym);
A = sprandn(n,n,5/n) + 10*speye(n);
B = randn(Float64,n)
Bs = randn(Float64,n,5)

X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-5

X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-5

X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-5

X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-5

X, = solveLinearSystem(A,Bs,sNonSym,1);
@test norm(A'*X - Bs)/norm(Bs) < 1e-5

println("===  Testing parallel backend with transposed  ===")
sNonSym  = getParallelJuliaSolver(Float64,Int64,numCores=2,backend=3);
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-8
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-8
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-8
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-8
X, = solveLinearSystem(A,Bs,sNonSym,1);
@test norm(A'*X - Bs)/norm(Bs) < 1e-8

println("===  Testing parallel backend with transposed  ===")
sNonSym  = getParallelJuliaSolver(Float32,UInt32,numCores=2,backend=3);
B = randn(Float32,n)
Bs = randn(Float32,n,5)
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,Bs,sNonSym,1);
@test norm(A'*X - Bs)/norm(Bs) < 1e-4

println("===  Testing parallel backend with complex transposed 64bit  ===")
sNonSym  = getParallelJuliaSolver(ComplexF64,Int64,numCores=2,backend=3);
B = randn(ComplexF64,n)
Bs = randn(ComplexF64,n,5)
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-5
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-5
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-5
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-5
X, = solveLinearSystem(A,Bs,sNonSym,1);
@test norm(A'*X - Bs)/norm(Bs) < 1e-5

println("===  Testing parallel backend with complex transposed 32 bit ===")
sNonSym  = getParallelJuliaSolver(ComplexF32,UInt32,numCores=2,backend=3);
B = randn(ComplexF32,n)
Bs = randn(ComplexF32,n,5)
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,B,sNonSym,1);
@test norm(A'*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,B,sNonSym,0);
@test norm(A*X - B)/norm(B) < 1e-4
X, = solveLinearSystem(A,Bs,sNonSym,1);
@test norm(A'*X - Bs)/norm(Bs) < 1e-4

println("===  End Test Parallel Julia Solver ====");


text:
function _read(cmd)
    out = Pipe()
    if VERSION > v"1.1-"
      procs = Base._spawn(cmd, Any[devnull,out,out])
    else
      procs = Base._spawn(cmd, (devnull,out,out))
    end
    close(out.in)
    return String(read(out)), success(procs)
end

function initialise(dir)
  (isfile(joinpath(dir, "Project.toml")) && isdir(joinpath(dir, "src"))) ||
    error("No Julia project found at $dir")
  tmp = joinpath(tempdir(), "vimes-$(rand(UInt64))")
  mkdir(tmp)
  for path in readdir(dir)
      if !startswith(path, ".")
          cp(joinpath(dir, path), joinpath(tmp, path))
      end
  end
  atexit(() -> isdir(tmp) && rm(tmp, recursive=true))
  return tmp
end

function runtests(dir)
  out, pass = _read(`$(Base.julia_cmd()) --project=$dir -e 'using Pkg; Pkg.test()'`)
end

function mutate(dir, idx)
  pick!(idx)
end

function reset(dir, tmp)
  rm(joinpath(tmp, "src"), recursive=true)
  cp(joinpath(dir, "src"), joinpath(tmp, "src"))
  return
end

diff(dir, tmp) = _read(`git diff $(joinpath(dir, "src")) $(joinpath(tmp, "src"))`)[1]

function logdiff(dir, tmp, dead = false)
  d = !dead ? ".vimes" : ".vimes-dead"
  mkpath(joinpath(dir, d))
  write(joinpath(dir, d, "$(now()).diff"), diff(dir, tmp))
end

function checktests(dir, tmp)
  out, pass = runtests(tmp)
  if !pass
    mkpath(joinpath(dir, ".vimes"))
    write(joinpath(dir, ".vimes", "log.txt"), out)
    error("Tests fail before starting. See `.vimes/log.txt`.")
  end
end

function test(dir, tmp, idx; dead = false)
  while isempty(diff(dir, tmp))
    mutate(tmp, idx)
  end
  _, pass = runtests(tmp)
  pass && logdiff(dir, tmp)
  !pass && dead && logdiff(dir, tmp, true)
  reset(dir, tmp)
  return !pass
end

function go(dir, ps = defaults; procs = 1, dead = false, tests = typemax(Int))::Int
  runs, pass = 0, 0
  function run(r, i)
    runs += 1
    pass += r
    @info("($i) Ran $runs tests, precision $(@sprintf("%.2f", pass/runs*100))%")
  end
  tmp = initialise(dir)
  idx = indices(joinpath(tmp, "src"), ps)
  checktests(dir, tmp)
  rm(tmp, recursive=true)
  @sync for i = 1:procs
    let tmp = initialise(dir), idx = indices(joinpath(tmp, "src"), ps)
      @async while runs < tests
        run(test(dir, tmp, idx, dead = dead), i)
      end
    end
  end
  return pass
end


text:
module NonlinearExtendedTest

using Test
using LinearAlgebra
using ForneyLab
using ForneyLab: outboundType, isApplicable, Extended, requiresBreaker, breakerParameters
using ForneyLab: SPNonlinearEOutNG, SPNonlinearEIn1GG, SPNonlinearEOutNGX, SPNonlinearEInGX, MNonlinearEInGX
using ForneyLab: concatenate, localLinearization, requiresBreaker, breakerParameters

f(x) = x

g(x::Number) = x^2 - 5.0
g(x::Vector) = x.^2 .- 5.0
g_inv(y::Number) = sqrt(y + 5.0)
g_inv(y::Vector) = sqrt.(y .+ 5.0)

h(x::Number, y::Number) = x^2 - y
h(x::Vector, y::Vector) = x.^2 .- y
h_inv_x(z::Number, y::Number) = sqrt(z + y)
h_inv_x(z::Vector, y::Vector) = sqrt.(z .+ y)

#--------
# Helpers
#--------

@testset "concatenate" begin
    @test concatenate([[1.0, 2.0], [3.0]]) == ([1.0, 2.0, 3.0], [2, 1])
end

@testset "localLinearization" begin
    (a, b) = localLinearization(Univariate, g, 1.0)
    @test a == 2.0
    @test b == -6.0

    (A, b) = localLinearization(Multivariate, g, [1.0])
    @test A == mat(2.0)
    @test b == [-6.0]

    (A, b) = localLinearization(Univariate, h, [1.0, 2.0])
    @test A == [2.0, -1.0]'
    @test b == -1.0

    (A, b) = localLinearization(Multivariate, h, [[1.0], [2.0]])
    @test A == [2.0 -1.0]
    @test b == [-1.0]
end

@testset "requiresBreaker and breakerParameters" begin
    # Without given inverse
    fg = FactorGraph()
    x = Variable()
    y = Variable()
    nd = GaussianMeanVariance(x, 0.0, 1.0)
    Nonlinear{Extended}(y, x, g=g)

    @test requiresBreaker(nd.i[:out])
    @test_throws Exception breakerParameters(nd.i[:out].partner)
    @test breakerParameters(nd.i[:out]) == (Message{GaussianMeanVariance, Univariate}, ())

    # With given inverse
    fg = FactorGraph()
    x = Variable()
    y = Variable()
    nd = GaussianMeanVariance(x, 0.0, 1.0)
    Nonlinear{Extended}(y, x, g=g, g_inv=g)

    @test !requiresBreaker(nd.i[:out])
    @test_throws Exception breakerParameters(nd.i[:out])
end

#-------------
# Update rules
#-------------

@testset "SPNonlinearEOutNG" begin
    @test SPNonlinearEOutNG <: SumProductRule{Nonlinear{Extended}}
    @test outboundType(SPNonlinearEOutNG) == Message{GaussianMeanVariance}
    @test isApplicable(SPNonlinearEOutNG, [Nothing, Message{GaussianMeanVariance}])

    @test ruleSPNonlinearEOutNG(g, nothing, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0)) == Message(Univariate, GaussianMeanVariance, m=-1.0, v=48.0)
    @test ruleSPNonlinearEOutNG(g, nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0))) == Message(Multivariate, GaussianMeanVariance, m=[-1.0], v=mat(48.0))
end

@testset "SPNonlinearEOutNGX" begin
    @test SPNonlinearEOutNGX <: SumProductRule{Nonlinear{Extended}}
    @test outboundType(SPNonlinearEOutNGX) == Message{GaussianMeanVariance}
    @test !isApplicable(SPNonlinearEOutNGX, [Nothing, Message{Gaussian}])
    @test isApplicable(SPNonlinearEOutNGX, [Nothing, Message{Gaussian}, Message{Gaussian}])
    @test !isApplicable(SPNonlinearEOutNGX, [Message{Gaussian}, Nothing, Message{Gaussian}])

    @test ruleSPNonlinearEOutNGX(h, nothing, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == Message(Univariate, GaussianMeanVariance, m=-1.0, v=49.0)
    @test ruleSPNonlinearEOutNGX(h, nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == Message(Multivariate, GaussianMeanVariance, m=[-1.0], v=mat(49.0))
end

@testset "SPNonlinearEIn1GG" begin
    @test SPNonlinearEIn1GG <: SumProductRule{Nonlinear{Extended}}
    @test outboundType(SPNonlinearEIn1GG) == Message{Gaussian}
    @test isApplicable(SPNonlinearEIn1GG, [Message{Gaussian}, Nothing])

    # Without given inverse
    @test ruleSPNonlinearEIn1GG(g, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=14.666666666666666, w=5.333333333333333)
    @test ruleSPNonlinearEIn1GG(g, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(1.0))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[14.666666666666671], w=mat(5.333333333333335))

    # With given inverse
    @test ruleSPNonlinearEIn1GG(g, g_inv, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), nothing) == Message(Univariate, GaussianMeanVariance, m=2.6457513110645907, v=0.10714285714285711)
    @test ruleSPNonlinearEIn1GG(g, g_inv, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), nothing) == Message(Multivariate, GaussianMeanVariance, m=[2.6457513110645907], v=mat(0.10714285714285711))
end

@testset "SPNonlinearEInGX" begin
    @test SPNonlinearEInGX <: SumProductRule{Nonlinear{Extended}}
    @test outboundType(SPNonlinearEInGX) == Message{Gaussian}
    @test !isApplicable(SPNonlinearEInGX, [Message{Gaussian}, Nothing])
    @test !isApplicable(SPNonlinearEInGX, [Nothing, Message{Gaussian}, Message{Gaussian}])
    @test isApplicable(SPNonlinearEInGX, [Message{Gaussian}, Nothing, Message{Gaussian}])

    # Without given inverse
    @test ruleSPNonlinearEInGX(h, 1, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0), Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=10.999999999999996, w=3.9999999999999982)
    @test ruleSPNonlinearEInGX(h, 1, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(1.0)), Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[10.999999999999996], w=mat(3.9999999999999982))

    # With given inverse
    @test ruleSPNonlinearEInGX(h, h_inv_x, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), nothing, Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == Message(Univariate, GaussianMeanVariance, m=1.3228756555322954, v=0.14285714285714282)
    @test ruleSPNonlinearEInGX(h, h_inv_x, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), nothing, Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == Message(Multivariate, GaussianMeanVariance, m=[1.3228756555322954], v=mat(0.14285714285714282))
end

@testset "MNonlinearEInGX" begin
    @test MNonlinearEInGX <: MarginalRule{Nonlinear{Extended}}
    @test isApplicable(MNonlinearEInGX, [Nothing, Message{Gaussian}, Message{Gaussian}])

    @test ruleMNonlinearEInGX(h, Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0), Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0), Message(Univariate, GaussianMeanVariance, m=5.0, v=1.0)) == ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[2.6, 4.85], v=[0.20000000000000007 0.19999999999999998; 0.19999999999999998 0.95])
    @test ruleMNonlinearEInGX(h, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(3.0)), Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(1.0)), Message(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(1.0))) == ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[2.6, 4.85], v=[0.20000000000000007 0.19999999999999998; 0.19999999999999998 0.95])
end

#------------
# Integration
#------------

@testset "Nonlinear integration via local linear approximation with given inverse" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    n = Nonlinear{Extended}(y, x, g=g, g_inv=g_inv)

    @test isa(n, Nonlinear{Extended})

    # Forward; g_inv should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearEOutNG(g, nothing, messages[2])", code)
    @test !occursin("g_inv", code)

    # Backward; g_inv should be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearEIn1GG(g, g_inv, messages[2], nothing)", code)
end

@testset "Multi-argument nonlinear integration via local linear approximation" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    @RV z ~ GaussianMeanVariance(5.0, 1.0)
    n = Nonlinear{Extended}(y, x, z, g=h, g_inv=[h_inv_x, nothing])

    # Forward; h_inv_x should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearEOutNGX(h, nothing, messages[3], messages[1])", code)
    @test !occursin("h_inv_x", code)

    # Backward with given inverse; h_inv_x should be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearEInGX(h, h_inv_x, messages[3], nothing, messages[1])", code)

    # Backward without given inverse
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(z)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearEInGX(h, 2, messages[3], messages[2], messages[1])", code)
    @test !occursin("h_inv_x", code)
    @test occursin("messages[1] = Message(vague(GaussianMeanVariance))", code)
end

@testset "Nonlinear integration via local linear approximation without given inverse" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    n = Nonlinear{Extended}(y, x, g=g)

    # Forward; g_inv should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearEOutNG(g, nothing, messages[1])", code)
    @test !occursin("g_inv", code)

    # Backward; g_inv should not be present in call,
    # both messages should be required, and initialization should take place
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearEIn1GG(g, messages[2], messages[1])", code)
    @test !occursin("g_inv", code)
    @test occursin("messages[1] = Message(vague(GaussianMeanVariance))", code)
end

end #module


text:
using MultiResolutionIterators
using Test

names = ["basic_functions", "functionality", "namedlevels", "demoscript", "customisability"]

@testset "$name" for name in names
    include("test_" * name *".jl")
end


text:
using DataWrangler
using Test
using Random
using Logging

const tests = [
    "boxcox",
    "d",
    "impute",
    "normality",
    "normalize",
    "outlie",
    "p",
]

printstyled("\nTest Summary List:\n", color=:underline)

Random.seed!(36)
Base.disable_logging(Base.CoreLogging.Error) # disable warnings

for t in tests
    @testset "Test $t" begin
        Random.seed!(36)
        include("$t.jl")
        println()
    end
end

println()


text:
function mountaincar_finite_model_class_main(; true_model = false, local_agent = false)
    finite_model_class_steps = []
    model = MountainCar(0.0)
    planners = Vector{MountainCarRTAAPlanner}()
    push!(planners, MountainCarRTAAPlanner(model, 1000, true_params))
    push!(
        planners,
        MountainCarRTAAPlanner(model, 1000, MountainCarParameters(-0.0025, 2.95)),
    )
    push!(planners, MountainCarRTAAPlanner(model, 1000, MountainCarParameters(-0.0024, 3)))
    push!(
        planners,
        MountainCarRTAAPlanner(model, 1000, MountainCarParameters(-0.0026, 3)),
    )
    push!(planners, MountainCarRTAAPlanner(model, 1000, MountainCarParameters(-0.0025, 3.05)))
    generateHeuristic!(planners)
    for rock_c in range_of_values
        println("Rock_c is ", rock_c)
        mountaincar = MountainCar(rock_c)
        if true_model
            planner = MountainCarRTAAPlanner(mountaincar, 1000, true_params)
            push!(planners, planner)
            generateHeuristic!(planners)
        end
        agent = MountainCarFiniteModelClassAgent(mountaincar, planners)
        if local_agent
            agent = MountainCarFiniteModelClassLocalAgent(mountaincar, planners)
        end
        n_steps = run(agent, debug = false)
        push!(finite_model_class_steps, n_steps)
        if true_model
            pop!(planners)
        end
        clearResiduals!(planners)
    end
    finite_model_class_steps
end


text:
# Function for running manifold Metropolis adjusted Langevin algorithm (MMALA)
function mmala(model::Model, opts::MmalaOpts)
  mcmc = Array(Float64, opts.mcmc.nPostBurnin, model.nPars)
  z = Array(Float64, opts.mcmc.nPostBurnin, model.nPars)
  invGDerivG = Array(Float64, model.nPars, model.nPars, model.nPars)
  traceInvGDerivG = Array(Float64, model.nPars)
  currentSecondTerm = Array(Float64, model.nPars, model.nPars)
  proposedSecondTerm = Array(Float64, model.nPars, model.nPars)
  proposed, accepted = 0., 0.

  println("Running burn-in iterations...")

  currentPars = model.randPrior()
  currentLogPosterior = model.logPosterior(currentPars)
  currentGradLogPosterior = model.gradLogPosterior(currentPars)
  currentG = model.tensor(currentPars)
  currentInvG = inv(currentG)
  derivG = model.derivTensor(currentPars)
  currentFirstTerm = currentInvG*currentGradLogPosterior
  for i = 1:model.nPars
    invGDerivG[:, :, i] = currentInvG*derivG[:, :, i]
    traceInvGDerivG[i] = trace(invGDerivG[:, :, i])
    currentSecondTerm[:, i] = invGDerivG[:, :, i]*currentInvG[:, i]
  end
  currentThirdTerm = currentInvG*traceInvGDerivG
  driftStep = opts.setDriftStep(1, 0., opts.mcmc.n, opts.mcmc.nBurnin, 0.)

  for i = 1:opts.mcmc.n
    proposed += 1

    proposedPars = copy(currentPars)
    parameterMeans = (proposedPars+(driftStep/2)*currentFirstTerm
      -driftStep*sum(currentSecondTerm, 2)[:]+(driftStep/2)*currentThirdTerm)
    cholCurrentInvG = chol(driftStep*currentInvG)
    proposedPars = parameterMeans+cholCurrentInvG'*randn(model.nPars)
    proposedLogPosterior = model.logPosterior(proposedPars)
    probNewGivenOld = (-sum(log(diag(cholCurrentInvG)))
      -(0.5*(parameterMeans-proposedPars)'
      *(currentG/driftStep)*(parameterMeans-proposedPars))[1])
    proposedG = model.tensor(proposedPars)
    proposedInvG = inv(proposedG)
    proposedGradLogPosterior = model.gradLogPosterior(proposedPars)
    proposedFirstTerm = proposedInvG*proposedGradLogPosterior
    derivG = model.derivTensor(proposedPars)
    for j = 1:model.nPars
      invGDerivG[:, :, j] = proposedInvG*derivG[:, :, j]
      traceInvGDerivG[j] = trace(invGDerivG[:, :, j])
      proposedSecondTerm[:, j] = invGDerivG[:, :, j]*proposedInvG[:, j]
    end
    proposedThirdTerm = proposedInvG*traceInvGDerivG
    parameterMeans = (proposedPars+(driftStep/2)*proposedFirstTerm
      -driftStep*sum(proposedSecondTerm, 2)[:]+(driftStep/2)*proposedThirdTerm)
    probOldGivenNew =
      (-sum(log(diag(chol(driftStep*eye(model.nPars)*proposedInvG))))
      -(0.5*(parameterMeans-currentPars)'
      *(proposedG/driftStep)*(parameterMeans-currentPars))[1])

    ratio =
      proposedLogPosterior+probOldGivenNew-currentLogPosterior-probNewGivenOld

    if ratio > 0 || (ratio > log(rand()))
      accepted += 1

      currentPars = copy(proposedPars)
      currentLogPosterior = copy(proposedLogPosterior)
      currentGradLogPosterior = copy(proposedGradLogPosterior)
      currentG = copy(proposedG)
      currentInvG = copy(proposedInvG)
      currentFirstTerm = copy(proposedFirstTerm)
      currentSecondTerm = copy(proposedSecondTerm)
      currentThirdTerm = copy(proposedThirdTerm)
    end

    if i > opts.mcmc.nBurnin
      mcmc[i-opts.mcmc.nBurnin, :] = currentPars;
      z[i-opts.mcmc.nBurnin, :] = -currentGradLogPosterior/2;
    end

    if mod(i, opts.mcmc.monitorRate) == 0
      acceptanceRatio = accepted/proposed

    driftStep = opts.setDriftStep(i, acceptanceRatio, opts.mcmc.n,
      opts.mcmc.nBurnin, driftStep)

      println("Iteration $i of $(opts.mcmc.n): ", round(100*acceptanceRatio, 2),
        " % acceptance ratio")

      proposed, accepted = 0., 0.
    end

    if i == opts.mcmc.nBurnin
      println("Burn-in completed...\n\nRunning post burn-in MCMC...");
    end
  end

  return mcmc, z
end