text:
# BLIS datatypes for uses within Julia and for calling backends.
#

# BLAS-like type characters.
cblas_typechar = Dict( Float32    => 's',
                       Float64    => 'd',
                       ComplexF32 => 'c',
                       ComplexF64 => 'z' )

BliCompatibleType = Union{keys(cblas_typechar)...}

# Int type should be consistent with system word size.
if Sys.WORD_SIZE == 64
    BliInt  = Int64
    BliUInt = UInt64
elseif Sys.WORD_SIZE == 32
    BliInt  = Int32
    BliUInt = UInt32
end

# Integer parameters.
BliDim  = BliInt
BliInc  = BliInt
BliDoff = BliInt
BliSiz  = BliUInt

# Integer params for object API.
BliObjBits = UInt32
BliAtomic  = ComplexF64

# Enumerates, defined as structures.
#--------------------------------------------------
# Enumerate types generated from list_blis_enums.c
struct BliTrans enum::BliInt end
const BLIS_NO_TRANSPOSE      = BliTrans(0)
const BLIS_TRANSPOSE         = BliTrans(8)
const BLIS_CONJ_NO_TRANSPOSE = BliTrans(16)
const BLIS_CONJ_TRANSPOSE    = BliTrans(24)

struct BliConj enum::BliInt end
const BLIS_NO_CONJUGATE      = BliConj(0)
const BLIS_CONJUGATE         = BliConj(16)

struct BliUpLo enum::BliInt end
const BLIS_ZEROS             = BliUpLo(0)
const BLIS_LOWER             = BliUpLo(192)
const BLIS_UPPER             = BliUpLo(96)
const BLIS_DENSE             = BliUpLo(224)

struct BliSide enum::BliInt end
const BLIS_LEFT              = BliSide(0)
const BLIS_RIGHT             = BliSide(1)

struct BliDiag enum::BliInt end
const BLIS_NONUNIT_DIAG      = BliDiag(0)
const BLIS_UNIT_DIAG         = BliDiag(256)

struct BliInvDiag enum::BliInt end
const BLIS_NO_INVERT_DIAG    = BliInvDiag(0)
const BLIS_INVERT_DIAG       = BliInvDiag(512)

struct BliStruc enum::BliInt end
const BLIS_GENERAL           = BliStruc(0)
const BLIS_HERMITIAN         = BliStruc(134217728)
const BLIS_SYMMETRIC         = BliStruc(268435456)
const BLIS_TRIANGULAR        = BliStruc(402653184)

struct BliNum enum::BliInt end
const BLIS_FLOAT             = BliNum(0)
const BLIS_DOUBLE            = BliNum(2)
const BLIS_SCOMPLEX          = BliNum(1)
const BLIS_DCOMPLEX          = BliNum(3)
const BLIS_INT               = BliNum(4)
const BLIS_CONSTANT          = BliNum(5)
const BLIS_DT_LO             = BliNum(0)
const BLIS_DT_HI             = BliNum(3)

const BLIS_CONJTRANS_BITS    = BliObjBits(24)
const BLIS_DATATYPE_BITS     = BliObjBits(7)
const BLIS_DOMAIN_BIT        = BliObjBits(1)
const BLIS_PRECISION_BIT     = BliObjBits(2)
const BLIS_CONJTRANS_BITS    = BliObjBits(24)
const BLIS_TRANS_BIT         = BliObjBits(8)
const BLIS_CONJ_BIT          = BliObjBits(16)
const BLIS_UPLO_BITS         = BliObjBits(224)
const BLIS_UPPER_BIT         = BliObjBits(32)
const BLIS_DIAG_BIT          = BliObjBits(64)
const BLIS_LOWER_BIT         = BliObjBits(128)
const BLIS_UNIT_DIAG_BIT     = BliObjBits(256)
const BLIS_INVERT_DIAG_BIT   = BliObjBits(512)
const BLIS_TARGET_DT_BITS    = BliObjBits(7168)
const BLIS_TARGET_DOMAIN_BIT = BliObjBits(1024)
const BLIS_TARGET_PREC_BIT   = BliObjBits(2048)
const BLIS_EXEC_DT_BITS      = BliObjBits(57344)
const BLIS_EXEC_DOMAIN_BIT   = BliObjBits(8192)
const BLIS_EXEC_PREC_BIT     = BliObjBits(16384)
const BLIS_STRUC_BITS        = BliObjBits(402653184)

#-------------------------
# Export enumerate types.
export BLIS_NO_TRANSPOSE, BLIS_TRANSPOSE, BLIS_CONJ_NO_TRANSPOSE, BLIS_CONJ_TRANSPOSE
export BLIS_NO_CONJUGATE, BLIS_CONJUGATE
export BLIS_ZEROS, BLIS_LOWER, BLIS_UPPER, BLIS_DENSE
export BLIS_LEFT, BLIS_RIGHT
export BLIS_NONUNIT_DIAG, BLIS_UNIT_DIAG
export BLIS_NO_INVERT_DIAG, BLIS_INVERT_DIAG
export BLIS_GENERAL, BLIS_HERMITIAN, BLIS_SYMMETRIC, BLIS_TRIANGULAR
export BLIS_FLOAT, BLIS_DOUBLE, BLIS_SCOMPLEX, BLIS_DCOMPLEX, BLIS_IN, BLIS_CONSTANT
export BLIS_DT_LO, BLIS_DT_HI

#------------------------------------------
# Lookup table for traditional BLAS chars.
char_to_trans = Dict( 'N' => BLIS_NO_TRANSPOSE,
                      'n' => BLIS_NO_TRANSPOSE,
                      'T' => BLIS_TRANSPOSE,
                      't' => BLIS_TRANSPOSE,
                      'C' => BLIS_CONJ_TRANSPOSE,
                      'c' => BLIS_CONJ_TRANSPOSE )

char_to_conj = Dict( 'N' => BLIS_NO_CONJUGATE,
                     'n' => BLIS_NO_CONJUGATE,
                     'C' => BLIS_CONJUGATE,
                     'c' => BLIS_CONJUGATE )

char_to_side = Dict( 'L' => BLIS_LEFT,
                     'l' => BLIS_LEFT,
                     'R' => BLIS_RIGHT,
                     'r' => BLIS_RIGHT )

char_to_uplo = Dict( 'U' => BLIS_UPPER,
                     'u' => BLIS_UPPER,
                     'L' => BLIS_LOWER,
                     'l' => BLIS_LOWER )

char_to_diag = Dict( 'N' => BLIS_NONUNIT_DIAG,
                     'n' => BLIS_NONUNIT_DIAG,
                     'U' => BLIS_UNIT_DIAG,
                     'u' => BLIS_UNIT_DIAG )

#-----------------------------------------------------------
# Lookup table from Julia datatype to BLIS type indicators.
ctype_to_bli_num = Dict( Float32    => BLIS_FLOAT,
                         Float64    => BLIS_DOUBLE,
                         ComplexF32 => BLIS_SCOMPLEX,
                         ComplexF64 => BLIS_DCOMPLEX,
                         BliInt     => BLIS_INT )


text:
const Bst_ulong = Culonglong

"Calls an xgboost API function and correctly reports errors."
macro xgboost(f, params...)
    return quote
        err = ccall(($f, libxgboost), Int64,
                    ($((esc(i.args[end]) for i in params)...),),
                    $((esc(i.args[end-1]) for i in params)...))
        if err != 0
            err_msg = unsafe_string(ccall((:XGBGetLastError, libxgboost), Cstring, ()))
            error("Call to XGBoost C function ", string($(esc(f))), " failed: ", err_msg)
        end
    end
end

function XGDMatrixCreateFromFile(fname::String, silent::Int32)
    out = Ref{Ptr{Nothing}}()
    @xgboost(:XGDMatrixCreateFromFile,
             fname => Cstring,
             silent => Cint,
             out => Ref{Ptr{Nothing}})
    return out[]
end

@deprecate XGDMatrixCreateFromCSC(data) XGDMatrixCreateFromCSCEx(data) false

function XGDMatrixCreateFromCSCEx(data::SparseMatrixCSC)
    out = Ref{Ptr{Nothing}}()
    @xgboost(:XGDMatrixCreateFromCSCEx,
             convert(Vector{UInt64}, data.colptr .- 1) => Ptr{Csize_t},
             convert(Vector{UInt32}, data.rowval .- 1) => Ptr{Cuint},
             convert(Vector{Float32}, data.nzval) => Ptr{Cfloat},
             convert(UInt64, size(data.colptr)[1]) => Csize_t,
             convert(UInt64, nnz(data)) => Csize_t,
             convert(UInt64, size(data)[1]) => Csize_t, #num_rows
             out => Ref{Ptr{Nothing}})
    return out[]
end

function XGDMatrixCreateFromCSCT(data::SparseMatrixCSC)
    handle = Ref{Ptr{Nothing}}()
    @xgboost(:XGDMatrixCreateFromCSREx,
             convert(Vector{UInt64}, data.colptr .- 1) => Ptr{Csize_t},
             convert(Vector{UInt32}, data.rowval .- 1) => Ptr{Cuint},
             convert(Vector{Float32}, data.nzval) => Ptr{Cfloat},
             convert(UInt64, size(data.colptr)[1]) => Csize_t,
             convert(UInt64, nnz(data)) => Csize_t,
             convert(UInt64, size(data)[2]) => Csize_t, #num_cols
             handle => Ref{Ptr{Nothing}})
    return handle[]
end

function XGDMatrixCreateFromMat(data::Matrix{Float32}, missing::Float32)
    XGDMatrixCreateFromMatT(Matrix(transpose(data)), missing)
end

function XGDMatrixCreateFromMatT(data::Matrix{Float32}, missing::Float32)
    ncol, nrow = size(data)
    handle = Ref{Ptr{Nothing}}()
    @xgboost(:XGDMatrixCreateFromMat,
             data => Ptr{Cfloat},
             nrow => Bst_ulong,
             ncol => Bst_ulong,
             missing => Cfloat,
             handle => Ref{Ptr{Nothing}})
    return handle[]
end

function XGDMatrixSliceDMatrix(handle::Ptr{Nothing}, idxset::Vector{Int32}, len::UInt64)
    ret = Ref{Ptr{Nothing}}()
    @xgboost(:XGDMatrixSliceDMatrix,
             handle => Ptr{Nothing},
             idxset => Ptr{Cint},
             len => Bst_ulong,
             ret => Ref{Ptr{Nothing}})
    return ret[]
end

function XGDMatrixFree(handle::Ptr{Nothing})
    @xgboost(:XGDMatrixFree,
             handle => Ptr{Nothing})
end

function XGDMatrixSaveBinary(handle::Ptr{Nothing}, fname::String, silent::Int32)
    @xgboost(:XGDMatrixSaveBinary,
             handle => Ptr{Nothing},
             fname => Cstring,
             silent => Cint)
end

function XGDMatrixSetFloatInfo(handle::Ptr{Nothing}, field::String, array::Vector{Float32},
                               len::UInt64)
    @xgboost(:XGDMatrixSetFloatInfo,
             handle => Ptr{Nothing},
             field => Cstring,
             array => Ptr{Cfloat},
             len => Bst_ulong)
end

function XGDMatrixSetUIntInfo(handle::Ptr{Nothing}, field::String, array::Vector{UInt32},
                              len::UInt64)
    @xgboost(:XGDMatrixSetUIntInfo,
             handle => Ptr{Nothing},
             field => Cstring,
             array => Ptr{Cuint},
             len => Bst_ulong)
end

function XGDMatrixSetGroup(handle::Ptr{Nothing}, array::Vector{UInt32}, len::UInt64)
    @xgboost(:XGDMatrixSetGroup,
             handle => Ptr{Nothing},
             array => Ptr{Cuint},
             len => Bst_ulong)
end

function XGDMatrixGetFloatInfo(handle::Ptr{Nothing}, field::String, out_len::Vector{Bst_ulong})
    out_dptr = Ref{Ptr{Cfloat}}()
    @xgboost(:XGDMatrixGetFloatInfo,
             handle => Ptr{Nothing},
             field => Cstring,
             out_len => Ptr{Bst_ulong},
             out_dptr =>  Ref{Ptr{Cfloat}})
    return out_dptr[]
end

function JLGetFloatInfo(handle::Ptr{Nothing}, field::String)
    len = Bst_ulong[1]
    ptr = XGDMatrixGetFloatInfo(handle, field, len)
    return unsafe_wrap(Array, ptr, len[1])
end

function JLGetUintInfo(handle::Ptr{Nothing}, field::String)
    len = Bst_ulong[1]
    ptr = XGDMatrixGetUIntInfo(handle, field, len)
    return unsafe_wrap(Array, ptr, len[1])
end

function XGDMatrixGetUIntInfo(handle::Ptr{Nothing}, field::String, out_len::Vector{Bst_ulong})
    out_dptr = Ref{Ptr{Cuint}}()
    @xgboost(:XGDMatrixGetUIntInfo,
             handle => Ptr{Nothing},
             field => Cstring,
             out_len => Ptr{Bst_ulong},
             out_dptr => Ref{Ptr{Cuint}})
    return out_dptr[]
end

function XGDMatrixNumRow(handle::Ptr{Nothing})
    out = Ref{Bst_ulong}()
    @xgboost(:XGDMatrixNumRow,
             handle => Ptr{Nothing},
             out => Ref{Bst_ulong})
    return out[]
end

function XGBoosterCreate(cachelist::Vector{Ptr{Nothing}}, len::Int64)
    out = Ref{Ptr{Nothing}}()
    @xgboost(:XGBoosterCreate,
             cachelist => Ptr{Ptr{Nothing}},
             len => Bst_ulong,
             out => Ref{Ptr{Nothing}})
    return out[]
end

function XGBoosterFree(handle::Ptr{Nothing})
    @xgboost(:XGBoosterFree,
             handle => Ptr{Nothing})
end

function XGBoosterSetParam(handle::Ptr{Nothing}, name::String, value::String)
    @xgboost(:XGBoosterSetParam,
             handle => Ptr{Nothing},
             name => Cstring,
             value => Cstring)
end

function XGBoosterUpdateOneIter(handle::Ptr{Nothing}, iter::Int32, dtrain::Ptr{Nothing})
    @xgboost(:XGBoosterUpdateOneIter,
             handle => Ptr{Nothing},
             iter => Cint,
             dtrain => Ptr{Nothing})
end

function XGBoosterBoostOneIter(handle::Ptr{Nothing}, dtrain::Ptr{Nothing}, grad::Vector{Float32},
                               hess::Vector{Float32}, len::UInt64)
    @xgboost(:XGBoosterBoostOneIter,
             handle => Ptr{Nothing},
             dtrain => Ptr{Nothing},
             grad => Ptr{Cfloat},
             hess => Ptr{Cfloat},
             len => Bst_ulong)
end

function XGBoosterEvalOneIter(handle::Ptr{Nothing}, iter::Int32, dmats::Vector{Ptr{Nothing}},
                              evnames::Vector{String}, len::UInt64)
    out_result = Ref{Cstring}()
    @xgboost(:XGBoosterEvalOneIter,
             handle => Ptr{Nothing},
             iter => Cint,
             dmats => Ptr{Ptr{Nothing}},
             evnames => Ptr{Cstring},
             len => Bst_ulong,
             out_result => Ref{Cstring})
    return unsafe_string(out_result[])
end

function XGBoosterPredict(handle::Ptr{Nothing}, dmat::Ptr{Nothing}, option_mask::Int32,
                          ntree_limit::UInt32, training::Int32, out_len::Vector{UInt64})
    out_result = Ref{Ptr{Float32}}()
    @xgboost(:XGBoosterPredict,
             handle => Ptr{Nothing},
             dmat => Ptr{Nothing},
             option_mask => Cint,
             ntree_limit => Cuint,
             training => Cint,
             out_len => Ptr{Bst_ulong},
             out_result => Ref{Ptr{Cfloat}})
    return out_result[]
end

function XGBoosterLoadModel(handle::Ptr{Nothing}, fname::String)
    @xgboost(:XGBoosterLoadModel,
             handle => Ptr{Nothing},
             fname => Cstring)
end

function XGBoosterSaveModel(handle::Ptr{Nothing}, fname::String)
    @xgboost(:XGBoosterSaveModel,
             handle => Ptr{Nothing},
             fname => Cstring)
end

function XGBoosterDumpModel(handle::Ptr{Nothing}, fmap::String, with_stats::Int64, dump_format::String="text")
    out_dump_array = Ref{Ptr{Cstring}}()
    out_len = Ref{Bst_ulong}(0)
    @xgboost(:XGBoosterDumpModelEx,
             handle => Ptr{Nothing},
             fmap => Cstring,
             with_stats => Cint,
             dump_format => Cstring,
             out_len => Ref{Bst_ulong},
             out_dump_array => Ref{Ptr{Cstring}})
    return unsafe_wrap(Array, out_dump_array[], out_len[])
end

ext: jl


text:
using DrakeVisualizer
using GeometryTypes
using CoordinateTransformations
using Meshing
using Base.Test
using IJulia
import Iterators: product

proc = DrakeVisualizer.new_window()

try
    @testset "robot_load" begin
        f = x -> norm(x)^2
        bounds = HyperRectangle(Vec(0.,0,0), Vec(1.,1,1))
        geom = contour_mesh(f, minimum(bounds), maximum(bounds))
        robot = convert(Robot, geom)
        vis = Visualizer(robot, 1)
    end

    @testset "link_list_load" begin
        links = Link[]
        levels = [0.5; 1]
        for i = 1:2
            f = x -> norm(x)^2
            geom = contour_mesh(f, Vec(0.,0,0), Vec(1.,1,1), levels[i])
            push!(links, Link(geom))
        end
        vis = Visualizer(links, 1)
    end

    @testset "geom_load" begin
        f = x -> norm(x)^2
        iso_level = 0.5
        geom = GeometryData(contour_mesh(f, [0.,0,0], [1.,1,1], iso_level))
        vis = Visualizer(geom)
    end

    @testset "robot_draw" begin
        links = Link[]
        link_lengths = [1.; 2; 3]
        for (i, l) in enumerate(link_lengths)
            geometry = HyperRectangle(Vec(0., -0.1, -0.1), Vec(l, 0.2, 0.2))
            geometry_data = GeometryData(geometry)
            push!(links, Link([geometry_data], "link$(i)"))
        end

        function link_origins(joint_angles)
            transforms = Array{Transformation}(length(link_lengths))
            transforms[1] = LinearMap(AngleAxis(joint_angles[1], 0, 0, 1.0))
            for i = 2:length(link_lengths)
                T = compose(Translation(link_lengths[i-1], 0, 0),
                            LinearMap(AngleAxis(joint_angles[i], 0, 0, 1.0))
                            )
                transforms[i] = compose(transforms[i-1], T)
            end
            transforms
        end

        robot = Robot(links)
        model = Visualizer(robot)

        for x in product([linspace(-pi, pi, 11) for i in 1:length(link_lengths)]...)
            origins = link_origins(reverse(x))
            draw(model, origins)
        end
    end

    @testset "ellipsoid" begin
        ellipsoid = DrakeVisualizer.HyperEllipsoid(Point(1.,0,0.1), Vec(0.3, 0.2, 0.1))
        Visualizer(ellipsoid)
    end

    @testset "cylinder" begin
        cylinder = DrakeVisualizer.HyperCylinder{3, Float64}(1.0, 0.5)
        Visualizer(cylinder)
    end

    @testset "multiple_robots" begin
        robot1 = convert(Robot, DrakeVisualizer.HyperCylinder{3, Float64}(1.0, 0.5))
        robot2 = convert(Robot, DrakeVisualizer.HyperRectangle(Vec(0., -0.1, -0.1), Vec(1, 0.2, 0.2)))
        vis1, vis2 = Visualizer([robot1, robot2])
        for vis in [vis1, vis2]
            draw(vis, [IdentityTransformation()])
        end
    end

    @testset "demo_notebook" begin
        if VERSION < v"0.6-dev"
            jupyter = IJulia.jupyter
            demo_file = "../demo.ipynb"
            run(`$jupyter nbconvert --to notebook --execute $(demo_file) --output $(demo_file)`)
        end
    end
finally
    kill(proc)
end


text:
module PandorasBoxLamp

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/lamp" Lamp(x::Integer, y::Integer, inverted::Bool=false, flag::String="", baseColor::String="White", lightColor::String="White", lightMode::String="Smooth", lightStartRadius::Int=48, lightEndRadius::Int=64)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

const placements = Ahorn.PlacementDict(
    "Lamp (Pandora's Box)" => Ahorn.EntityPlacement(
        Lamp,
        "point"
    )
)

const lightModes = String[
    "Smooth",
    "Instant"
]

Ahorn.editingOptions(entity::Lamp) = Dict{String, Any}(
    "baseColor" => colors,
    "lightColor" => colors,
    "lightMode" => lightModes
)

function getColor(color)
    if haskey(Ahorn.XNAColors.colors, color)
        return Ahorn.XNAColors.colors[color]

    else
        try
            return ((Ahorn.argb32ToRGBATuple(parse(Int, replace(color, "#" => ""), base=16))[1:3] ./ 255)..., 1.0)

        catch

        end
    end

    return (1.0, 1.0, 1.0, 1.0)
end

gemTexture = "objects/pandorasBox/lamp/idle0"
baseTexture = "objects/pandorasBox/lamp/base"

function Ahorn.selection(entity::Lamp)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(gemTexture, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Lamp, room::Maple.Room)
    colorRaw = get(entity, "baseColor", "White")
    color = getColor(colorRaw)

    Ahorn.drawSprite(ctx, baseTexture, 0, 0, tint=color)
    Ahorn.drawSprite(ctx, gemTexture, 0, 0)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
module AnimTest

using SDLGamesLib
using SDLGamesLib.Animations: Animation, loopanim_callback, render
using SDLGamesLib: UIPixelPos, Camera

using Test

# Mock out `render` function for tests
struct MockFrame
    x
end
function SDLGamesLib.render(f::MockFrame, pos, camera, renderer; kwargs...)
    f.x
end

@testset "Animation render & update!" begin
    @testset "Simple animation" begin
        every_sec = Animation([MockFrame(x) for x in 1:3], [1,1,1])
        @test 1 == render(every_sec, UIPixelPos(0,0), Camera(), nothing; size=1)

        # After 0.5 secs, still on first frame
        update!(every_sec, 0.5)
        @test 1 == render(every_sec, UIPixelPos(0,0), Camera(), nothing; size=1)

        # After 0.6 secs, into the second frame
        update!(every_sec, 0.6)
        @test 2 == render(every_sec, UIPixelPos(0,0), Camera(), nothing; size=1)

        # At *exactly* 3 seconds, onto the third frame  (from the `>=` in update!)
        update!(every_sec, 0.9)
        @test 3 == render(every_sec, UIPixelPos(0,0), Camera(), nothing; size=1)

        # After all subsequent updates, it's still emitting the last frame.
        update!(every_sec, 2)
        @test 3 == render(every_sec, UIPixelPos(0,0), Camera(), nothing; size=1)
    end

    @testset "skip frames" begin
        a = Animation([MockFrame(x) for x in 1:10], ones(10))
        update!(a, 1)  # start with an update
        @test 2 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        update!(a, 5)  # skip frames
        @test 7 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        update!(a, 5)  # skip past end (10)
        @test 10 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
    end

    @testset "variable timings" begin
        a = Animation([MockFrame(x) for x in 1:3], [1, 10, 1])
        @test 1 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        update!(a, 1)
        @test 2 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        update!(a, 1)
        @test 2 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        update!(a, 9)
        @test 3 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
    end
end

@testset "Animation donecallback" begin
    let _rem_dt = nothing
        _test_callback(a, rem_dt) = (_rem_dt = rem_dt)

        a = Animation([MockFrame(x) for x in 1:2], [1,1], _test_callback)
        # During last frame, callback hasn't triggered
        update!(a, 1)
        @test 2 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        @test _rem_dt == nothing

        # After finishing last frame, callback triggered.
        update!(a, 1)  # Jumping _exactly_, so no remainder
        @test 2 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        @test _rem_dt == 0

        # Now, further updates don't trigger the callback
        _test_callback(nothing, nothing)  # Reset
        update!(a, 1)  # Ignored
        @test 2 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        @test _rem_dt == nothing
    end

    let _rem_dt = nothing
        _test_callback(a, rem_dt) = (_rem_dt = rem_dt)

        # Jumping _over_ the last frame sets a remainder
        a = Animation([MockFrame(x) for x in 1:2], [1,1], _test_callback)
        update!(a, 2.5)
        @test 2 == render(a, UIPixelPos(0,0), Camera(), nothing; size=1)
        @test _rem_dt == 0.5
    end
end

@testset "looping" begin
    looping = Animation([MockFrame(x) for x in 1:3], [1,1,1], loopanim_callback)
    # Jump to last frame
    update!(looping, 2)
    @test 3 == render(looping, UIPixelPos(0,0), Camera(), nothing; size=1)

    # After finishing the last frame, it loops back to the beginning
    update!(looping, 1)
    @test 1 == render(looping, UIPixelPos(0,0), Camera(), nothing; size=1)

    # updates longer than the entire animation wrap
    update!(looping, 3)
    @test 1 == render(looping, UIPixelPos(0,0), Camera(), nothing; size=1)
    update!(looping, 5)
    @test 3 == render(looping, UIPixelPos(0,0), Camera(), nothing; size=1)
end

end


text:
# Practical Machine learning
# Deep learning - Sparse autoencoder example
# Chapter 11

using SparseAutoencoder
using MAT

vars = matread("./data/mnist-images.mat")
data = vars["images"]

visiblesize = 28*28
hiddensize = 196
sparsityparameter = 0.1
lambda = 3e-3
beta = 3.0
patches = data[:,1:10000]

minf,W1,W2,b1,b2 = autoencode(patches,hiddensize,visiblesize,lambda=lambda,beta=beta,rho=sparsityparameter)

using HDF5, JLD
@save "./digits-results.jld"


text:
using PyCall

check() = print("Docker is Up and Running and Python work")
check()

function run()
    while true
    end
end

run()


text:
# Named Tuples
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_named_tuples
x = (a = 1, b = 2)
x.a
x.b

# Functions: Compact form
f(x, y) = x + y
f(3,2)

# Anonymous Functions
x -> x^2 + 2x + 1

# Keyword Arguments
function myplot(x, y; style = "solid", width = 10)
    println("style: $style")
    println("width: $width")

    nothing
end

myplot(10, 3, style = "dash", width = 13)

# Closures
# Function uses variable defined outside of calling scope
foo(x) = () -> x^2
bar = foo(2) # closure
bar() # 4

# -------
# Blocks
# -------
begin
    bar = foo(3)
    bar()
end

# let Blocks
x, y, z = 1, 2, 3
let x = 1, z
    @show x y z
end
@show x y z

# do Blocks
# Auto close files

# Style 1
open("testfile.txt", "w") do io
    println(io, "Hello World!")
end

# Style 2
open(io->read(io, String), "testfile.txt")

# -----------------
# Control Flow: Ternary Operator
# -----------------

# -------------------
# Control Flow: Coroutines or Task
# -------------------
function generator(ch::Channel)
    while true
        put!(ch, rand(Float64))
    end
end

randgen = Channel(generator)
take!(randgen)
take!(randgen)

# ------
# Types
# ------

# Parametric Types
struct Point{T<:Real}
    x::T
    y::T
end

pt = Point(0.0, 0.0)
pt = Point(1,2)

# --------
# Methods
# --------
# Parametric Methods
# Short form: isintpoint(pt::Point{T}) where {T} = (T === Int64)
function isintpoint(pt::Point{T}) where {T}
    (T === Int64)
end

isintpoint(pt)


text:
#=
step distance search helpers for s,z
=#

mutable struct StepSearcher{T <: Real}
    min_prox::T
    prox_bound::T
    use_max_prox::Bool
    alpha_sched::Vector{T}

    szk::Vector{T}
    nup1::T
    cone_times::Vector{Float64}
    cone_order::Vector{Int}
    prev_sched::Int
    prox::T

    function StepSearcher{T}(
        model::Models.Model{T};
        min_prox::T = T(0.01),
        prox_bound::T = T(0.99),
        use_max_prox::Bool = true,
        alpha_sched::Vector{T} = default_alpha_sched(T),
        ) where {T <: Real}
        cones = model.cones
        searcher = new{T}()
        searcher.min_prox = min_prox
        searcher.prox_bound = prox_bound
        searcher.use_max_prox = use_max_prox
        searcher.alpha_sched = alpha_sched
        searcher.szk = zeros(T, length(cones))
        searcher.nup1 = T(model.nu + 1)
        searcher.cone_times = zeros(length(cones))
        searcher.cone_order = collect(1:length(cones))
        searcher.prev_sched = 0
        searcher.prox = zero(T)
        return searcher
    end
end

default_alpha_sched(T::Type{<:Real}) = T[
    0.9999, 0.999, 0.99, 0.97, 0.95, 0.9, 0.85, 0.8, 0.7, 0.6, 0.5,
    0.3, 0.1, 0.05, 0.01, 0.005, 0.001, 0.0005]

# backwards search on alphas in alpha schedule
function search_alpha(
    point::Point{T},
    model::Models.Model{T},
    stepper::Stepper{T};
    sched::Int = start_sched(stepper, stepper.searcher)
    ) where {T <: Real}
    searcher = stepper.searcher

    while sched <= length(searcher.alpha_sched)
        alpha = searcher.alpha_sched[sched]
        # update ztsk only in stepper.temp
        update_stepper_points(alpha, point, stepper, true)

        # NOTE updates cone points and grad
        if check_cone_points(model, stepper)
            searcher.prev_sched = sched
            return alpha
        end
        sched += 1
    end

    searcher.prev_sched = sched
    return zero(T)
end

# fallback starts at first alpha in schedule
start_sched(stepper::Stepper, searcher::StepSearcher) = 1

function check_cone_points(
    model::Models.Model{T},
    stepper::Stepper{T};
    ) where {T <: Real}
    searcher = stepper.searcher
    cand = stepper.temp
    szk = searcher.szk
    cones = model.cones
    min_prox = searcher.min_prox
    use_max_prox = searcher.use_max_prox
    proxsqr_bound = abs2(searcher.prox_bound)

    taukap = cand.tau[] * cand.kap[]
    (min(cand.tau[], cand.kap[], taukap) < eps(T)) && return false

    for k in eachindex(cones)
        szk[k] = dot(cand.primal_views[k], cand.dual_views[k])
        (szk[k] < eps(T)) && return false
    end
    mu = (sum(szk) + taukap) / searcher.nup1
    (mu < eps(T)) && return false

    taukap_rel = taukap / mu
    (taukap_rel < min_prox) && return false
    taukap_proxsqr = abs2(taukap_rel - 1)
    (taukap_proxsqr > proxsqr_bound) && return false

    for k in eachindex(cones)
        nu_k = Cones.get_nu(cones[k])
        sz_rel_k = szk[k] / (mu * nu_k)
        if (sz_rel_k < min_prox) || (nu_k * abs2(sz_rel_k - 1) > proxsqr_bound)
            return false
        end
    end

    # order the cones by how long it takes to check proximity condition and
    # iterate in that order, to improve efficiency
    cone_order = searcher.cone_order
    sortperm!(cone_order, searcher.cone_times, initialized = true) # stochastic

    irtmu = inv(sqrt(mu))
    agg_proxsqr = taukap_proxsqr
    aggfun = (use_max_prox ? max : +)

    for k in cone_order
        cone_k = cones[k]
        start_time = time()
        Cones.load_point(cone_k, cand.primal_views[k], irtmu)
        Cones.load_dual_point(cone_k, cand.dual_views[k])
        Cones.reset_data(cone_k)

        in_prox_k = false
        if Cones.is_feas(cone_k) && Cones.is_dual_feas(cone_k) &&
            Cones.check_numerics(cone_k)
            proxsqr_k = Cones.get_proxsqr(cone_k, irtmu, use_max_prox)
            agg_proxsqr = aggfun(agg_proxsqr, proxsqr_k)
            in_prox_k = (agg_proxsqr < proxsqr_bound)
        end
        searcher.cone_times[k] = time() - start_time
        in_prox_k || return false
    end

    searcher.prox = sqrt(agg_proxsqr)
    return true
end


text:
function predictionGSBR(x,degree,maxiter,burnin,seed=1,T=0,filename = "/results")
  Random.seed!(seed)
  savelocation = string(pwd(), filename, "/gsbr_seed$seed")
  mkpath(savelocation)

  n = length(x)
  if T > 0
    append!(x, zeros(T))
    xpred = fill(Float64[],T)
    for t in 1:1:T
      xpred[t] = zeros(maxiter-burnin)
    end
    n_ = n + T
  else
    n_ = n
  end

  degree += 1
  theta = zeros(degree) # polynomial coeeficients vector
  sampledTheta = zeros(maxiter-burnin,degree) # matrix to store sample thetas
  thetaprev = copy(theta)
  epsj = 20.
  x0 = 0.0
  sampledX0 = zeros(maxiter-burnin)
  zp = zeros(maxiter-burnin)
  clusters = zeros(Int64,maxiter)
  L = zeros(maxiter)
  # MCMC OPTIONS
  ### specify geometric probability and the hyperparameters for the transformed gamma prior
  lambda = 0.5
  alambda,blambda= 1.,1.

  ### shape and rate for the gamma prior on precisions
  at, bt = 3,0.001

  # GSB meausre latent variables
  d = ones(Int64,n_)  # clustering variables
  N = copy(d)         # geometric slice variables

  ### begin MCMC
  for its in 1:1:maxiter

    M = maximum(N)
    w = zeros(M)
    tau = zeros(M)
    for j in 1:1:M
      w[j] = lambda * (1-lambda)^(j-1)
    end

    # sample precisions
    for j in 1:1:M
      term,counts = 0.,0.
      if d[1].==j
        counts += 1
        term += (x[1] - polyMap(theta,x0))^2
      end
      for i in 2:1:n_
        if d[i].==j
          counts += 1
          term += (x[i] - polyMap(theta,x[i-1]))^2
        end
      end
      tau[j] = rand(Gamma(at + 0.5counts,(bt + 0.5term)^(-1)))
    end

    # sample clustering variables
    nc0 = 0.
    for j in 1:1:N[1]
      nc0 += tau[j]^(0.5) * exp(-0.5*tau[j]*(x[1]-polyMap(theta,x0))^2)
    end
    rd0 = rand()
    prob0 = 0.
    for j in 1:1:N[1]
      prob0 += tau[j]^(0.5) * exp(-0.5*tau[j]*(x[1]-polyMap(theta,x0))^2) / nc0
      if rd0.<prob0
        d[1] = j
        break
      end
    end

    for i in 2:1:n_
      nc = 0.
      for j in 1:1:N[i]
        nc += tau[j]^(0.5) * exp(-0.5*tau[j]*(x[i]-polyMap(theta,x[i-1]))^2)
      end
      prob = 0.
      rd = rand()
      for j in 1:1:N[i]
        prob += tau[j]^(0.5) * exp(-0.5*tau[j]*(x[i]-polyMap(theta,x[i-1]))^2) / nc
        if rd.<prob
          d[i] = j
          break
        end
      end
    end

    clusters[its] = length(unique(d))
    ### update the geometric slice variables
    for i in 1:1:n_
      N[i] = tgeornd(lambda, d[i])
    end

    ### update probability lambda
    lambda = rand(Beta(alambda + 2n_, blambda + sum(N) - n_)) # Beta conjugate
    #lambda = synchronized(alambda, blambda, lambda, sum(N), n_)  # transformed gamma prior
    L[its] = lambda

    if its > burnin
      # sample from noise predictive
      W = cumsum(w)
      rp = rand()
      if rp.>W[end]
        zp[its-burnin] = rand(Normal(0, sqrt(1/rand(Gamma(at,1/bt)))))
      else
        for i in 1:1:length(W)
          if rp.<W[i]
            zp[its-burnin] = rand(Normal(0, sqrt(1/tau[i])))
            break
          end
        end
      end

      # sample the vector with the coefficients
      for j in 1:1:length(theta)
        thetaj = copy(theta)
        thetaj[j] = 0.

        tauj = tau[d[1]]*x0^(2(j-1))
        meanj = tau[d[1]]*x0^(j-1) * (x[1]-polyMap(thetaj,x0))
        for i in 2:1:n_
          tauj += tau[d[i]]*x[i-1]^(2(j-1))
          meanj += tau[d[i]]*x[i-1]^(j-1) * (x[i] - polyMap(thetaj,x[i-1]))
        end
        meanj /= tauj
        vj = -(2/tauj)*log(rand()) + (thetaprev[j]-meanj)^2
        theta[j] = rand(Uniform(max(-sqrt(vj)+meanj,-epsj), min(sqrt(vj)+meanj,epsj)))
      end
      thetaprev = copy(theta)
      sampledTheta[its-burnin,:] = copy(thetaprev)
      # sample x0
      aux = -(2.0/tau[d[1]])*log(rand()) + (x[1]-polyMap(theta,x0))^2
      poly_right = copy(theta)
      poly_right[1] = poly_right[1] - (x[1] - sqrt(aux))
      poly_right = poly_right./poly_right[end]
      roots_right = Polynomials.roots(Polynomial(poly_right))
      poly_left = copy(theta)
      poly_left[1] = poly_left[1] - (x[1] + sqrt(aux))
      poly_left = poly_left./poly_left[end]
      roots_left = Polynomials.roots(Polynomial(poly_left))
      allroots = [roots_left roots_right]
      idx = abs.(imag(allroots)).<1e-08
      realroots = sort(real(allroots[idx]))
      intervals = range_intersection(realroots,[-20. 20.])
      x0 = unifmixrnd(intervals)
      sampledX0[its-burnin] = x0

      # predict future unobserved observations directly
      # good results only for the noise process f₂
      if T > 0
        for j in 1:1:T-1
          ax1 = -(2/tau[d[n+j]])*log(rand()) + (x[n+j]-polyMap(theta,x[n+j-1]))^2
          ax2 = -(2/tau[d[n+j+1]])*log(rand()) + (x[n+j+1]-polyMap(theta,x[n+j]))^2

          xnm = polyMap(theta,x[n+j-1]) - sqrt(ax1)
          xnp = polyMap(theta,x[n+j-1]) + sqrt(ax1)

          poly_right_x = copy(theta)
          poly_right_x[1] = poly_right_x[1] - (x[n+j+1]-sqrt(ax2))
          poly_right_x = poly_right_x./poly_right_x[end]
          roots_right_x = Polynomials.roots(Polynomial(poly_right_x))

          poly_left_x = copy(theta)
          poly_left_x[1] = poly_left_x[1] - (x[n+j+1]+sqrt(ax2))
          poly_left_x = poly_left_x./poly_left_x[end]
          roots_left_x = Polynomials.roots(Polynomial(poly_left_x))

          allrootsx = [roots_left_x roots_right_x]
          idxx = abs.(imag(allrootsx)).<1e-06
          realrootsx = sort(real(allrootsx[idxx]))
          intervalsx = range_intersection(realrootsx,[xnm xnp])
          x[n+j] = unifmixrnd(intervalsx)
          xpred[j][its-burnin] = x[n+j]
        end
        new_mu = polyMap(theta,x[n+T-1])
        new_var = 1/tau[d[n+T]]
        x[n+T] = rand(Normal(new_mu, sqrt(new_var)))
        xpred[T][its-burnin] = x[n+T]
      end

    end

    if mod(its,10000)==0
      println("MCMC iterations: $its")
    end
  # end gibbs
  end

  ### write outputs to data files
  writedlm(string(savelocation, "/thetas.txt"), sampledTheta)
  writedlm(string(savelocation, "/X0.txt"), sampledX0)
  writedlm(string(savelocation, "/clusters.txt"), clusters)
  writedlm(string(savelocation, "/Prob.txt"), L)
  writedlm(string(savelocation, "/noise.txt"), zp)
  writedlm(string(savelocation, "/data.txt"), x)

  if T > 0
    for i in 1:1size(xpred)[1]
      writedlm(string(savelocation, "/xpred$i.txt"), xpred[i])
    end
    return sampledTheta, zp, sampledX0, xpred, clusters
  else
    return sampledTheta, zp, sampledX0, clusters
  end
end

max_stars_repo_name: cmerkatas/GSBR

max_issues_repo_name: cmerkatas/GSBR

max_forks_repo_name: cmerkatas/GSBR


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -8.232182974914569 <= q <= 6.070970235703473)

                     add_NL_constraint(m, :(softplus(-0.17705543456567296 + 0.5634748900178774*softplus(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])) + 0.6106029608527352*softplus(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])) + 0.5384388086098779*softplus(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))) + softplus(0.08053767290922664 + 0.911477397877237*softplus(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])) + 0.7134526118957916*softplus(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])) + -0.34097520379571256*softplus(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))) + softplus(0.5992621446321471 + 0.15304126986204158*softplus(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])) + -0.37134199071135*softplus(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])) + 0.8747239793986465*softplus(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
    clip(layer::T, p1::Point, p2::Point) where {T <: SimpleSDMLayer}

Return a raster by defining a bounding box through two points. The order of the
points (in terms of bottom/top/left/right) is not really important, as the
correct coordinates will be extracted.
"""
function clip(layer::T, p1::Point, p2::Point) where {T <: SimpleSDMLayer}
    latextr = extrema([p1[2], p2[2]])
    lonextr = extrema([p1[1], p2[1]])
    pmin = _point_to_cartesian(layer, Point(minimum(lonextr), minimum(latextr)); side=:bottomleft)
    pmax = _point_to_cartesian(layer, Point(maximum(lonextr), maximum(latextr)); side=:topright)
    R = T <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
    return R(
        layer.grid[pmin:pmax],
        longitudes(layer)[last(pmin.I)]-stride(layer, 1),
        longitudes(layer)[last(pmax.I)]+stride(layer, 1),
        latitudes(layer)[first(pmin.I)]-stride(layer, 2),
        latitudes(layer)[first(pmax.I)]+stride(layer, 2)
    )
end

"""
    clip(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer}

Clips a raster by giving the (optional) limites `left`, `right`, `bottom`, and `top`.
"""
function clip(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer}
    p1 = Point(
        isnothing(left) ? layer.left : left,
        isnothing(bottom) ? layer.bottom : bottom
    )
    p2 = Point(
        isnothing(right) ? layer.right : right,
        isnothing(top) ? layer.top : top
    )
    return clip(layer, p1, p2)
end

"""
    clip(origin::T1, destination::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}

Clips a layer by another layer, *i.e.* subsets the first layer so that it has
the dimensions of the second layer. This operation applies a very small
displacement to the limits (`5eps()`) to ensure that if the coordinate to cut at
falls exactly on a cell boundary, the correct cell will be used in the return
layer.
"""
function clip(origin::T1, destination::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
    _d = 5eps()
    err = false
    destination.right > origin.right && (err = true)
    destination.left < origin.left && (err = true)
    destination.bottom < origin.bottom && (err = true)
    destination.top > origin.top && (err = true)
    err && throw(ArgumentError("The two layers are not compatible"))
    return clip(origin, Point(destination.left+_d, destination.bottom+_d), Point(destination.right-_d, destination.top-_d))
end


text:
mutable struct InventorySystem <: RLBase.AbstractEnv
    t::Int
    T::Int
    bom::Vector{AbstractItem}
    constraints::Vector{AbstractConstraint}
    reward::Float64
end

function InventorySystem(T, bom::Vector{<:AbstractItem}, constraints = AbstractConstraint[])
    constraint_names = unique(map(c -> getfield(c, :name), constraints))
    if length(constraint_names) < length(constraints) || "" in constraint_names
        for (i,c) in enumerate(constraints)
            if c.name == ""
                c.name = "constraint_$i"
            else
                c.name *= "$i"
            end
        end
    end
    item_names = unique(map(it -> getfield(it, :name), bom))
    if length(item_names) < length(bom) || "" in item_names
        for (i,it) in enumerate(bom)
            if it.name == ""
                it.name = "item_$i"
            else
                it.name *= "$i"
            end
        end
    end
    maxT = T == Inf ? typemax(Int) : Int(T)
    bom_to = topological_order(bom)
    InventorySystem(1, maxT, bom_to, constraints, 0.)
end

RLBase.state(is::InventorySystem) = reduce(vcat, map(state,is.bom))
state_size(is::InventorySystem) = sum(map(state_size,is.bom))
action_size(is::InventorySystem)::Int = sum(map(action_size, is.bom))
RLBase.is_terminated(is::InventorySystem) = is.t > is.T
RLBase.reward(is::InventorySystem) = is.reward
print_state(is::InventorySystem) = reduce(vcat, map(print_state, is.bom))
print_action(is::InventorySystem) = reduce(vcat, map(print_action, is.bom))
RLBase.state_space(is::InventorySystem) = state_size(is)
RLBase.action_space(is::InventorySystem) = action_size(is)

function compute_quantities(is,action)
    actions = Iterators.Stateful(action)
    quantity = IdDict{AbstractItem, Vector{Float64}}()
    for item in is.bom
        act_size = action_size(item)
        Qs = Float64[]
        for polparams in partition(Iterators.take(actions, act_size), action_size(item.policy))
            push!(Qs, item.policy(item, polparams...))
        end
        quantity[item] = Qs
    end
    return quantity
end

function (is::InventorySystem)(action::AbstractVector)
    @assert !is_terminated(is) "InventorySystem is at terminal state, please use reset!(::InventorySystem)"
    @assert action_size(is) == length(action) "action must be of length $(action_size(is))"
    quantity = compute_quantities(is, action)
    #Stage: PreActStage
    for item in is.bom
        item(quantity[item])
    end
    #Stage: PreConsStage
    for cons in is.constraints
        cons()
    end
    #Stage: PreDispatchStage
    for item in Iterators.reverse(is.bom)
        dispatch!(item)
    end
    #Stage: PreRewardStage
    is.t += 1
    is.reward = sum(reward!.(is.bom))
end

(is::InventorySystem)(action::AbstractMatrix) = is(vec(action))

function RLBase.reset!(is::InventorySystem)
    reset!.(is.bom)
    reset!.(is.constraints)
    is.t = 1
end

function topological_order(bom)
    L = eltype(bom)[]
    unmarked = copy(bom)
    tempmarked = Set()
    permarked = Set()
    function visit(node)
        if node in permarked
            return nothing
        elseif node in tempmarked
            throw("Node $node was visited multiple times during topological sorting. The BOM contains a cycle.")
        end
        push!(tempmarked, node)
        for child in children(node)
            visit(child)
        end
        push!(permarked, pop!(tempmarked, node))
        push!(L, node)
    end
    while !(isempty(unmarked) && isempty(tempmarked))
        node = pop!(unmarked)
        visit(node)
    end
    return reverse(L)
end

Base.show(io::IO, is::InventorySystem) = print(io, "InventorySystem(", is.T," periods, ", length(is.bom) ," items, (1x$(state_size(is))) state, ($(action_size(is))) action)")


text:
using Configurations
using Configurations: OptionDef, to_dict, to_toml, from_kwargs, from_dict, alias,
    from_toml, no_default, field_defaults, field_default, field_alias, field_aliases,
    option_print, resolve_defaults, PartialDefault
using OrderedCollections
using Test

"Option A"
@option "option_a" struct OptionA
    name::String
    int::Int = 1
end

"Option B"
@option "option_b" struct OptionB
    opt::OptionA = OptionA(;name = "Sam")
    float::Float64 = 0.3
end

dict1 = OrderedDict{String, Any}(
    "opt" => OrderedDict{String, Any}(
        "name" => "Roger",
        "int" => 2,
    ),
    "float" => 0.33
)

dict2 = OrderedDict{String, Any}(
    "float" => 0.33
)

option = from_dict(OptionB, dict1)

@testset "options" begin
    @test option == OptionB(;
        opt = OptionA(;
            name = "Roger",
            int = 2,
        ),
        float = 0.33,
    )

    @test from_toml(OptionB, "option.toml") == option
end

@testset "to_dict" begin
    @test_throws ErrorException to_dict("aaa")
    @test to_dict(option) == dict1
    @test to_dict(from_dict(OptionB, dict2)) == dict2
end

@testset "to_toml" begin
    @test to_toml(option) == "float = 0.33\n\n[opt]\nname = \"Roger\"\nint = 2\n"
    to_toml("test.toml", option)
    @test read("test.toml", String) == "float = 0.33\n\n[opt]\nname = \"Roger\"\nint = 2\n"
end

@testset "default reflection" begin
    @test field_defaults(OptionA) == [no_default, 1]
    @test field_defaults(OptionB) == [OptionA("Sam", 1), 0.3]
    @test field_default(OptionA, :name) === no_default
    @test field_default(OptionA, :int) == 1
end

@testset "from_kwargs" begin
    @test Configurations.from_kwargs(OptionB; opt_name="Roger", opt_int=2, float=0.33) == option
end

@option struct OptionC
    num::Float64

    function OptionC(num::Float64)
        num > 0 || error("not positive")
        new(num)
    end
end

@testset "inside constructor" begin
    @test_throws ErrorException OptionC(-1.0)
end

@option struct OptionD
    opt::Union{OptionA, OptionB}
end

@testset "test multi option type" begin
    d1 = OrderedDict{String, Any}(
        "opt" => OrderedDict{String, Any}(
            "option_b" => dict1
        )
    )

    @test from_dict(OptionD, d1).opt isa OptionB

    d2 = OrderedDict{String, Any}(
        "opt" => OrderedDict{String, Any}(
            "option_a" => OrderedDict{String, Any}(
                "name" => "Roger",
                "int" => 2,
            ),
        )
    )

    @test from_dict(OptionD, d2).opt isa OptionA

    d3 = OrderedDict{String, Any}(
        "opt" => OrderedDict{String, Any}(
            OrderedDict{String, Any}(
                "name" => "Roger",
                "int" => 2,
            ),
        )
    )

    @test_throws ErrorException from_dict(OptionD, d3)
end

@option struct OptionE
    field::Union{Nothing, OptionA} = nothing
end

@testset "optional field" begin
    d = OrderedDict{String, Any}(
        "float" => 0.33
    )

    @test from_dict(OptionB, d) == OptionB(;
        opt = OptionA(;
            name = "Sam",
            int = 1,
        ),
        float = 0.33,
    )

    @test from_kwargs(OptionE) == OptionE(;
        field = nothing,
    )
end

@testset "validate keys" begin
    @test_throws ArgumentError Configurations.validate_keywords(OptionA; abc=2)
    @test Configurations.validate_keywords(OptionB; opt_name="AAA") === nothing
    @test_throws ArgumentError Configurations.validate_keywords(OptionB; opt_abc="AAA")
end

@testset "dict override" begin
    d = OrderedDict{String, Any}(
        "opt" => OrderedDict{String, Any}(
            "name" => "Roger",
            "int" => 2,
        ),
        "float" => 0.33
    )

    @test from_dict(OptionB, d; opt_name = "AAA") == OptionB(;
        opt = OptionA(;
            name = "AAA",
            int = 1,
        ),
        float = 0.33,
    )
end

@option mutable struct Julia
    active::Union{String, Nothing} = nothing
    stable::Union{String, Nothing} = nothing
    # NOTE: we store nightly here too
    versions::Dict{VersionNumber, String} = Dict{VersionNumber, String}()
end

@option mutable struct Ion
    julia::Julia = Julia()
end

@testset "non-dict value conversion" begin
    d = Dict{String, Any}(
        "julia" => Dict{String, Any}(
            "active" => "some/path/to/active",
            "stable" => "1.5.3",
            "versions" => Dict{String, String}(
                "1.5.3" => "some/path/to/1.5.3",
            )
        ),
    )

    @test_throws MethodError from_dict(Ion, d)

    function Configurations.option_convert(::Type{Julia}, ::Type{Dict{VersionNumber,String}}, x::Dict{String,String})
        d = Dict{VersionNumber, String}()
        for (k, v) in x
            d[VersionNumber(k)] = v
        end
        return d
    end

    @test from_dict(Ion, d) == Ion(;
        julia = Julia(;
            active = "some/path/to/active",
            stable = "1.5.3",
            versions = Dict(v"1.5.3" => "some/path/to/1.5.3"),
        ),
    )
end

@testset "to_dict nothing conversion" begin
    d = to_dict(Julia())
    @test !haskey(d, "versions")
    @test !haskey(d, "active")
end

function foo(::Type{T}) where T
    return 1.0
end

function foo(::Type{Int})
    return 2
end

@option struct Inferrable{A, B}
    a::A
    b::B = 1.0
end

@option struct NotInferrable1{T}
    a::Float64
    b::Int
end

@option struct NotInferrable2{T}
    a::T = foo(T)
    b::Int
end

@option struct NotInferrable3{T}
    a::Float64 = foo(T)
    b::Int
end

@testset "parametric types" begin
    @test Inferrable(;a = 1) == Inferrable(1, 1.0)
    # NOTE: we probably should just let this work
    # @test_throws ErrorException field_defaults(Inferrable)
    @test field_defaults(Inferrable{Float64, Float64}) == Any[no_default, 1.0]
    @test field_default(Inferrable{Float64, Float64}, :a) === no_default
    # @test_throws ErrorException field_default(Inferrable, :a) === no_default
    # @test_throws ErrorException field_default(Inferrable, :a)
    @test_throws MethodError NotInferrable1(;a = 1.0, b = 1)
    @test NotInferrable1{Int}(;a = 1.0, b = 1) == NotInferrable1{Int64}(1.0, 1)
    @test_throws MethodError NotInferrable2(;a = 1.0, b = 2)
    @test NotInferrable2{Float64}(;b=2) == NotInferrable2(1.0, 2)
    @test_throws MethodError NotInferrable3(;a = 1.0, b = 1)
    @test field_default(NotInferrable2{Int}, :a) == 2
    @test field_default(NotInferrable2{Float64}, :a) == 1.0
    @test NotInferrable3{Int}(;b = 2) == NotInferrable3{Int}(2.0, 2)
end

@option struct Empty end

@testset "emptry struct" begin
    @test_throws ErrorException field_default(Empty, :name)
end

@option struct CustomKwFn
    a::Float64

    CustomKwFn(a) = new(a)

    function CustomKwFn(; a = 1.0, b = 2.0)
        @test b == 2.0
        new(a)
    end
end

@testset "custom kwfn" begin
    @test CustomKwFn() == CustomKwFn(1.0)
end

@option struct ExtraKwFn
    a::Int = 2

    ExtraKwFn(a) = new(a)

    function ExtraKwFn(; extra=true, kwargs...)
        @test extra
        return Configurations.create(ExtraKwFn; kwargs...)
    end
end

@testset "kwargs forward" begin
    @test ExtraKwFn(;a = 3) == ExtraKwFn(3)
end

@option struct DefaultResolve
    a::Int = 1
    b::Float64 = sin(a)
end

@testset "default resolve" begin
    @test field_default(DefaultResolve, :b) == sin(1)
end

@option struct FieldAlias
    "alpha"
    α::Float64 = 1
    "beta"
    β::Float64
end

@testset "field alias" begin
    d = Dict("alpha" => 2, "beta" => 3)
    @test from_dict(FieldAlias, d) == FieldAlias(;α=2, β=3)
    @test field_aliases(FieldAlias) == ["alpha", "beta"]
end

@testset "non-option type handling" begin
    @test_throws ErrorException field_default(Int, :a)
    @test_throws ErrorException field_alias(Int, :a)
    @test_throws ErrorException alias(Int)
end

@testset "printings" begin
    ex = :(struct OptionA
        name::String
        int::Int = 1
    end)
    def = OptionDef(ex)
    print(def)

    ex = :(struct Inferrable{A, B}
        a::A
        b::B = 1.0
    end)
    def = OptionDef(ex)
    print(def)

    show(stdout, MIME"text/plain"(), FieldAlias(;β=2.0))
    show(stdout, MIME"text/plain"(), OptionB())
    option_print(stdout, MIME"text/plain"(), 1)
    option_print(stdout, MIME"text/plain"(), Dict("a"=>1))
    option_print(stdout, MIME"text/plain"(), [1, 2, 3])
    option_print(stdout, MIME"text/plain"(), rand(2, 2))
end

@option struct UnionToDict
    option::Union{Nothing, OptionA, OptionB} = nothing
end

@testset "multi-union to_dict/from_dict" begin
    d = OrderedDict{String, Any}(
        "option" => OrderedDict{String, Any}(
            "option_a" => to_dict(OptionA(; name="Name"))
        )
    )

    @test from_dict(UnionToDict, d) == UnionToDict(;
        option = OptionA(;
            name = "Name",
        ),
    )

    d = OrderedDict{String, Any}(
        "option" => OrderedDict{String, Any}(
        )
    )

    @test from_dict(UnionToDict, d) == UnionToDict()

    x = UnionToDict(;option=OptionA(; name="Name"))
    d = to_dict(x)
    @test from_dict(UnionToDict, d) == x
end

@option struct DefaultFunction
    a::Float64
    b::Float64 = sin(a)
end

@testset "partial default" begin
    x = DefaultFunction(;a=1.0, b=2.0)
    @test field_default(DefaultFunction, :b)(x) == sin(1.0)
end


text:
function NSDEBase.solve!(solution::MovingWindowSolution, problem, solver::MovingWindowSolver)
    @↓ u0, (t0, tN) ← tspan = problem
    @↓ 𝒫, τ, Δτ = solver
    @↓ 𝒢, P = 𝒫
    for m = 1:length(solution)
        solution[m] = TimeParallelSolution(problem, 𝒫)
        @↓ U, T = solution[m]
        if m == 1
            TimeParallel.coarseguess!(solution[m], problem, u0, t0, t0 + τ, 𝒫)
        else
            ΔP = trunc(Int, P * Δτ / τ)
            N = P - ΔP + 1
            for n = 1:length(T)
                T[n] = solution[m-1].T[n] + Δτ
            end
            for n = 1:N
                U[n] = solution[m-1].U[ΔP+n]
            end
            for n = N:P
                chunk = 𝒢(problem, U[n], T[n], T[n+1])
                U[n+1] = chunk.u[end]
            end
        end
        𝒫(solution[m], problem)
    end
    solution
end

function NSDEBase.solve(problem, solver::MovingWindowSolver)
    solution = MovingWindowSolution(problem, solver)
    solve!(solution, problem, solver)
    solution
end


text:
abstract AbstractPointwise

immutable Identity <: AbstractPointwise end

(ρ::Identity)(x::Any) = x

immutable Log1P{T<:AbstractFloat} <: AbstractPointwise
    threshold::T
    function (::Type{Log1P}){T}(threshold::T)
        (threshold<0) && throw(DomainError)
        new{T}(threshold)
    end
end

(ρ::Log1P)(blob_in::ScatteredBlob) = map(ρ, ScatteredBlob(blob_in.nodes))

function Base.abs{T,N}(nodes::DataStructures.SortedDict{
        Path,Node{T,N},Base.Order.ForwardOrdering})
    realT = real(T)
    absnodes =
        DataStructures.SortedDict{Path,Node{realT,N},Base.Order.ForwardOrdering}()
    for (path, nodevalue) in nodes
        absnodes[path] = Node{realT,N}(abs.(nodevalue.data), nodevalue.ranges)
    end
    return absnodes::DataStructures.SortedDict{
        Path,Node{realT,N},Base.Order.ForwardOrdering}
end

Base.abs(Wx::ScatteredBlob) = ScatteredBlob(abs(Wx.nodes))

function Base.map(ρ::AbstractPointwise, innodes::SortedDict)
    outnodes =
        DataStructures.SortedDict{Path,Node,Base.Order.ForwardOrdering}()
    for path in keys(innodes)
        outnodes[path] = ρ(innodes[path])
    end
    return outnodes
end

function Base.map!(
        ρ::AbstractPointwise,
        blob_out::ScatteredBlob,
        blob_in::ScatteredBlob)
    @inbounds for id in eachindex(blob_in.nodes)
        map!(ρ, blob_out.nodes[id].data, blob_in.nodes[id].data)
    end
end


text:
module TestMisc

import LinearAlgebra
using BenchmarkCI
using BenchmarkCI: as_https, as_target_url, metadata_from
using PkgBenchmark: BenchmarkConfig
using Test

@testset "versionof" begin
    @test (@test_logs (:error,) BenchmarkCI.versionof(Base)) === nothing
    @test (@test_logs (:error,) BenchmarkCI.versionof(LinearAlgebra)) === nothing
    @test (@test_logs (:error,) BenchmarkCI.versionof(LinearAlgebra.BLAS)) === nothing
end

@testset "format_period" begin
    @test BenchmarkCI.format_period(3) == "3 seconds"
    @test BenchmarkCI.format_period(125) == "2 minutes 5 seconds"
end

@testset "metadata" begin
    if !Sys.iswindows()
        metadata = mktempdir() do path
            fakejulia = joinpath(path, "julia")
            symlink(Base.julia_cmd().exec[1], fakejulia)
            metadata_from(;
                target = BenchmarkConfig(),
                baseline = BenchmarkConfig(juliacmd = `$fakejulia`),
                pkgdir = :dummy,
                script = :dummy,
                project = :dummy,
            )
        end
        @test metadata[:target_julia_info] != nothing
        @test metadata[:baseline_julia_info] != nothing
    end
end

@testset "as_https" begin
    @test as_https("git@github.com:USER/REPO.git") === "https://github.com/USER/REPO"
    @test as_https("git@github.com:USER/REPO") === "https://github.com/USER/REPO"
    @test as_https("git@github.com:USER/REPO.jl.git") === "https://github.com/USER/REPO.jl"
    @test as_https("git@github.com:USER/REPO.jl") === "https://github.com/USER/REPO.jl"
    @test as_https("spam") === nothing
end

@testset "as_target_url" begin
    @test as_target_url(
        "git@github.com:USER/REPO.git",
        "benchmark-results",
        "2020/07/27/023806",
    ) === "https://github.com/USER/REPO/blob/benchmark-results/2020/07/27/023806/result.md"
    @test as_target_url("spam", "benchmark-results", "2020/07/27/023806") === nothing
end

@testset "error_on_missing_github_token" begin
    err = nothing
    @test try
        BenchmarkCI.error_on_missing_github_token()
        false
    catch err
        true
    end
    @test occursin("`GITHUB_TOKEN` is not set", sprint(showerror, err))
end

@testset "compress_tar" begin
    mktempdir() do dir
        srcdir = joinpath(dir, "src")
        mkdir(srcdir)
        write(joinpath(srcdir, "README.md"), "hello")

        tarfile = joinpath(dir, "dest.tar.zst")
        BenchmarkCI.compress_tar(tarfile, srcdir)

        @test isfile(tarfile)

        destdir = joinpath(dir, "dest")
        BenchmarkCI.decompress_tar(destdir, tarfile)
        @test read(joinpath(destdir, "README.md"), String) == "hello"
    end
end

end  # module


text:
# need conversions for things like:
# a = BInt{1:10,Int8}[ 1, 2, 3, 4 ]
# a = BInt{1:10}[ 1, 2, 3, 4 ]
# a[2] = 4

@generated function Base.convert{R}( ::Type{ Bool }, x::BInt{R} )
        l,u = get_range_bounds( R )
        if l > 0 || u < 0
                # range does not include 0 - so must be true
                return true
        elseif l == u == 0
                return false
        else
                return :( Bool(x.n) )
        end
end

# conversion from Integer to BInt:
@generated function Base.convert{R,T1<:Integer}( ::Type{ BInt{R} }, x::T1 )
        T = default_int_type( R )
        Tworking = default_int_type( range_union( T, T1 ) )
        l,u = get_range_bounds( R, Tworking )
        return :(
                if ( $Tworking(x) >= $l ) ;
                        if ( $Tworking(x) <= $u ) ;
                                return UncheckedBInt{R,$T}( x )
                        end;
                end;
                throw( InexactError() );
        )
end

@generated function Base.convert{R,T,T1<:Integer}( ::Type{ BInt{R,T} }, x::T1 )
        l,u = get_range_bounds( R, T )
        return :(
                if ( $T(x) >= $l ) ;
                        if ( $T(x) <= $u ) ;
                                return UncheckedBInt{R,T}( x )
                        end;
                end;
                throw( InexactError() );
        )
end

# conversion from BInt to Integer:
Base.convert{R1,T1,T<:Integer}( ::Type{ T }, x::BInt{R1,T1} ) = T( x.n )

# for conversion from BInt to BInt, only allow narrower ranges to
# be assigned to wider ranges
@generated function Base.convert{R1,R2,T2}( ::Type{ BInt{R2} }, x::BInt{R1,T2} )
        if !range_within_range( R1, R2 )
                throw( InexactError() );
        end
        :( UncheckedBInt{ R2 }( x.n ) )
end
@generated function Base.convert{R1,T1,R2,T2}( ::Type{ BInt{R2,T2} }, x::BInt{R1,T1} )
        if !range_within_range( R1, R2 )
                throw( InexactError() );
        end
        :( UncheckedBInt{ R2, T2 }( x.n ) )
end

# converting from a BInt to a narrower range BInt
# is a special case - want to make the runtime check
# be explicit
# so allow
#           Int -> BInt conversion
#           BInt -> wider BInt conversion
#           BInt -> narrower BInt - only with explicit tighten_range fn.

@generated function tighten_range{R,R2,T2}( ::Type{BInt{R}}, x::BInt{R2,T2} )
        T = default_int_type( R )
        l,u = get_range_bounds( R, T )

        :(
                if ( $T(x.n) < $l ) ; throw( InexactError() ); end;
                if ( $T(x.n) > $u ) ; throw( InexactError() ); end;
                BInt{R,$T}( $T(x.n), 0 )
        )
end

# this is a little lazy - there are corner cases where it will fail
# but Bint * Int * Int cannot be protected from overflow, so its a lost
# cause - stay within the BInt type, or else risk overflow

# this is because of things like promote_type( Int64, UInt64 ) which
# gives a UInt64.  While it is true that an Int64 can be losslessly converted
# to a UInt64 and back again, it doesn't mean that arithmetic will work
Base.promote_rule{R1,T1,T2<:Integer }(::Type{BInt{R1,T1}}, ::Type{T2} ) = promote_type( T1,T2)

Base.promote_rule{R1,R2}(::Type{ BInt{R1} },::Type{ BInt{R2} } ) = begin
        R = range_union( R1, R2 )
        T = default_int_type( R )
        return BInt{R,T}
end
Base.promote_rule{R1,T1,R2}(::Type{ BInt{R1,T1} },::Type{ BInt{R2} } ) = begin
        R = range_union( R1, R2 )
        T = default_int_type( R )
        return BInt{R,T}
end
Base.promote_rule{R1,T1,R2,T2}(::Type{ BInt{R1,T1} },::Type{ BInt{R2,T2} } ) = begin
        R = range_union( R1, R2 )
        T = default_int_type( R )
        return BInt{R,T}
end


text:
using ScottishTaxBenefitModel
using Documenter

makedocs(;
    modules=[ScottishTaxBenefitModel],
    authors="Graham Stark",
    repo="https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/{commit}{path}#L{line}",
    sitename="ScottishTaxBenefitModel.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://grahamstark.github.io/ScottishTaxBenefitModel.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Module Listing" => "module_listing.md",
        "Weighting" => "weighting.md",
        "Annotated Bibliography" => "annotated_bibliography.md",
        "Test Suite" => "test_suite.md",
        "Validation in Aggregate" => "validation_in_aggregate.md",
        "Coding Notes" => "coding_notes.md",
        "TODO" => "TODO.md"
    ],
)

deploydocs(;
    repo="github.com/grahamstark/ScottishTaxBenefitModel.jl",
)

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using OPCSPs
using MCTS
using ProfileView
using POMDPs
import POMDPToolbox

p = gen_informative(noise_2norm=100.0,
                n_nodes=10,
                n_connected=4,
                rbar=10.0,
                rng=MersenneTwister(1))
is = rand!(MersenneTwister(1543), create_state(p), initial_belief(p))
s = GurobiExactSolver(multithreaded=false)
# s = HeuristicSolver()
solver = AgUCTSolver(
    aggregator=VoronoiOPCSPAg(50.0),
    rollout_solver=FeedbackSolver(s),
    exploration_constant=100.0,
    n_iterations=10000,
    rng=MersenneTwister(1947)
)

policy = MCTSAdapter(POMDPs.solve(solver, OPCSPBeliefMDP(p)))
sim = POMDPToolbox.HistoryRecorder(rng=MersenneTwister(2011), initial_state=is)
simr = simulate(sim, p, policy, OPCSPUpdater(p), initial_belief(p))

Profile.clear()

policy = MCTSAdapter(POMDPs.solve(solver, OPCSPBeliefMDP(p)))
@profile simulate(sim, p, policy, OPCSPUpdater(p), initial_belief(p))

policy = MCTSAdapter(POMDPs.solve(solver, OPCSPBeliefMDP(p)))
@time simulate(sim, p, policy, OPCSPUpdater(p), initial_belief(p))

ProfileView.view()


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

# ## 12 steps to Navier–Stokes
# =====
# ***

# The final two steps in this interactive module teaching beginning [CFD with Python](https://bitbucket.org/cfdpython/cfd-python-class) will both solve the Navier–Stokes equations in two dimensions, but with different boundary conditions.

# The momentum equation in vector form for a velocity field $\vec{v}$ is:

# $$\frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v}=-\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v}$$

# This represents three scalar equations, one for each velocity component $(u,v,w)$. But we will solve it in two dimensions, so there will be two scalar equations.

# Remember the continuity equation? This is where the [Poisson equation](./13_Step_10.ipynb) for pressure comes in!

# Step 11: Cavity Flow with Navier–Stokes
# ----
# ***

# Here is the system of differential equations: two equations for the velocity components $u,v$ and one equation for pressure:

# $$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right) $$

# $$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) $$

# $$\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)$$

# From the previous steps, we already know how to discretize all these terms. Only the last equation is a little unfamiliar. But with a little patience, it will not be hard!

# ### Discretized equations

# First, let's discretize the $u$-momentum equation, as follows:

# $$
# \begin{split}
# & \frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{u_{i,j}^{n}-u_{i-1,j}^{n}}{\Delta x}+v_{i,j}^{n}\frac{u_{i,j}^{n}-u_{i,j-1}^{n}}{\Delta y} = \\
# & \qquad -\frac{1}{\rho}\frac{p_{i+1,j}^{n}-p_{i-1,j}^{n}}{2\Delta x}+\nu\left(\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{\Delta x^2}+\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{\Delta y^2}\right)
# \end{split}
# $$

# Similarly for the $v$-momentum equation:

# $$
# \begin{split}
# &\frac{v_{i,j}^{n+1}-v_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{v_{i,j}^{n}-v_{i-1,j}^{n}}{\Delta x}+v_{i,j}^{n}\frac{v_{i,j}^{n}-v_{i,j-1}^{n}}{\Delta y} = \\
# & \qquad -\frac{1}{\rho}\frac{p_{i,j+1}^{n}-p_{i,j-1}^{n}}{2\Delta y}
# +\nu\left(\frac{v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}}{\Delta x^2}+\frac{v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}}{\Delta y^2}\right)
# \end{split}
# $$

# Finally, the discretized pressure-Poisson equation can be written thus:

# $$
# \begin{split}
# & \frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2} = \\
# & \qquad \rho \left[ \frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right) -\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} - 2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x} - \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
# \end{split}
# $$

# You should write these equations down on your own notes, by hand, following each term mentally as you write it.

# As before, let's rearrange the equations in the way that the iterations need to proceed in the code. First, the momentum equations for the velocity at the next time step.

# The momentum equation in the $u$ direction:

# $$
# \begin{split}
# u_{i,j}^{n+1} = u_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(u_{i,j}^{n}-u_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(u_{i,j}^{n}-u_{i,j-1}^{n}\right) \\
# & - \frac{\Delta t}{\rho 2\Delta x} \left(p_{i+1,j}^{n}-p_{i-1,j}^{n}\right) \\
# & + \nu \left(\frac{\Delta t}{\Delta x^2} \left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)\right)
# \end{split}
# $$

# The momentum equation in the $v$ direction:

# $$
# \begin{split}
# v_{i,j}^{n+1} = v_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(v_{i,j}^{n}-v_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(v_{i,j}^{n}-v_{i,j-1}^{n})\right) \\
# & - \frac{\Delta t}{\rho 2\Delta y} \left(p_{i,j+1}^{n}-p_{i,j-1}^{n}\right) \\
# & + \nu \left(\frac{\Delta t}{\Delta x^2} \left(v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}\right)\right)
# \end{split}
# $$

# Almost there! Now, we rearrange the pressure-Poisson equation:

# $$
# \begin{split}
# p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2\left(\Delta x^2+\Delta y^2\right)} \\
# & -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \\
# & \times \left[\frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right)-\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} -2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x}-\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
# \end{split}
# $$

# The initial condition is $u, v, p = 0$ everywhere, and the boundary conditions are:

# $u=1$ at $y=2$ (the "lid");

# $u, v=0$ on the other boundaries;

# $\frac{\partial p}{\partial y}=0$ at $y=0$;

# $p=0$ at $y=2$

# $\frac{\partial p}{\partial x}=0$ at $x=0,2$

# Implementing Cavity Flow
# ----

using Plots

# TODO Unicode
nx = 41
ny = 41
nt = 500
nit = 50
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
x = LinRange(0, 2, nx)
y = LinRange(0, 2, ny)
X, Y = numpy.meshgrid(x, y)

rho = 1
nu = .1
dt = .001

u = zeros((ny, nx))
v = zeros((ny, nx))
p = zeros((ny, nx))
b = zeros((ny, nx))

# The pressure Poisson equation that's written above can be hard to write out without typos.
# The function `build_up_b` below represents the contents of the square brackets, so that the entirety of the PPE is slightly more manageable.

function build_up_b(b, rho, dt, u, v, dx, dy)

    b[2:end, 2:end] = (rho * (1 / dt *
                    ((u[2:end, 3:end] - u[2:end, begin:end-1]) /
                     (2 * dx) + (v[2:, 2:end] - v[begin:end-1, 2:end]) / (2 * dy)) -
                    ((u[2:end, 3:end] - u[2:end, begin:end-1]) / (2 * dx))^2 -
                      2 * ((u[2:, 2:end] - u[begin:end-1, 2:end]) / (2 * dy) *
                           (v[2:end, 3:end] - v[2:end, begin:end-1]) / (2 * dx))-
                          ((v[2:, 2:end] - v[begin:end-1, 2:end]) / (2 * dy))^2))

    return b
end

# The function `pressure_poisson` is also defined to help segregate the different rounds of calculations.
# Note the presence of the pseudo-time variable `nit`.
# This sub-iteration in the Poisson calculation helps ensure a divergence-free field.

function pressure_poisson(p, dx, dy, b)
    pn = empty(p)
    pn = copy(p)

    for q in 1:nit
        pn = copy(p)
        p[2:end, 2:end] = (((pn[2:end, 3:end] + pn[2:end, begind:end-1]) * dy^2 +
                          (pn[2:, 2:end] + pn[begin:end-1, 2:end]) * dx^2) /
                          (2 * (dx^2 + dy^2)) -
                          dx^2 * dy^2 / (2 * (dx^2 + dy^2)) *
                          b[2:end,2:end])

        p[:, end] = p[:, end-1] # dp/dx = 0 at x = 2
        p[begin, :] = p[2, :]   # dp/dy = 0 at y = 0
        p[:, end] = p[:, 2]   # dp/dx = 0 at x = 0
        p[begin, :] = 0        # p = 0 at y = 2
    end

    return p
end

# Finally, the rest of the cavity flow equations are wrapped inside the function `cavity_flow`,
# allowing us to easily plot the results of the cavity flow solver for different lengths of time.

function cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)
    un = empty(u)
    vn = empty(v)
    b = zeros((ny, nx))

    for n in range(nt):
        un = copy(u)
        vn = copy(v)

        b = build_up_b(b, rho, dt, u, v, dx, dy)
        p = pressure_poisson(p, dx, dy, b)

        u[2:end, 2:end] = (un[2:end, 2:end]-
                         un[2:end, 2:end] * dt / dx *
                        (un[2:end, 2:end] - un[2:end, begin:end-1]) -
                         vn[2:end, 2:end] * dt / dy *
                        (un[2:end, 2:end] - un[begin:end-1, 2:end]) -
                         dt / (2 * rho * dx) * (p[2:end, 3:end] - p[2:end, begin:end-1]) +
                         nu * (dt / dx^2 *
                        (un[2:end, 3:end] - 2 * un[2:end, 2:end] + un[2:end, begin:end-1]) +
                         dt / dy^2 *
                        (un[3:end, 2:end] - 2 * un[2:end, 2:end] + un[begin:end-1, 2:end])))

        v[2:end,2:end] = (vn[2:end, 2:end] -
                        un[2:end, 2:end] * dt / dx *
                       (vn[2:end, 2:end] - vn[2:end, begin:end-1]) -
                        vn[2:end, 2:end] * dt / dy *
                       (vn[2:end, 2:end] - vn[begin:end-1, 2:end]) -
                        dt / (2 * rho * dy) * (p[3:end, 2:end] - p[begin:end-1, 2:end]) +
                        nu * (dt / dx^2 *
                       (vn[2:end, 3:end] - 2 * vn[2:end, 2:end] + vn[2:end, begin:end-1]) +
                        dt / dy^2 *
                       (vn[3:end, 2:end] - 2 * vn[2:end, 2:end] + vn[begin:end-1, 2:end])))

        u[begin, :]  = 0
        u[:, begin]  = 0
        u[:,   end] = 0
        u[end,   :] = 1    # set velocity on cavity lid equal to 1
        v[begin, :]  = 0
        v[end,   :] = 0
        v[:, begin]  = 0
        v[:,   end] = 0
    end

    return u, v, p
end

# Let's start with `nt = 100` and see what the solver gives us:

u = zeros((ny, nx))
v = zeros((ny, nx))
p = zeros((ny, nx))
b = zeros((ny, nx))
nt = 100
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)

# TODO PLOT
fig = pyplot.figure(figsize=(11,7), dpi=100)
# plotting the pressure field as a contour
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
pyplot.colorbar()
# plotting the pressure field outlines
pyplot.contour(X, Y, p, cmap=cm.viridis)
# plotting velocity field
pyplot.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2])
pyplot.xlabel('X')
pyplot.ylabel('Y');

# You can see that two distinct pressure zones are forming and that the spiral pattern expected from lid-driven cavity flow is beginning to form.  Experiment with different values of `nt` to see how long the system takes to stabilize.

u = zeros((ny, nx))
v = zeros((ny, nx))
p = zeros((ny, nx))
b = zeros((ny, nx))
nt = 700
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)

# TODO PLOT
fig = pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
pyplot.colorbar()
pyplot.contour(X, Y, p, cmap=cm.viridis)
pyplot.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2])
pyplot.xlabel('X')
pyplot.ylabel('Y');

# The quiver plot shows the magnitude of the velocity at the discrete points in the mesh grid we created.
# (We're actually only showing half of the points because otherwise it's a bit of a mess.
# TODO FIX
# The `X[::2, ::2]` syntax above is a convenient way to ask for every other point.)

# Another way to visualize the flow in the cavity is to use a `streamplot`:

# TODO PLOT
fig = pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
pyplot.colorbar()
pyplot.contour(X, Y, p, cmap=cm.viridis)
pyplot.streamplot(X, Y, u, v)
pyplot.xlabel('X')
pyplot.ylabel('Y');

## Learn More

# The interactive module **12 steps to Navier–Stokes** is one of several components of the Computational Fluid Dynamics class taught by Prof. Lorena A. Barba in Boston University between 2009 and 2013.

# For a sample of what the other components of this class are, you can explore the **Resources** section of the Spring 2013 version of [the course's Piazza site](https://piazza.com/bu/spring2013/me702/resources).

# ***

max_stars_repo_name: miguelraz/CFDPython

max_issues_repo_name: miguelraz/CFDJulia

max_forks_repo_name: miguelraz/CFDJulia


text:
using Documenter
using CyclicArrays

makedocs(
    sitename = "CyclicArrays",
    format = Documenter.HTML(),
    modules=[CyclicArrays],
    repo="git@github.com:udistr/CyclicArrays.jl.git",
    authors="udistr <strobach.ehud@gmail.com>",
    pages = [
    "Home" => "index.md",
    "Main" => "main.md",
            ]
)

deploydocs( repo = "github.com:udistr/CyclicArrays.jl" )


text:
#    Unit signalProcessing.jl, part of PosDefManifold Package for julia language
#   v 0.3.5 - last update 6th of November 2019
#
#    MIT License
#    Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#    https://sites.google.com/site/marcocongedo/home
#
#    DESCRIPTION
#    This Unit implements Signal Processing methods and prcedures
#    useful in relation to the Riemannian Geometry of the manifold
#    of Symmetric Positive Definite (SPD) or Hermitian matrices

"""
    randChi²(df::Int)

 **alias**: `randχ²`

 Generate a random variable distributed as a *chi-squared* with `df`
 degrees of freedom.

 It uses the *Wilson–Hilferty transformation* for `df`>=20 -
 see [chi-squared distribution](https://en.wikipedia.org/wiki/Chi-squared_distribution).

 ## Examples
    using Plots, PosDefManifold
    chi=[randχ²(2) for i=1:10000]
    histogram(chi) # needs Plots package. Check your plots back-end.
"""
randChi²(df::Int) =
    df<20 ? sum(randn()^2 for i=1:df) : df*(1.0-2.0/(9.0*df)+randn()*sqrt2/sqrt(9.0*df))^3

randχ²=randChi²

"""
```
    randEigvals(n::Int;
    <
    df::Int=2,
    eigvalsSNR::Real=10e3 >)
```

 **alias**: `randλ`

 Generate an ``n``-vector of random real positive eigenvalues.
 The eigenvalues are generated as in function `randΛ`([`randEigvalsMat`](@ref)),
 the syntax of which is used.

 **See also**: `randU` ([`randUnitaryMat`](@ref)), `randP` ([`randPosDefMat`](@ref)).

 ## Examples
    using Plots, PosDefManifold
    λ=sort(randλ(10), rev=true)
    σ=sort(randλ(10, eigvalsSNR=10), rev=true)
    plot(λ) # needs Plots package. Check your plots back-end.
    plot!(σ) # needs Plots package. Check your plots back-end.

"""
randEigvals(n::Int;
            df::Int=2,
            eigvalsSNR::Real=10e3) =
    eigvalsSNR==Inf ? [randχ²(df) for i in 1:n] : [randχ²(df)+(df/eigvalsSNR) for i in 1:n]

randλ=randEigvals

"""
```
    (1) randEigvalsMat(n::Int;
    <
    df::Int=2,
    eigvalsSNR::Real=10e3 >)

    (2) randEigvalsMat(n::Int, k::Int;
    < same keyword arguments as in (1) >)
```

 **alias**: `randΛ`

 (1) Generate an ``n⋅n`` diagonal matrix of random real positive eigenvalues.
 (2) An array 1d (of [𝔻Vector type](@ref)) of ``k`` matrices of the kind in (1)

 The eigenvalues are generated according to model

 ``λ_i=χ_{df}^2+η,\\hspace{6pt}\\textrm{for}\\hspace{2pt}i=1:n,``

 where
 - ``χ_{df}^2`` (signal term) is randomly distributed as a [chi-square](https://bit.ly/1IXkulE) with `df` degrees of freedom,
 - ``η`` is a [white noise](https://bit.ly/2TN8472) term, function of *<keyword argument>* `eigvalsSNR`, such that

 ``\\textrm{eigenvalues SNR}=\\mathbb{E}\\big(\\sum_{i=1}^{n}λ_i\\big)\\big/nη.``

 The expected sum ``\\mathbb{E}\\big(\\sum_{i=1}^{n}λ_i\\big)`` here above is the
 expected variance of the signal term, i.e., ``n(df)``, since the expectation
 of a random chi-squared variable is equal to its degrees of freedom.

 If `eigvalsSNR=Inf` is passed as argument, then ``η`` is set to zero, *i.e.*,
 no white noise is added. In any case `eigvalsSNR` must be positive.

 Note that with the default value of *<keyword argument>* `df` (`df=2`)
 the generating model assumes that the eigenvalues
 have exponentially decaying variance, which is often observed on real data.

!!! note "Nota Bene"
    The *<keyword argument>* `eigvalsSNR` expresses the expected eigenvalues SNR
    ([signal-to-noise ratio](https://bit.ly/1VvpvnQ)),
    not the real one, and is not expressed in decibels,
    but as the expected SNR variance ratio.

 This function is used by function `randP` ([`randPosDefMat`](@ref)) to generate
 random positive definite matrices with added white noise in order
 to emulate eigenvalues observed in real data and to
 improve the conditioning of the generated matrices with respect to inversion.

 **See also**: `randλ` ([`randEigvals`](@ref)), `randU` ([`randUnitaryMat`](@ref)),
 `randP` ([`randPosDefMat`](@ref)), `randχ²` ([`randChi²`](@ref)).

 ## Examples
    using PosDefManifold
    # (1)
    n=3;
    U=randU(n);
    Λ=randΛ(n, eigvalsSNR=100)
    P=U*Λ*U' # generate an SPD matrix
    using LinearAlgebra
    Q=ℍ(U*Λ*U') # generate an SPD matrix and flag it as 'Hermitian'

    # (2) generate an array of 10 matrices of simulated eigenvalues
    Dvec=randΛ(n, 10)

"""
randEigvalsMat(n::Int;
                df::Int=2,
                eigvalsSNR::Real=10e3) =
        𝔻(randλ(n, df=df, eigvalsSNR=eigvalsSNR))

randEigvalsMat(n::Int, k::Int;
                df::Int=2,
                eigvalsSNR::Real=10e3) =
        𝔻Vector([randΛ(n, df=df, eigvalsSNR=eigvalsSNR) for j=1:k])

randΛ=randEigvalsMat

"""
    (1) randUnitaryMat(n::Int)
    (2) randUnitaryMat(::Type{Complex{T}}, n::Int)

 **aliases**: `randOrthMat`, `randU`

 Generate a random ``n⋅n``
 - (1) [orthogonal](https://bit.ly/2vrr0wU) matrix (real)
 - (2) [unitary](https://bit.ly/2JCHbmC) matrix (complex)

 The matrices are generated running the modified (stabilized)
 [Gram-Schmidt orthogonalization](https://bit.ly/2YE6zvy)
 procedure ([`mgs`](@ref)) on an ``n⋅n`` matrix filled with random Gaussian elements.

 **See also**: `randΛ` ([`randEigvals`](@ref)), `randP` ([`randPosDefMat`](@ref)).

 ## Examples
    using PosDefManifold
    n=3;
    X=randU(n)*sqrt(randΛ(n))*randU(n)'  # (1) generate a random square real matrix

    U=randU(ComplexF64, n);
    V=randU(ComplexF64, n);
    Y=U*sqrt(randΛ(n))*V' # (2) generate a random square complex matrix

"""
randUnitaryMat(n::Int)=mgs(randn(Float64, n, n))

randOrthMat(n::Int)=mgs(randn(Float64, n, n))

randUnitaryMat(::Type{Complex{T}}, n::Int) where {T<:AbstractFloat} = mgs(randn(ComplexF64, n, n))

randU=randUnitaryMat

"""
```
    (1) randPosDefMat(n::Int;
    <
    df::Int=2,
    eigvalsSNR::Real=10e3 >)

    (2) randPosDefMat(::Type{Complex{T}}, n:: Int;
    < same keyword arguments as in (1) >)

    (3) randPosDefMat(n::Int, k::Int;
    <
    df::Int=2,
    eigvalsSNR::Real=10e3,
    SNR::Real=100,
    commuting=false >)

    (4) randPosDefMat(::Type{Complex{T}}, n::Int, k::Int;
    < same keyword arguments as in (3) >)
```

 **alias**: `randP`

 Generate
 - (1) one random `Hermitian` positive definite matrix (real) of size ``n⋅n``
 - (2) one random `Hermitian` positive definite matrix (complex) of size ``n⋅n``
 - (3) an array 1d (of [ℍVector type](@ref)) of ``k`` matrices of the kind in (1)
 - (4) an array 1d (of [ℍVector type](@ref)) of ``k`` matrices of the kind in (2).

 Methods (3) and (4) are **multi-threaded**. See [Threads](@ref).

 For (1) and (2) the matrix is generated according to model

 ``UΛU^H+ηI``,

 where ``U`` is a random orthogonal (1) or unitary (2) matrix generated by
 function `randU`([`randUnitaryMat`](@ref)) and ``Λ``, ``η`` are a positive definite
 diagonal matrix and a non-negative scalar depending on *<optional keywords arguments>*
 `df` and `eigvalsSNR` randomly generated calling function
 `randΛ`([`randEigvalsMat`](@ref)).

 For (3) and (4), if the *<optional keyword argument>* `commuting=true` is passed,
 the ``k`` matrices are generated according to model

  ``UΛ_iU^H+ηI,\\hspace{8pt}``, for ``i``=1:``k``

 otherwise they are generated accrding to model

 ``(UΛ_iU^H+ηI)+φ(V_iΔ_iV_i^H+ηI),\\hspace{8pt}``, for ``i``=1:``k``  Eq.[1]

 where
 - ``U`` and the ``V_i`` are random (3) orthogonal/(4) unitary matrices,
 - ``Λ_i`` and ``Δ_i`` are positive definite diagonal matrices
 - ``η`` is a non-negative scalar.
  All variables here above are randomly generated as in (1) and (2)
 - ``φ`` is adjusted so as to obtain a desired output `SNR` ([signal-to-noise ratio](https://bit.ly/1VvpvnQ)), which is also an
 *<optional keywords arguments>*, such as

 ``SNR=\\frac{\\displaystyle\\sum_{i=1}^{k}\\textrm{tr}(UΛ_iU^H+ηI)}{\\displaystyle\\sum_{i=1}^{k}\\textrm{tr}φ(V_iΔ_iV_i^H+ηI)}``.

!!! note "Nota Bene"
    The keyword arguments `SNR` is not expressed in decibels,
    but as the expected SNR variance ratio. It must be a positive number.

 A slightly different version of this model for generating positive definite
 matrices has been proposed in (Congedo et *al.*, 2017b)[🎓];
 in the model of Eq. [1]
 - ``UΛ_iU^H`` is the signal term, where the signal is supposed sharing the same coordinates for all matrices,
 - ``φ(V_iΔ_iV_i^H)`` is a structured noise term, which is different for all matrices
 - ``ηI`` is a [white noise](https://bit.ly/2TN8472) term, with same variance for all matrices.

 **See also**: the aforementioned paper and `randΛ` ([`randEigvalsMat`](@ref)).

 ## Examples
    using PosDefManifold
    R=randP(10, df=10, eigvalsSNR=1000) # 1 SDP Matrix of size 10x10 #(1)
    H=randP(ComplexF64, 5, eigvalsSNR=10) # 1 Hermitian Matrix of size 5x5 # (2)
    ℛ=randP(10, 1000, eigvalsSNR=100) # 1000 SPD Matrices of size 10x10 # (3)
    using Plots
    heatmap(Matrix(ℛ[1]), yflip=true, c=:bluesreds)
    ℋ=randP(ComplexF64, 20, 1000) # 1000 Hermitian Matrices of size 20x20 # (4)

"""
randPosDefMat(n::Int;
              df::Int=2,
              eigvalsSNR::Real=10e3) =
    congruence(randU(n), randΛ(n, df=df, eigvalsSNR=eigvalsSNR), ℍ)

randPosDefMat(::Type{Complex{T}}, n::Int;
             df::Int=2,
             eigvalsSNR::Real=10e3)     where {T<:AbstractFloat} =
    congruence(randU(ComplexF64, n), randΛ(n, df=df, eigvalsSNR=eigvalsSNR), ℍ)

function randPosDefMat(n::Int, k::Int;
             df::Int=2,
             eigvalsSNR::Real=10e3,
             SNR::Real=100,
             commuting=false)
    U=randU(n)
    𝐏=ℍVector(undef, k)
    φ=1/SNR
    if commuting
        @threads for j=1:k 𝐏[j]=ℍ(cong(U, randΛ(n, df=df, eigvalsSNR=eigvalsSNR), 𝕄)) end
    else
        𝐐=𝕄Vector(undef, k)
        @threads for j=1:k 𝐐[j]=randU(n) end
        @threads for j=1:k 𝐏[j]=ℍ(cong(U, randΛ(n, df=df, eigvalsSNR=eigvalsSNR), 𝕄)
                                +cong(𝐐[j], (φ*randΛ(n, df=df, eigvalsSNR=eigvalsSNR)), 𝕄)) end
    end
    return 𝐏
end

function randPosDefMat(::Type{Complex{T}}, n::Int, k::Int;
            df::Int=2,
            eigvalsSNR::Real=10e3,
            SNR::Real=100,
            commuting=false)            where {T<:AbstractFloat}
    U=randU(ComplexF64, n)
    𝐏=ℍVector(undef, k)
    φ=1/SNR
    if commuting
        @threads for j=1:k 𝐏[j]=ℍ(cong(U, randΛ(n, df=df, eigvalsSNR=eigvalsSNR), 𝕄)) end
    else
        𝐐=𝕄Vector(undef, k)
        @threads for j=1:k 𝐐[j]=randU(ComplexF64, n) end
        @threads for j=1:k 𝐏[j]=ℍ(cong(U, randΛ(n, df=df, eigvalsSNR=eigvalsSNR), 𝕄)
                            +cong(𝐐[j], (φ*randΛ(n, df=df, eigvalsSNR=eigvalsSNR)), 𝕄)) end
    end
    return 𝐏
end

randP=randPosDefMat

"""
    (1) regularize!(P::ℍ; <SNR=10e3>)
    (2) regularize!(𝐏::ℍVector; <SNR=10e3>)

 Add [white noise](https://bit.ly/2TN8472) to either
 - (1) a positive definite matrix ``P`` of size ``n⋅n``, or
 - (2) a 1d array ``𝐏`` of ``k`` positive definite matrices of size ``n⋅n``, of [ℍVector type](@ref).

 The added noise improves the matrix conditioning with respect to inversion.
 This is used to avoid numerical errors when decomposing these matrices
 or when evaluating some functions of their eigevalues such as the log.

 A constant value is added to all diagonal elements of (1) ``P``
 or (2) af all matrices in ``𝐏``,
 that is, on output:

 ``\\textrm{(1)}\\hspace{2pt}P\\leftarrow P+ηI``

 ``\\textrm{(2)}\\hspace{2pt}𝐏_i\\leftarrow 𝐏_i+ηI, \\hspace{2pt}\\textrm{for}\\hspace{2pt} i=1:k.``

 The amount of added noise ``η`` is determined by the `SNR`
 *<keyword argument>*, which by default is 10000. This is
 such that

 ``\\textrm{(1)}\\hspace{2pt}SNR=\\frac{\\displaystyle\\textrm{tr}(P)}{\\displaystyle\\textrm{tr}(ηI)}.``

 ``\\textrm{(2)}\\hspace{2pt}SNR=\\frac{\\displaystyle\\sum_{i=1}^{k}\\textrm{tr}(𝐏_i)}{\\displaystyle k\\hspace{1pt}\\textrm{tr}(ηI)}.``

 ``P`` in (1) must be flagged as Hermitian. See [typecasting matrices](@ref).

!!! note "Nota Bene"
    The keyword argument `SNR` expresses a SNR ([signal-to-noise ratio](https://bit.ly/1VvpvnQ)),
    and is not expressed in decibels,  but as the SNR variance ratio.
    It must be a positive number. Differently from function
    `randΛ`[`randEigvalsMat`](@ref), `randλ`[`randEigvals`](@ref) and
    `randP`[`randPosDefMat`](@ref), the SNR here is not the expected SNR,
    but the actual SNR.

**See also**: `randP` ([`randPosDefMat`](@ref)).

 ## Examples
    # (1)
    using LinearAlgebra, Plots, PosDefManifold
    n=3
    U=randU(n)
    # in Q we will write two matrices,
    # the unregularized and regularized matrix side by side
    Q=Matrix{Float64}(undef, n, n*2)
    P=ℍ(U*Diagonal(randn(n).^2)*U') # generate a real 3x3 positive matrix
    for i=1:n, j=1:n Q[i, j]=P[i, j] end
    regularize!(P, SNR=5)
    for i=1:n, j=1:n Q[i, j+n]=P[i, j] end # the regularized matrix is on the right
    heatmap(Matrix(Q), yflip=true, c=:bluesreds)

    # (2)
    𝐏=[ℍ(U*Diagonal(randn(3).^2)*U') for i=1:5] # 5 real 3x3 positive matrices
    regularize!(𝐏, SNR=1000)

 ## Run a test
    using LinearAlgebra
    𝐏=randP(10, 100, SNR=1000); # 100 real Hermitian matrices
    signalVar=sum(tr(P) for P in 𝐏);
    regularize!(𝐏, SNR=1000);
    signalPlusNoiseVar=sum(tr(P) for P in 𝐏);
    output_snr=signalVar/(signalPlusNoiseVar-signalVar)
    # output_snr should be approx. equal to 1000

"""
function regularize!(P::ℍ;
                    SNR=10e3)
    n=size(P, 1)
    η=tr(P)/(SNR*n)
    for i in 1:n P[i, i]+=η  end
end

function regularize!(𝐏::ℍVector;
                    SNR=10e3)
    k=length(𝐏)
    n=size(𝐏[1], 1)
    η=sum(tr(P) for P in 𝐏)/(SNR*n*k)
    for l in 1:k, i in 1:n 𝐏[l][i, i]+=η  end
end

"""
    gram(X::Matrix{T}) where T<:RealOrComplex

 Given a generic data matrix ``X``, comprised of real or complex elements,
 return the normalized [Gram matrix](https://bit.ly/2I0FQn2), that is,
 the covariance matrix of ``X``
 corrected by sample size, but without subtracting the mean.

 The result is flagged as `Hermitian`.
 See [typecasting matrices](@ref).

!!! note "Nota Bene"
    If ``X`` is wide or square (r<=c) return ``XX^H/c``.
    If ``X`` is tall (r>c)            return ``X^HX/r``.

 ## Examples
    using PosDefManifold
    X=randn(5, 150);
    G=gram(X) # => G=X*X'/150
    X=randn(100, 2);
    F=gram(X); # => G=X'*X/100
"""
function gram(X::𝕄{T}) where T<:RealOrComplex
    (r, c)=size(X)
    r<c ? ℍ((X*X')/c) : ℍ((X'*X)/r)
end # function gram

"""
trade(P::ℍ{T}) where T<:RealOrComplex

 Given a positive definite matrix `P`, return as a 2-tuple the
 *trace* and the *determinant* of `P`.
 This is used to plot positive matrices in two dimensions
 (*TraDe plots*: log(trace/n) vs. log(determinant), see exemple here below).

 `P` must be flagged by julia as `Hermitian`.
  See [typecasting matrices](@ref).

 ### Examples
    using PosDefManifold
    P=randP(3)
    t, d=trade(P)  # equivalent to (t, d)=trade(P)

    # TraDe plot
    using Plots
    k=100
    n=10
    𝐏=randP(n, k, SNR=1000); # 100 real Hermitian matrices
    x=Vector{Float64}(undef, k)
    y=Vector{Float64}(undef, k)
    for i=1:k
        x[i], y[i] = trade(𝐏[i])
    end
    x=log.(x./n)
    y=log.(y)
    plot(x, y, seriestype=:scatter)
"""
trade(P::ℍ{T}) where T<:RealOrComplex = (tr(P), det(P))


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(The BTG Program)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#################################################################################
#copy this box into a file to use the BTG program                               #
using LinearAlgebra                                                             #
using Test                                                                      #
using Distributions                                                             #
using DataFrames                                                                #
using CSV                                                                       #
using Polynomials                                                               #
                                                                                #
include("computation/finitedifference.jl")                                      #
include("quadrature/quadrature.jl")                                             #
include("transforms/transforms.jl")                                             #
include("priors/priors.jl")                                                     #
include("various_mle_things/incremental.jl")                                              #
include("kernels/kernel.jl")                                                    #
include("datastructs.jl")                                                       #
include("computation/buffers0.jl") #datastruct, kernel, incremental, quadrature #
include("model0.jl") #buffers, datastructs, several auxiliary                   #
include("computation/tdist.jl") #model0 and buffer0                             #
#################################################################################

############################### Load in data
df = DataFrame(CSV.File("datasets/abalone.csv"))
data = convert(Matrix, df[:,2:8]) #length, diameter, height, whole weight, shucked weight, viscera weight, shell weight
target = convert(Array, df[:, 9]) #age
normalizing_constant = maximum(target)
target = target/normalizing_constant #normalization

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(Problem Setup)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ind = 11:200
posx = 1:7 #locations
posc = 1:1 #covariates
x = data[ind, posx]
#choose a subset of variables to be regressors for the mean
Fx = data[ind, posc]
y = float(target[ind])
pind = 10:10 #prediction index
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
#prediction data
d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
Fx0 = reshape(data[pind, posc], 1, length(posc)) #Pind associate with test data
x0 = reshape(data[pind, posx], 1, length(posx))
#rangeθ = [10.0 20]
#rangeθ = [10.0 1000.0; 10.0 1000]
rangeθ = [1000.0 10000.0]
rangeλ = [-3 3.0] #we will always used 1 range scale for lambda
#################################################################################
#init btg object for testing various functions
btg1 = btg(trainingData1, rangeθ, rangeλ)

θ1 = btg1.nodesWeightsθ.nodes[1, 6] #pick some theta value which doubles as quadrature node

if true #test solve_btg
    #output = @capture_err begin
        (dpdf, pdf, cdf) = solve(btg1)
    #end;
    dpdf_fixed = y0 -> dpdf(x0, Fx0, y0)
    pdf_fixed = y0 -> pdf(x0, Fx0, y0)
    cdf_fixed = y0 -> cdf(x0, Fx0, y0)
    a = dpdf_fixed
    b = pdf_fixed
    c = cdf_fixed

    #(r1, r2, plt1, pol) = checkDerivative(c, b, 0.5, nothing, 1, 17, 20) #function first, derivative second
    #display(plt1)
    #println(pol)
    #display(plt1)
    if true
    #plt(a, .01, 1, 150)
    println("Plotting pdf...")
    (x, y) = plt(b, .001, 1, 200, label = "pdf")
    println("Plotting cdf...")
    (x1, y1) = plt!(c, 0.001, 1, 200, label = "cdf")
    Plots.plot!(target[pind], seriestype = :vline, title = "PDF and CDF of Bayesian Predictive Distribution", label = "actual")
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
function compare_stack(stack, nums; vm = 3, dvm = 2, cmap = "seismic_wide_y")

    fig, axes = plt.subplots(2, length(nums), sharex = true, sharey = true)
    ref = stack[:, :, end-nums[1]]

    for (idx, n) in enumerate(nums)
        axes[1, idx].imshow(stack[:, :, end-n], vmax = vm, vmin = -vm, cmap = cmap)
        axes[2, idx].imshow(stack[:, :, end-n] .- ref, vmax = dvm, vmin = -dvm, cmap = cmap)
    end

    return fig, axes
end