text:
facts("Lots of tests") do
    context("First group") do
        @fact 1   => 1
        @fact 2*2 => 4
        @fact uppercase("foo") => "FOO"
    end

    context("Second group") do
        @fact_throws 2^-1
        @fact_throws DomainError 2^-1
        @fact_throws DomainError 2^-1 "a nifty message"
        @fact 2*[1,2,3] => [2,4,6]
    end
end


text:
using NESSie.Format
using JSON
using LightXML: parse_string, root, name, child_elements, attribute, content

@testset "writexml3d_json (nodes)" begin
    for T in testtypes
        # empty system
        js = JSON.parse(readback(fh -> writexml3d_json(fh, Vector{T}[])))
        @test js["format"] == "xml3d-json"
        @test haskey(js, "version") == true
        @test js["data"]["position"]["type"] == "float3"
        @test length(js["data"]["position"]["seq"]) == 1
        @test js["data"]["position"]["seq"][1]["value"] == []
        # small system
        nodes = Vector{T}[T[0, 0, 0], T[0, 0, 3], T[0, 3, 0], T[1, -3, 3]]
        js = JSON.parse(readback(fh -> writexml3d_json(fh, nodes)))
        @test js["format"] == "xml3d-json"
        @test haskey(js, "version") == true
        @test js["data"]["position"]["type"] == "float3"
        @test length(js["data"]["position"]["seq"]) == 1
        @test js["data"]["position"]["seq"][1]["value"] == [0, 0, 0, 0, 0, 3, 0, 3, 0, 1, -3, 3]
    end
end

@testset "writexml3d_json (surface model)" begin
    for T in testtypes
        # empty system
        js = JSON.parse(readback(fh -> writexml3d_json(fh, Model{T, Triangle{T}}())))
        @test js["format"] == "xml3d-json"
        @test haskey(js, "version") == true
        @test js["data"]["index"]["type"] == "int"
        @test length(js["data"]["index"]["seq"]) == 1
        @test js["data"]["index"]["seq"][1]["value"] == []
        @test js["data"]["position"]["type"] == "float3"
        @test length(js["data"]["position"]["seq"]) == 1
        @test js["data"]["position"]["seq"][1]["value"] == []
        @test js["data"]["normal"]["type"] == "float3"
        @test length(js["data"]["normal"]["seq"]) == 1
        @test js["data"]["normal"]["seq"][1]["value"] == []
        # small system
        nodes    = Vector{T}[T[0, 0, 0], T[0, 0, 3], T[0, 3, 0], T[1, -3, 3]]
        elements = [Triangle(nodes[1], nodes[2], nodes[3]),
                    Triangle(nodes[1], nodes[4], nodes[2])]
        js = JSON.parse(readback(fh -> writexml3d_json(fh, Model(nodes, elements))))
        @test js["format"] == "xml3d-json"
        @test haskey(js, "version") == true
        @test js["data"]["index"]["type"] == "int"
        @test length(js["data"]["index"]["seq"]) == 1
        @test js["data"]["index"]["seq"][1]["value"] == [0, 1, 2, 0, 3, 1]
        @test js["data"]["position"]["type"] == "float3"
        @test length(js["data"]["position"]["seq"]) == 1
        @test js["data"]["position"]["seq"][1]["value"] == [0, 0, 0, 0, 0, 3, 0, 3, 0, 1, -3, 3]
        @test js["data"]["normal"]["type"] == "float3"
        @test length(js["data"]["normal"]["seq"]) == 1
        @test map(T, js["data"]["normal"]["seq"][1]["value"]) ≈ map(T,
            [√360 \ [-9 - √90, -3, 0]; √360 \ [-9 - √90, -3, 0]; [-1, 0, 0]; √90 \ [-9, -3, 0]]
        )
    end
end

@testset "writexml3d_xml" begin
    for T in testtypes
        # empty system
        xroot = root(parse_string(readback(fh -> writexml3d_xml(fh, Vector{T}[]))))
        @test name(xroot) == "xml3d"
        xchildren = collect(child_elements(xroot))
        @test length(xchildren) == 1
        @test name(xchildren[1]) == "data"
        @test attribute(xchildren[1], "id") == "mesh"
        xchildren = collect(child_elements(xchildren[1]))
        @test length(xchildren) == 1
        @test name(xchildren[1]) == "float3"
        @test strip(content(xchildren[1])) == ""
        # small system
        nodes    = Vector{T}[T[0, 0, 0], T[0, 0, 3], T[0, 3, 0], T[1, -3, 3]]
        xroot = root(parse_string(readback(fh -> writexml3d_xml(fh, nodes))))
        @test name(xroot) == "xml3d"
        xchildren = collect(child_elements(xroot))
        @test length(xchildren) == 1
        @test name(xchildren[1]) == "data"
        @test attribute(xchildren[1], "id") == "mesh"
        xchildren = collect(child_elements(xchildren[1]))
        @test length(xchildren) == 1
        @test name(xchildren[1]) == "float3"
        @test T[parse(T, e) for e in split(strip(content(xchildren[1])))] == T[0, 0, 0, 0, 0, 3, 0, 3, 0, 1, -3, 3]
    end
end

ext: jl


text:
function json(response::HODSyncResponse, client)

end

function get_result(response::HODAsyncResponse, client)

end


text:
using Example1

@show Point(1,2) + Point(3,4)

q = MPoint(1,2)
@show q

q.x = 9
@show q

@show Polynomial([1,2,3])
@show evaluate(Polynomial([1,2,3]),2)

@show 3* Polynomial([1,2,4,42])

#@show deriv(Polynomial([1,1,1]))

#@show deriv(Polynomial([0,0,1]))

Example1.greet()


text:
# Functions:
# reportfeasible
# subout
# back

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

function reportfeasible(logicset::LogicalCombo;
    command=""::String,
    verbose=true::Bool)::LogicalCombo

    !verbose && return logicset
    filler = repeat("\t", max(1, 3-Integer(floor(length(command)/8))))

    Nfeas = nfeasible(logicset)
    Perceived = percievedfeasible(logicset)

    (Nfeas == 0)  && (check = "X ")
    (Nfeas  > 1)  && (check = "✓ ")
    (Nfeas == 1)  && (check = "✓✓")

    ender = (Nfeas>0) ? ":" *
      join(logicset[rand(1:Nfeas),:,:], " ") : " [empty set]"

    println(" $filler Feasible Outcomes: $Nfeas \t Perceived Outcomes: $Perceived $check \t $ender")
    return logicset
end

reportfeasible(command::String, verbose::Bool)::Function =
  (x -> reportfeasible(x, command=command, verbose=verbose))

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

function subout(txtcmd, i, arg, mykeys)
  lookup(vect, i) = i ∈ 1:length(vect) ? vect[i] : "~~OUTOFBOUNDS~~"

  (arg[end] ∈ ['i', 'j', 'J'])  && return replace(txtcmd, "{{$arg}}"=>lookup(mykeys,i))

  mod = integer(match(r"([0-9]+$)", arg).match)

  occursin("+", arg) && return replace(txtcmd, "{{$arg}}"=>lookup(mykeys,i+mod))
  occursin("-", arg) && return replace(txtcmd, "{{$arg}}"=>lookup(mykeys,i-mod))

  txtcmd
end

subout("{{i + 1}}", 1, "i + 1", [:a,:b,:c,:d])

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

function stringnumparse(x)
    !occursin(r"\+|\-",x) && return integer(x)

    x2 = split(x, r"\+|\-")

    occursin("+", x) && return +(integer.(x2)...)
    occursin("-", x) && return -(integer.(x2)...)

    throw("$x not interpreted!")
end

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

Base.:+(x::String, y::String) = x * " " * y

function skipthischeck(txt, i, j; verbose=false)
    txt = replace(txt, " "=>"")

    occursin(r"!(i|j|J)", txt) && return i==j

    occursin(r"![0-9]+$", txt) &&
     return integer(match(r"^!([0-9])+$", txt).captures[1]) == j

    occursin(r"^!(i|j|J)$", txt) && (i  == j) && return true

      m = match(r"^\*?([>=<]{0,2})(.*)([>=<]{0,2})\*?$", txt)
      if m.match != nothing
        centernum = stringnumparse(strip.(m.captures[2]))
        verbose && print(" (" +
         (m.captures[1] == "" ? "" : string(j) + string(m.captures[1])) *
           " $centernum " *
         (m.captures[3] == "" ? "" : string(m.captures[3]) + string(j)) * ") ")
        # * [operator] center
        (m.captures[1] == ">")  && !(j >  centernum    ) && return true
        (m.captures[1] == ">=") && !(j >= centernum    ) && return true
        (m.captures[1] == ">>") && !(j >  centernum + 1) && return true
        (m.captures[1] == "<")  && !(j <  centernum    ) && return true
        (m.captures[1] == "<=") && !(j <= centernum    ) && return true
        (m.captures[1] == "<<") && !(j <  centernum - 1) && return true

        # center [operator] *
        (m.captures[3] == "<")  && !(centernum     <  j) && return true
        (m.captures[3] == "<=") && !(centernum     <= j) && return true
        (m.captures[3] == "<<") && !(centernum - 1 << j) && return true
        (m.captures[3] == ">")  && !(centernum     >  j) && return true
        (m.captures[3] == ">=") && !(centernum     >= j) && return true
        (m.captures[3] == ">>") && !(centernum + 1 >> j) && return true

        (m.captures[1] ∈ ["=","=="]) && (centernum != j) && return true
        (m.captures[3] ∈ ["=","=="]) && (centernum != j) && return true
      end
    false
end

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

function splitthenskipcheck(txt, i, j; verbose=false)
    rangematchessplit =
       strip.(split(replace(txt, r"i|j|J"=>i), r"(,|\|)"))

    occursin("|", txt) && return all(skipthischeck.(rangematchessplit, i, j, verbose=verbose))

    !occursin("|", txt) && return any(skipthischeck.(rangematchessplit, i, j, verbose=verbose))

    false
end

# txt = "!1,!3"
#
# for i in 1:5, j in 1:5
#     print("i = $i, j = $j skipthis {{$txt}} Skip? ")
#     println(splitthenskipcheck(txt, i, j, verbose=true))
# end
#

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

function back(logicset::LogicalCombo; verbose=false)
  (length(logicset.commands)==0) && (println("Nothing to go back to!"); return logicset)
  commandlist = copy(logicset.commands)
  pop!(commandlist)
  return logicalparse(commandlist, verbose=verbose)
end


text:
using PseudoBoolean
using Test

include("variable_utils.jl")
include("qubo.jl")


text:
using MotifPvalue
using Test

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


text:
using JuliaStudies
using Test

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


text:
macro pretty(json)
    return esc(:(JSON3.pretty($json)))
end

pretty(str; kw...) = pretty(stdout, str; kw...)
pretty(out::IO, x) = pretty(out, JSON3.write(x))
function pretty(out::IO, str::String, indent=0, offset=0; kw...)
    buf = codeunits(str)
    len = length(buf)
    if len == 0
        return
    end
    pos = 1
    b = getbyte(buf, pos)
    @wh
    # printing object?
    if b == UInt8('{')
        Base.write(out, b)
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b == UInt8('}')
            Base.write(out, b)
            return
        end
        Base.write(out, '\n')
        indent += 1
        keys = []
        vals = []
        # loop thru all key-value pairs, keeping track of longest key to pad others
        while b != UInt8('}')
            pos, str = JSON3.read(StructTypes.StringType(), buf, pos, len, b, String; kw...)
            push!(keys, str)
            @eof
            b = getbyte(buf, pos) # ':'
            @wh
            pos += 1
            @eof
            b = getbyte(buf, pos)
            @wh
            pos, x = JSON3.read(StructTypes.Struct(), buf, pos, len, b, Any; kw...)
            push!(vals, x)
            @eof
            b = getbyte(buf, pos) # ',' or '}'
            @wh
            if b == UInt8('}')
                break
            end
            pos += 1
            @eof
            b = getbyte(buf, pos)
            @wh
        end
        maxlen = maximum(map(sizeof, keys)) + 5
        # @show maxlen
        for i = 1:length(keys)
            Base.write(out, "  "^indent)
            Base.write(out, lpad("\"$(keys[i])\"" * ": ", maxlen + offset, ' '))
            pretty(out, JSON3.write(vals[i]; kw...), indent, maxlen + offset; kw...)
            if i == length(keys)
                indent -= 1
                Base.write(out, "\n" * ("  "^indent * " "^offset) * "}")
            else
                Base.write(out, ",\n")
            end
        end

    # printing array?
    elseif b == UInt8('[')
        Base.write(out, b)
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b == UInt8(']')
            Base.write(out, b)
            return
        end
        Base.write(out, '\n')
        indent += 1
        vals = []
        while b != UInt8(']')
            pos, x = JSON3.read(StructTypes.Struct(), buf, pos, len, b, Any; kw...)
            push!(vals, x)
            @eof
            b = getbyte(buf, pos) # ',' or ']'
            @wh
            if b == UInt8(']')
                break
            end
            pos += 1
            @eof
            b = getbyte(buf, pos)
            @wh
        end
        for (i, val) in enumerate(vals)
            Base.write(out, "  "^indent * " "^offset)
            pretty(out, JSON3.write(vals[i]; kw...), indent, offset; kw...)
            if i == length(vals)
                indent -= 1
                Base.write(out, "\n" * ("  "^indent * " "^offset) * "]")
            else
                Base.write(out, ",\n")
            end
        end

    # printing constant?
    else
        Base.write(out, str)
    end
    return
@label invalid
    Base.error("error pretty-fying json: $error")
end


text:
function jt11a()
#
#  test polynomial suggested by Jenkins and Traub
#
    m = 15;
    y = exp.((im*pi) .* (1-m:m-1) ./ (2*m))
    z = 0.9*exp.((im*pi) .*(m:3*m) ./ (2*m))
    p1 = poly(y)
    p2 = poly(z)
    p = reverse(conv(p1.a, p2.a))
    z = [[y; z] ones(4*m)]
    p, PolyZeros(z)
end


text:
include("Functions.jl")

abstract type AbstractRobot
end

import HorizonSideRobots: move!, isborder, putmarker!, ismarker

move!(robot::AbstractRobot, side::HorizonSide) = HorizonSideRobots.move!(get(robot), side)
isborder(robot::AbstractRobot,  side::HorizonSide) = HorizonSideRobots.isborder(get(robot), side)
putmarker!(robot::AbstractRobot) = HorizonSideRobots.putmarker!(get(robot))
ismarker(robot::AbstractRobot) = HorizonSideRobots.ismarker(get(robot))
get(robot::AbstractRobot) = nothing

function snake(robot::AbstractRobot)
    side = Ost
    movements!(robot,side)
    while !isborder(robot, Nord)
        move!(robot,Nord)
        side = reverse_side(side)
        movements!(robot,side)
    end
end

#-------------------------------------------

mutable struct ChessRobot <: AbstractRobot
    robot::Robot
    flag::Int
    ChessRobot(robot::Robot)=new(robot, 1)
end

function move!(robot::ChessRobot, side::HorizonSide)
    move!(get(robot), side)
    robot.flag = (robot.flag - 1) * (-1) # 1 -> 0, 0 -> 1
    if robot.flag == 1
        putmarker!(robot)
    end
end

function movements!(robot::ChessRobot, side::HorizonSide)
    while !isborder(robot, side)
        move!(robot, side)
    end
end

get(robot::ChessRobot) = robot.robot

#-------------------------------------------

mutable struct Coord
    x::Int
    y::Int
    Coord() = new(0,0)
    Coord(x::Int,y::Int) = new(x,y)
end

function move!(coord::Coord, side::HorizonSide)
    if side==Nord
        coord.y += 1
    elseif side==Sud
        coord.y -= 1
    elseif side==Ost
        coord.x += 1
    else
        coord.x -= 1
    end
end

get_coord(coord::Coord) = (coord.x, coord.y)

#-------------------------------------------

mutable struct CoordRobot <: AbstractRobot
    robot::Robot
    coord::Coord
    CoordRobot(r,(x, y))=new(r,Coord(0,0))
    CoordRobot(r)=new(r,Coord())
end

get(robot::CoordRobot) = robot.robot

function move!(robot::CoordRobot, side::HorizonSide)
    move!(get(robot), side)
    move!(robot.coord, side)
end

#-------------------------------------------

abstract type AbstractCoordRobot <: AbstractRobot
end

function move!(robot::AbstractCoordRobot, side::HorizonSide)
    move!(get(robot), side)
    move!(robot.coord, side)
end

#------------------------------

struct CounterMarkersRobot <: AbstractCoordRobot
    robot::Robot
    num_markers::Int
    CounterMarkersRobot(robot::Robot)=new(robot, if ismarker(robot) 1 else 0 end)
end

get(robot::CounterMarkersRobot)=robot.robot

function move!(robot::CounterMarkersRobot, side::HorizonSide)
    move!(get(robot), side)
    if ismarker(robot)
        robot.num_markers += 1
    end
end

#------------------------------

abstract type AbstractBorderRobot <: AbstractRobot end

try_move!(robot::AbstractBorderRobot,side::HorizonSide) = try_move!(get(robot), side)

function movements!(robot::AbstractBorderRobot, side::HorizonSide)
    n = 0
    while try_move!(robot,side)
        n += 1
    end
    return n
end

function movements!(robot::AbstractBorderRobot, side::HorizonSide, num_steps::Integer)
    for _ in 1:num_steps
        try_move!(robot,side)
    end
end

function snake(robot::AbstractBorderRobot)
    side = Ost
    movements!(robot,side)
    while !isborder(robot, Nord)
        try_move!(robot,Nord)
        # putmarker!(robot )
        side = reverse_side(side)
        movements!(robot,side)
    end
end

#------------------------------

struct PutmarkersRobot <: AbstractBorderRobot
    robot::Robot
end

function try_move!(robot::PutmarkersRobot, side::HorizonSide)
    if try_move!(get(robot),side)
        putmarker!(robot)
        return true
    end
    return false
end

get(robot::PutmarkersRobot) = robot.robot

#------------------------------

struct MarkerBorderRobot <: AbstractBorderRobot
    robot::Robot
end

function marker_to_edge!(robot::MarkerBorderRobot, side::HorizonSide)
    path = movements!(robot, side)
    putmarker!(robot)
    rside = reverse_side(side)
    movements!(robot, rside, path)
end

get(robot::MarkerBorderRobot) = robot.robot

#------------------------------

mutable struct  CountbordersRobot <: AbstractBorderRobot
    robot::Robot
    state::Bool
    count::Int
    border_side::HorizonSide
    CountbordersRobot(r, border_side=Nord) = new(r,0,0, border_side)
end

get(robot::CountbordersRobot) = robot.robot

function try_move!(robot::CountbordersRobot, side::HorizonSide)
    result = try_move!(get(robot), side)
    if isborder(robot.robot, robot.border_side)
        if robot.state == 0
            robot.state = 1
            robot.count += 1
        end
    else
        robot.state = 0
    end
    return result
end

get_num_borders(robot::CountbordersRobot) = robot.count

#------------------------

abstract type BoundlessRobot <: AbstractRobot end

try_move!(robot::BoundlessRobot,side::HorizonSide) = try_move!(get(robot), side)

function moves!(r::BoundlessRobot, side::HorizonSide, n::Int)
    while (n > 0 && isborder(r, reverse_side(left(side))))
        move!(r, side)
        n -= 1
    end
end

function move_through_hole!(r::BoundlessRobot, side::HorizonSide)
    step = 1
    cur_step = step
    bypass_side = left(side)
    while isborder(r, side)
        moves!(r, bypass_side, cur_step)
        cur_step += step
        bypass_side = reverse_side(bypass_side)
    end
    move!(r, side)
    while isborder(r, bypass_side)
        move!(r, side)
    end
    if cur_step != 1
        for i in 0:div((cur_step + 1), 2)
            move!(r, bypass_side)
        end
    end
end

#------------------------------

struct MarkerBoundlessRobot <: BoundlessRobot
    robot::Robot
end

function movements!(r::MarkerBoundlessRobot, side::HorizonSide, n::Int)
    while (n > 0 && !ismarker(r))
        move_through_hole!(r, side)
        n -= 1
    end
end

function move!(r::MarkerBoundlessRobot, side::HorizonSide)
    if !ismarker(r)
        move!(get(r), side)
    end
end

get(robot::MarkerBoundlessRobot) = robot.robot

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
# Diffusion maps
# --------------
# Diffusion maps,
# Coifman, R. & Lafon, S., Applied and Computational Harmonic Analysis, Elsevier, 2006, 21, 5-30

#### DiffMap type
struct DiffMap{T <: AbstractFloat} <: SpectralResult
    t::Int
    ɛ::Float64
    K::AbstractMatrix{T}
    proj::Projection{T}

    DiffMap{T}(t::Int, ɛ::T, K::AbstractMatrix{T}, proj::Projection{T})  where T = new(t, ɛ, K, proj)
end

## properties
outdim(M::DiffMap) = size(M.proj, 1)
projection(M::DiffMap) = M.proj
kernel(M::DiffMap) = M.K

## show & dump
function show(io::IO, M::DiffMap)
    print(io, "Diffusion Maps(outdim = $(outdim(M)), t = $(M.t), ɛ = $(M.ɛ))")
end

function Base.dump(io::IO, M::DiffMap)
    println(io, "Dimensionality:")
    show(io, outdim(M))
    print(io, "\n\n")
    println(io, "Timesteps:")
    show(io, M.t)
    print(io, "\n\n")
    println(io, "Kernel: ")
    Base.showarray(io, M.K, false, header=false)
    println(io)
    println(io, "Embedding:")
    Base.showarray(io, projection(M), false, header=false)
end

## interface functions
function transform(::Type{DiffMap}, X::DenseMatrix{T};
                   d::Int=2, t::Int=1, ɛ::T=1.0) where T<:AbstractFloat
    transform!(fit(UnitRangeTransform, X), X)

    sumX = sum(X.^ 2, dims=1)
    K = exp.(( transpose(sumX) .+ sumX .- 2*transpose(X) * X ) ./ ɛ)

    p = transpose(sum(K, dims=1))
    K ./= ((p * transpose(p)) .^ t)
    p = transpose(sqrt.(sum(K, dims=1)))
    K ./= (p * transpose(p))

    U, S, V = svd(K, full=true)
    U ./= U[:,1]
    Y = U[:,2:(d+1)]

    return DiffMap{T}(t, ɛ, K, transpose(Y))
end


text:
module AlgebraicTraits

# useful because Complex is a leaf type in Julia
typealias AbstractComplex Union{Real, Complex}

include("associativity.jl")
include("closure.jl")
include("identity.jl")

include("monoid.jl")

end  # module AlgebraicTraits


text:
export StockTradingEnv, StockTradingEnvWithTurbulence

using Pkg.Artifacts
using DelimitedFiles
using LinearAlgebra:dot
using IntervalSets

function load_default_stock_data(s)
    if s == "prices.csv" || s == "features.csv"
        data, _ = readdlm(joinpath(artifact"stock_trading_data", s), ',', header=true)
        collect(data')
    elseif s == "turbulence.csv"
        readdlm(joinpath(artifact"stock_trading_data", "turbulence.csv")) |> vec
    else
        @error "unknown dataset $s"
    end
end

mutable struct StockTradingEnv{F<:AbstractMatrix{Float64}, P<:AbstractMatrix{Float64}} <: AbstractEnv
    features::F
    prices::P
    HMAX_NORMALIZE::Float32
    TRANSACTION_FEE_PERCENT::Float32
    REWARD_SCALING::Float32
    initial_account_balance::Float32
    state::Vector{Float32}
    total_cost::Float32
    day::Int
    first_day::Int
    last_day::Int
    daily_reward::Float32
end

_n_stocks(env::StockTradingEnv) = size(env.prices, 1)
_prices(env::StockTradingEnv) = @view(env.state[2:1+_n_stocks(env)])
_holds(env::StockTradingEnv) = @view(env.state[2+_n_stocks(env):_n_stocks(env)*2+1])
_features(env::StockTradingEnv) = @view(env.state[_n_stocks(env)*2+2:end])
_balance(env::StockTradingEnv) = @view env.state[1]
_total_asset(env::StockTradingEnv) = env.state[1] + dot(_prices(env), _holds(env))

"""
    StockTradingEnv(;kw...)

This environment is originally provided in [Deep Reinforcement Learning for Automated Stock Trading: An Ensemble Strategy](https://github.com/AI4Finance-LLC/Deep-Reinforcement-Learning-for-Automated-Stock-Trading-Ensemble-Strategy-ICAIF-2020)

# Keyword Arguments

- `initial_account_balance=1_000_000`.
"""
function StockTradingEnv(;
    initial_account_balance=1_000_000f0,
    features=nothing,
    prices=nothing,
    first_day=nothing,
    last_day=nothing,
    HMAX_NORMALIZE = 100f0,
    TRANSACTION_FEE_PERCENT = 0.001f0,
    REWARD_SCALING = 1f-4
)
    prices = isnothing(prices) ? load_default_stock_data("prices.csv") : prices
    features = isnothing(features) ? load_default_stock_data("features.csv") : features

    @assert size(prices, 2) == size(features, 2)

    first_day = isnothing(first_day) ? 1 : first_day
    last_day = isnothing(last_day) ? size(prices, 2) : last_day
    day = first_day

    # [balance, stock_prices..., stock_holds..., features...]
    state = zeros(Float32, 1 + size(prices, 1) * 2 + size(features, 1))

    env = StockTradingEnv(
        features,
        prices,
        HMAX_NORMALIZE,
        TRANSACTION_FEE_PERCENT,
        REWARD_SCALING,
        initial_account_balance,
        state,
        0f0,
        day,
        first_day,
        last_day,
        0f0
    )

    _balance(env)[] = initial_account_balance
    _prices(env) .= @view prices[:, day]
    _features(env) .= @view features[:, day]

    env
end

function (env::StockTradingEnv)(actions)
    init_asset = _total_asset(env)

    # sell first
    for (i, s) in enumerate(actions)
        if s < 0
            sell = min(-env.HMAX_NORMALIZE * s, _holds(env)[i])
            _holds(env)[i] -= sell
            gain = _prices(env)[i] * sell
            cost = gain * env.TRANSACTION_FEE_PERCENT
            _balance(env)[] += gain - cost
            env.total_cost += cost
        end
    end

    # then buy
    # better to shuffle?
    for (i,b) in enumerate(actions)
        if b > 0
            max_buy = div(_balance(env)[], _prices(env)[i])
            buy = min(b*env.HMAX_NORMALIZE, max_buy)
            _holds(env)[i] += buy
            deduction = buy * _prices(env)[i]
            cost = deduction * env.TRANSACTION_FEE_PERCENT
            _balance(env)[] -= deduction + cost
            env.total_cost += cost
        end
    end

    env.day += 1
    _prices(env) .= @view env.prices[:, env.day]
    _features(env) .= @view env.features[:, env.day]

    env.daily_reward = _total_asset(env) - init_asset
end

RLBase.reward(env::StockTradingEnv) = env.daily_reward * env.REWARD_SCALING
RLBase.is_terminated(env::StockTradingEnv) = env.day >= env.last_day
RLBase.state(env::StockTradingEnv) = env.state

function RLBase.reset!(env::StockTradingEnv)
    env.day = env.first_day
    _balance(env)[] = env.initial_account_balance
    _prices(env) .= @view env.prices[:, env.day]
    _features(env) .= @view env.features[:, env.day]
    env.total_cost = 0.
    env.daily_reward = 0.
end

RLBase.state_space(env::StockTradingEnv) = Space(fill(-Inf32..Inf32, length(state(env))))
RLBase.action_space(env::StockTradingEnv) = Space(fill(-1f0..1f0, length(_holds(env))))

RLBase.ChanceStyle(::StockTradingEnv) = DETERMINISTIC

# wrapper

struct StockTradingEnvWithTurbulence{E<:StockTradingEnv} <: AbstractEnvWrapper
    env::E
    turbulences::Vector{Float64}
    turbulence_threshold::Float64
end

function StockTradingEnvWithTurbulence(;
    turbulence_threshold=140.,
    turbulences=nothing,
    kw...
)
    turbulences = isnothing(turbulences) && load_default_stock_data("turbulence.csv")

    StockTradingEnvWithTurbulence(
        StockTradingEnv(;kw...),
        turbulences,
        turbulence_threshold
    )
end

function (w::StockTradingEnvWithTurbulence)(actions)
    if w.turbulences[w.env.day] >= w.turbulence_threshold
        actions .= ifelse.(actions .< 0, -Inf32, 0)
    end
    w.env(actions)
end


text:
function recsolve!(x, y, lu, rd, memo)
    any((x, y) .> size(memo[])) && return 0
    ismissing(memo[][x, y]) || return memo[][x, y]

    if x == size(memo[], 1)
        memo[][x, y] = recsolve!(x, y + 1, lu, rd, memo)
    elseif y == size(memo[], 2)
        memo[][x, y] = recsolve!(x + 1, y, lu, rd, memo)
    else
        rs = recsolve!(x, y + 1, lu, rd, memo) + recsolve!(x + 1, y, lu, rd, memo)
        memo[][x, y] = rs / 2
    end

    memo[][x, y]
end

function solve(case, wh, lu, rd)
    memo = Array{Union{Missing,Rational}}(missing, wh)
    memo[lu[1]:rd[1], lu[2]:rd[2]] .= Rational(0)
    memo[wh...] = Rational(1)

    p = recsolve!(1, 1, lu, rd, Ref(memo)) |> Float64
    println("Case #$(case): $p")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        W, H, L, U, R, D = parse.(Int, split(readline(stdin)))
        solve(case, (W, H), (L, U), (R, D))
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
module chapters_installation

using Dash, DashHtmlComponents, DashCoreComponents
include("../../reusable_components/Header.jl")

app = dash()

app.layout = html_div() do
    Header("Dash Installation"),
    dcc_markdown("""

    In order to use Dash, please ensure that you are using a version of Julia >= 1.2.

    To install the most recent version:

    ```
        pkg> add Dash DashCoreComponents DashHtmlComponents DashTable
    ```

    To install the latest (stable) development version of Dash instead:

    ```
        using Pkg;
        Pkg.add(PackageSpec(url="https://github.com/plotly/DashBase.jl.git"))
        Pkg.add(PackageSpec(url="https://github.com/plotly/dash-html-components.git", rev="master"))
        Pkg.add(PackageSpec(url="https://github.com/plotly/dash-core-components.git", rev="master"))
        Pkg.add(PackageSpec(url="https://github.com/plotly/dash-table.git", rev="master"))
        Pkg.add(PackageSpec(url="https://github.com/plotly/Dash.jl.git", rev="dev"))
    ```

    Once you have installed Dash, you are ready to [make your first Dash app](/getting-started).

    A quick note on checking your versions and on upgrading.
    These docs are run using the versions listed above and these versions should
    be the latest versions available.
    To check which version that you have installed, you can run:

    ```
    Pkg.status("Dash")
    Pkg.status("DashHtmlComponents")
    Pkg.status("DashCoreComponents")
    ```

    To see the latest changes of any package, check the CHANGELOG.MD file in
    the appropriate GitHub repo.

    * [Dash changelog](https://github.com/plotly/Dash.jl/blob/dev/CHANGELOG.md)
    * [DashCoreComponents changelog](https://github.com/plotly/dash-core-components/blob/master/CHANGELOG.md)
    * [DashHTMLComponents changelog](https://github.com/plotly/dash-html-components/blob/master/CHANGELOG.md)
    * [DashTable changelog](https://github.com/plotly/dash-table/blob/master/CHANGELOG.md)
    """),

    html_hr(),
    dcc_markdown("Ready? Now, [let's make your first Dash app.](/getting-started)"),
    dcc_markdown("[Back to the table of contents](/)")

end

end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
function heffmapLapply(A::AbstractTensor{<:Any, N}, L::AbstractTensor{<:Any,3}) where N
    IA = ntuple(i -> ifelse(i == 1, -1, i), N)
    IL = (-1,1,N+1)
    IB = ntuple(identity, N+1)
    B  = tensorcontract(A, IA, L, IL, IB)
    return B
end

function heffmapCapply(A::AbstractTensor{<:Any, N}, C::AbstractTensor{<:Any,4}) where N
    IA = ntuple(i -> ifelse(i == 1, -1, ifelse(i == 2, -2, i-1)), N)
    IC = (-2,-1,1,N)
    IB = ntuple(identity, N)
    B  = tensorcontract(A, IA, C, IC, IB)
    return B
end

function heffmapRapply(A::AbstractTensor{<:Any,N}, R::AbstractTensor{<:Any,3}) where N
    IA = ntuple(i -> ifelse(i == 1, -1, ifelse(i == 2, -2, i-2)), N)
    IR = (-2,-1,N-1)
    IB = ntuple(identity, N-1)
    B  = tensorcontract(A, IA, R, IR, IB)
    return B
end

function heffmapapply(A, heff)
    L, Cs, R = heff[1], heff[2:end-1], heff[end]
    A2 = heffmapLapply(A, L)
    for C in Cs
        A2 = heffmapCapply(A2, C)
    end
    A3 = heffmapRapply(A2, R)
    return A3
end

"""
    heffmap(heff)
returns a function that applies the effective hamiltonian to a tensor.
`heff` is an indexable with entries of type `AbstractTensor` with
`h[1]` being the left-environment or boundary mpo-tensor (rank-3),
`h[end]` being the right-environment or boundary mpo-tensor (rank-3) and
`h[2:end-1]` being the central mpo-tensors (rank-4).
Note that heff needs to have at least 2 elements.
"""
heffmap(heff) = x -> heffmapapply(x, heff)

function heffmapconst(mpo::AbstractMPO{L}, envs, i) where L
    if i == 1
        heff = (mpo[1],envs[1])
    elseif i == L
        heff = (envs[L-1], mpo[L])
    else
        heff = (envs[i-1], mpo[i], envs[i])
    end
    return heffmap(heff)
end


text:
#=
Some internal documentation (especially for people familiar with the nfft)

- The window is precomputed during construction of the NFFT plan
  When performing the nfft convolution, the LUT of the window is used to
  perform linear interpolation. This approach is reasonable fast and does not
  require too much memory. There are, however alternatives known that are either
  faster or require no extra memory at all.

The non-exported functions apodization and convolve are implemented
using Cartesian macros, that may not be very readable.
This is a conscious decision where performance has outweighed readability.
More readable versions can be (and have been) written using the CartesianRange approach,
but at the time of writing this approach require *a lot* of memory.
=#

#=
D is the number of dimensions of the array to be transformed.
DIM is the dimension along which the array is transformed.
DIM == 0 is the ordinary NFFT, i.e., where all dimensions are transformed.
DIM is a type parameter since it allows the @generated macro to
compile more efficient methods.
=#
mutable struct NFFTPlan{D,DIM,T} <: AbstractNFFTPlan{D,DIM,T}
    N::NTuple{D,Int64}
    M::Int64
    x::Matrix{T}
    m::Int64
    sigma::T
    n::NTuple{D,Int64}
    K::Int64
    windowLUT::Vector{Vector{T}}
    windowHatInvLUT::Vector{Vector{T}}
    forwardFFT::FFTW.cFFTWPlan{Complex{T},-1,true,D}
    backwardFFT::FFTW.cFFTWPlan{Complex{T},1,true,D}
    tmpVec::Array{Complex{T},D}
    B::SparseMatrixCSC{T,Int64}
end

function Base.copy(p::NFFTPlan{D,0,T}) where {D,T}
    tmpVec = similar(p.tmpVec)

    FP = plan_fft!(tmpVec; flags=p.forwardFFT.flags)
    BP = plan_bfft!(tmpVec; flags=p.backwardFFT.flags)

    return  NFFTPlan{D,0,T}(p.N, p.M, p.x, p.m, p.sigma, p.n, p.K, p.windowLUT,
                    p.windowHatInvLUT, FP, BP , tmpVec, p.B)
end

function Base.copy(p::NFFTPlan{D,DIM,T}) where {D,DIM,T}
    tmpVec = similar(p.tmpVec)

    FP = plan_fft!(tmpVec, DIM; flags=p.forwardFFT.flags)
    BP = plan_bfft!(tmpVec, DIM; flags=p.backwardFFT.flags)

    return  NFFTPlan{D,DIM,T}(p.N, p.M, p.x, p.m, p.sigma, p.n, p.K, p.windowLUT,
                    p.windowHatInvLUT, FP, BP , tmpVec, p.B)
end

@inline dim(::NFFTPlan{D,DIM}) where {D, DIM} = DIM

##############
# constructors
##############
function NFFTPlan(x::Matrix{T}, N::NTuple{D,Int}, m=4, sigma=2.0,
                window=:kaiser_bessel, K=2000;
                precompute::PrecomputeFlags=LUT, kwargs...) where {D,T}

    if D != size(x,1)
        throw(ArgumentError())
    end

    n = ntuple(d->round(Int,sigma*N[d]), D)

    tmpVec = zeros(Complex{T}, n)

    M = size(x,2)

    FP = plan_fft!(tmpVec; kwargs...)
    BP = plan_bfft!(tmpVec; kwargs...)

    # Create lookup table
    win, win_hat = getWindow(window)

    windowLUT = Vector{Vector{T}}(undef,D)
    windowHatInvLUT = Vector{Vector{T}}(undef,D)
    for d=1:D
        windowHatInvLUT[d] = zeros(T, N[d])
        for k=1:N[d]
            windowHatInvLUT[d][k] = 1. / win_hat(k-1-N[d]/2, n[d], m, sigma)
        end
    end

    if precompute == LUT
        Z = round(Int,3*K/2)
        for d=1:D
            windowLUT[d] = zeros(T, Z)
            for l=1:Z
                y = ((l-1) / (K-1)) * m/n[d]
                windowLUT[d][l] = win(y, n[d], m, sigma)
            end
        end

        B = sparse([],[],Float64[])
    else
        U1 = ntuple(d-> (d==1) ? 1 : (2*m+1)^(d-1), D)
        U2 = ntuple(d-> (d==1) ? 1 : prod(n[1:(d-1)]), D)
        B = precomputeB(win, x, n, m, M, sigma, T, U1, U2)
    end

    NFFTPlan{D,0,T}(N, M, x, m, sigma, n, K, windowLUT, windowHatInvLUT, FP, BP, tmpVec, B )
end

# Directional NFFT
"""
        NFFTPlan(x, d, N, ...) -> plan

Compute *directional* NFFT plan:
A 1D plan that is applied along dimension `d` of a `D` dimensional array of size `N` with sampling locations `x` (a vector).

It takes as optional keywords all the keywords supported by `plan_fft` function (like
`flags` and `timelimit`).  See documentation of `plan_fft` for reference.
"""
function NFFTPlan(x::AbstractVector{T}, dim::Integer, N::NTuple{D,Int64}, m=4,
                       sigma=2.0, window=:kaiser_bessel, K=2000; kwargs...) where {D,T}
    n = ntuple(d->round(Int, sigma*N[d]), D)

    sz = [N...]
    sz[dim] = n[dim]
    tmpVec = Array{Complex{T}}(undef,sz...)

    M = length(x)

    FP = plan_fft!(tmpVec, dim; kwargs...)
    BP = plan_bfft!(tmpVec, dim; kwargs...)

    # Create lookup table
    win, win_hat = getWindow(window)

    windowLUT = Vector{Vector{T}}(undef,1)
    Z = round(Int, 3*K/2)
    windowLUT[1] = zeros(T, Z)
    for l = 1:Z
        y = ((l-1) / (K-1)) * m/n[dim]
        windowLUT[1][l] = win(y, n[dim], m, sigma)
    end

    windowHatInvLUT = Vector{Vector{T}}(undef,1)
    windowHatInvLUT[1] = zeros(T, N[dim])
    for k = 1:N[dim]
        windowHatInvLUT[1][k] = 1. / win_hat(k-1-N[dim]/2, n[dim], m, sigma)
    end

    B = sparse([],[],Float64[])

    NFFTPlan{D,dim,T}(N, M, reshape(x,1,M), m, sigma, n, K, windowLUT,
                      windowHatInvLUT, FP, BP, tmpVec, B)
end

function Base.show(io::IO, p::NFFTPlan{D,0}) where D
    print(io, "NFFTPlan with ", p.M, " sampling points for ", p.N, " array")
end

function Base.show(io::IO, p::NFFTPlan{D,DIM}) where {D,DIM}
    print(io, "NFFTPlan with ", p.M, " sampling points for ", p.N, " array along dimension ", DIM)
end

size(p::NFFTPlan) = p.N
numFourierSamples(p::NFFTPlan) = p.M

################
# nfft functions
################
"""
        nfft!(p, f, fHat) -> fHat

Calculate the NFFT of `f` with plan `p` and store the result in `fHat`.

Both `f` and `fHat` must be complex arrays.
"""
function nfft!(p::NFFTPlan{D,DIM,T}, f::AbstractArray, fHat::StridedArray, verbose=false) where {D,DIM,T}
    consistencyCheck(p, f, fHat)

    fill!(p.tmpVec, zero(Complex{T}))
    t1 = @elapsed @inbounds apodization!(p, f, p.tmpVec)
    if nprocs() == 1
        t2 = @elapsed p.forwardFFT * p.tmpVec # fft!(p.tmpVec) or fft!(p.tmpVec, dim)
    else
        t2 = @elapsed dim(p) == 0 ? fft!(p.tmpVec) : fft!(p.tmpVec, dim(p))
    end
    t3 = @elapsed @inbounds convolve!(p, p.tmpVec, fHat)
    if verbose
      @info "Timing: apod=$t1 fft=$t2 conv=$t3"
    end
    return fHat
end

"""
        nfft_adjoint!(p, fHat, f) -> f

Calculate the adjoint NFFT of `fHat` and store the result in `f`.

Both `f` and `fHat` must be complex arrays.
"""
function nfft_adjoint!(p::NFFTPlan, fHat::AbstractArray, f::StridedArray, verbose=false)
    consistencyCheck(p, f, fHat)

    t1 = @elapsed @inbounds convolve_adjoint!(p, fHat, p.tmpVec)
    if nprocs() == 1
        t2 = @elapsed p.backwardFFT * p.tmpVec # bfft!(p.tmpVec) or bfft!(p.tmpVec, dim)
    else
        t2 = @elapsed dim(p) == 0 ? bfft!(p.tmpVec) : bfft!(p.tmpVec, dim(p))
    end
    t3 = @elapsed @inbounds apodization_adjoint!(p, p.tmpVec, f)
    if verbose
      @info "Timing: conv=$t1 fft=$t2 apod=$t3"
    end
    return f
end


text:
using Test
using Xpress

println(Xpress.getbanner())
println("Optimizer version: $(Xpress.getversion())")

@testset "$(folder)" for folder in [
    "MathOptInterface",
    "xprs_callbacks",
    "Derivative",
]
    @testset "$(file)" for file in readdir(folder)
        include(joinpath(folder, file))
    end
end

@testset "Xpress tests" begin

    prob = Xpress.XpressProblem()

    @test Xpress.getcontrol(prob, "HEURTHREADS") == 0

    vXpress_major = Int(Xpress.getversion().major)
    file_extension = ifelse(vXpress_major <= 38, ".mps","")
    msg = "Unable to call `Xpress.readprob`:\n\n85 Error: File not found: $(file_extension).\n"
    @test_throws Xpress.XpressError(32, msg) Xpress.readprob(prob,"","")
end


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

mutable struct IoK8sApiExtensionsV1beta1IngressSpec <: SwaggerModel
    backend::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1IngressBackend } # spec name: backend
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1IngressRule} } # spec name: rules
    tls::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1IngressTLS} } # spec name: tls

    function IoK8sApiExtensionsV1beta1IngressSpec(;backend=nothing, rules=nothing, tls=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1IngressSpec, Symbol("backend"), backend)
        setfield!(o, Symbol("backend"), backend)
        validate_property(IoK8sApiExtensionsV1beta1IngressSpec, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        validate_property(IoK8sApiExtensionsV1beta1IngressSpec, Symbol("tls"), tls)
        setfield!(o, Symbol("tls"), tls)
        o
    end
end # type IoK8sApiExtensionsV1beta1IngressSpec

const _property_map_IoK8sApiExtensionsV1beta1IngressSpec = Dict{Symbol,Symbol}(Symbol("backend")=>Symbol("backend"), Symbol("rules")=>Symbol("rules"), Symbol("tls")=>Symbol("tls"))
const _property_types_IoK8sApiExtensionsV1beta1IngressSpec = Dict{Symbol,String}(Symbol("backend")=>"IoK8sApiExtensionsV1beta1IngressBackend", Symbol("rules")=>"Vector{IoK8sApiExtensionsV1beta1IngressRule}", Symbol("tls")=>"Vector{IoK8sApiExtensionsV1beta1IngressTLS}")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1IngressSpec }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1IngressSpec))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1IngressSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiExtensionsV1beta1IngressSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1IngressSpec }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1IngressSpec[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1IngressSpec)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1IngressSpec }, name::Symbol, val)
end


text:
using Base.Test, GPUArrays
using GPUArrays: mapidx, gpu_sub2ind
using GPUArrays.TestSuite

function cartesian_iter(state, A, res, Asize)
    for i in CartesianRange(CartesianIndex(Asize))
        idx = gpu_sub2ind(Asize, i.I)
        res[idx] = A[idx]
    end
    return
end

function clmap!(state, f, out, b)
    i = linear_index(state) # get the kernel index it gets scheduled on
    out[i] = f(b[i])
    return
end

function ntuple_test(state, result, ::Val{N}) where N
    result[1] = ntuple(Val{N}) do i
        Float32(i) * 77f0
    end
    return
end

function ntuple_closure(state, result, ::Val{N}, testval) where N
    result[1] = ntuple(Val{N}) do i
        Float32(i) * testval
    end
    return
end

function run_base(Typ)
    @testset "base functionality" begin
        @testset "mapidx" begin
            a = rand(Complex64, 77)
            b = rand(Complex64, 77)
            A = Typ(a)
            B = Typ(b)
            off = UInt32(1)
            mapidx(A, (B, off, UInt32(length(A)))) do i, a, b, off, len
                x = b[i]
                x2 = b[min(i+off, len)]
                a[i] = x * x2
            end
            foreach(1:length(a)) do i
                x = b[i]
                x2 = b[min(i+off, length(a))]
                a[i] = x * x2
            end
            @test Array(A) ≈ a
        end

        @testset "copy!" begin
            x = zeros(Float32, 10, 10)
            y = rand(Float32, 20, 10)
            a = Typ(x)
            b = Typ(y)
            r1 = CartesianRange(CartesianIndex(1, 3), CartesianIndex(7, 8))
            r2 = CartesianRange(CartesianIndex(4, 3), CartesianIndex(10, 8))
            copy!(x, r1, y, r2)
            copy!(a, r1, b, r2)
            @test x == Array(a)

            x2 = zeros(Float32, 10, 10)
            copy!(x2, r1, b, r2)
            @test x2 == x

            fill!(a, 0f0)
            copy!(a, r1, y, r2)
            @test Array(a) == x
        end
        GPUArrays.allowslow(true)
        # right now in CLArrays we fallback to geindex since on some hardware
        # somehow the vcat kernel segfaults -.-
        @testset "vcat + hcat" begin
            x = zeros(Float32, 10, 10)
            y = rand(Float32, 20, 10)
            a = Typ(x)
            b = Typ(y)
            @test vcat(x, y) == Array(vcat(a, b))
            z = rand(Float32, 10, 10)
            c = Typ(z)
            @test hcat(x, z) == Array(hcat(a, c))

            against_base(hcat, Typ{Float32}, (3, 3), (3, 3))
            against_base(vcat, Typ{Float32}, (3, 3), (3, 3))
        end
        GPUArrays.allowslow(false)

        @testset "reinterpret" begin
            a = rand(Complex64, 22)
            A = Typ(a)
            af0 = reinterpret(Float32, a)
            Af0 = reinterpret(Float32, A)
            @test Array(Af0) == af0

            a = rand(Complex64, 10 * 10)
            A = Typ(a)
            af0 = reinterpret(Float32, a, (20, 10))
            Af0 = reinterpret(Float32, A, (20, 10))
            @test Array(Af0) == af0
        end

        @testset "ntuple test" begin
            result = Typ(Vector{NTuple{3, Float32}}(1))
            gpu_call(ntuple_test, result, (result, Val{3}()))
            @test Array(result)[1] == (77, 2*77, 3*77)
            x = 88f0
            gpu_call(ntuple_closure, result, (result, Val{3}(), x))
            @test Array(result)[1] == (x, 2*x, 3*x)
        end

        @testset "cartesian iteration" begin
            Ac = rand(Float32, 32, 32)
            A = Typ(Ac)
            result = zeros(A)
            gpu_call(cartesian_iter, result, (A, result, size(A)))
            Array(result) == Ac
        end

        @testset "Custom kernel from Julia function" begin
            x = Typ(rand(Float32, 100))
            y = Typ(rand(Float32, 100))
            gpu_call(clmap!, x, (-, x, y))
            jy = Array(y)
            @test map!(-, jy, jy) ≈ Array(x)
        end
        T = Typ{Float32}
        @testset "map" begin
            against_base((a, b)-> map(+, a, b), T, (10,), (10,))
            against_base((a, b)-> map!(-, a, b), T, (10,), (10,))
            against_base((a, b, c, d)-> map!(*, a, b, c, d), T, (10,), (10,), (10,), (10,))
        end

        @testset "repmat" begin
            against_base(a-> repmat(a, 5, 6), T, (10,))
            against_base(a-> repmat(a, 5), T, (10,))
            against_base(a-> repmat(a, 5), T, (5, 4))
            against_base(a-> repmat(a, 4, 3), T, (10, 15))
        end
    end
end


text:
abstract type ExpAlgorithm end

struct Exact <: ExpAlgorithm end

struct Trotter{Order} <: ExpAlgorithm
  nsteps::Int
end
one(::Trotter{Order}) where {Order} = Trotter{Order}(1)

function exp(o::∑; alg::ExpAlgorithm=Exact())
  return exp(alg, o)
end

function exp(::Exact, o::∑)
  return ∏([Applied(exp, o)])
end

function exp_one_step(trotter::Trotter{1}, o::∑)
  # TODO: Customize broadcast of `∏`.
  exp_o = ∏(exp.(o))
  return exp_o
end

function exp_one_step(trotter::Trotter{N}, o::∑) where {N}
  exp_o_order_1 = exp_one_step(Trotter{Int(N / 2)}(1), o / 2)
  exp_o = exp_o_order_1 * reverse(exp_o_order_1)
  return exp_o
end

function exp(trotter::Trotter, o::∑)
  expδo = exp_one_step(one(trotter), o / trotter.nsteps)
  return expδo^trotter.nsteps
end


text:
int_rules_4_5_4_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.5.4.1*(a+b*sec)^m*(A+B*sec+C*sec^2) =#
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (1 / (~b) ^ 2) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp((~b * ~B - ~a * ~C) + ~b * ~C * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && EqQ((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C, 0)
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (~C / (~b) ^ 2) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp(-(~a) + ~b * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m], ~x) && EqQ(~A * (~b) ^ 2 + (~a) ^ 2 * ~C, 0)
    @apply_utils Antiderivative((~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~A + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (~A * cot(~e + ~f * ~x) * (~b * csc(~e + ~f * ~x)) ^ ~m) / (~f * ~m) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m], ~x) && EqQ(~C * ~m + ~A * (~m + 1), 0)
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) ^ ~(m') * (~A + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => Antiderivative((~C + ~A * sin(~e + ~f * ~x) ^ 2) / sin(~e + ~f * ~x) ^ (~m + 2), ~x) <-- FreeQ([~e, ~f, ~A, ~C], ~x) && (NeQ(~C * ~m + ~A * (~m + 1), 0) && ILtQ((~m + 1) / 2, 0))
    @apply_utils Antiderivative((~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~A + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (~A * cot(~e + ~f * ~x) * (~b * csc(~e + ~f * ~x)) ^ ~m) / (~f * ~m) + ((~C * ~m + ~A * (~m + 1)) / ((~b) ^ 2 * ~m)) * Antiderivative((~b * csc(~e + ~f * ~x)) ^ (~m + 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C], ~x) && (NeQ(~C * ~m + ~A * (~m + 1), 0) && LeQ(~m, -1))
    @apply_utils Antiderivative((~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~A + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cot(~e + ~f * ~x) * (~b * csc(~e + ~f * ~x)) ^ ~m) / (~f * (~m + 1)) + ((~C * ~m + ~A * (~m + 1)) / (~m + 1)) * Antiderivative((~b * csc(~e + ~f * ~x)) ^ ~m, ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m], ~x) && (NeQ(~C * ~m + ~A * (~m + 1), 0) && Not(LeQ(~m, -1)))
    @apply_utils Antiderivative((~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (~B / ~b) * Antiderivative((~b * csc(~e + ~f * ~x)) ^ (~m + 1), ~x) + Antiderivative((~b * csc(~e + ~f * ~x)) ^ ~m * (~A + ~C * csc(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x)
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~b) * ~C * csc(~e + ~f * ~x) * cot(~e + ~f * ~x)) / (2 * ~f) + (1 / 2) * Antiderivative(Simp(2 * ~A * ~a + (2 * ~B * ~a + ~b * (2 * ~A + ~C)) * csc(~e + ~f * ~x) + 2 * (~a * ~C + ~B * ~b) * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x)
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~b) * ~C * csc(~e + ~f * ~x) * cot(~e + ~f * ~x)) / (2 * ~f) + (1 / 2) * Antiderivative(Simp(2 * ~A * ~a + ~b * (2 * ~A + ~C) * csc(~e + ~f * ~x) + 2 * ~a * ~C * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x)
    @apply_utils Antiderivative((~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2) / (~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => (~C / ~b) * Antiderivative(csc(~e + ~f * ~x), ~x) + (1 / ~b) * Antiderivative((~A * ~b + (~b * ~B - ~a * ~C) * csc(~e + ~f * ~x)) / (~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x)
    @apply_utils Antiderivative((~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2) / (~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => (~C / ~b) * Antiderivative(csc(~e + ~f * ~x), ~x) + (1 / ~b) * Antiderivative((~A * ~b - ~a * ~C * csc(~e + ~f * ~x)) / (~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x)
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(((~a * ~A - ~b * ~B) + ~a * ~C)) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~a * ~f * (2 * ~m + 1)) + (1 / (~a * ~b * (2 * ~m + 1))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp(~A * ~b * (2 * ~m + 1) + (~b * ~B * (~m + 1) - ~a * (~A * (~m + 1) - ~C * ~m)) * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && LtQ(~m, -1 / 2))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~a) * (~A + ~C) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~a * ~f * (2 * ~m + 1)) + (1 / (~a * ~b * (2 * ~m + 1))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp(~A * ~b * (2 * ~m + 1) - ~a * (~A * (~m + 1) - ~C * ~m) * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && LtQ(~m, -1 / 2))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~f * (~m + 1)) + (1 / (~b * (~m + 1))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ ~m * Simp(~A * ~b * (~m + 1) + (~a * ~C * ~m + ~b * ~B * (~m + 1)) * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(LtQ(~m, -1 / 2)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~f * (~m + 1)) + (1 / (~b * (~m + 1))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ ~m * Simp(~A * ~b * (~m + 1) + ~a * ~C * ~m * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(LtQ(~m, -1 / 2)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~f * (~m + 1)) + (1 / (~m + 1)) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m - 1) * Simp(~a * ~A * (~m + 1) + ((~A * ~b + ~a * ~B) * (~m + 1) + ~b * ~C * ~m) * csc(~e + ~f * ~x) + (~b * ~B * (~m + 1) + ~a * ~C * ~m) * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(2 * ~m, 0))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~f * (~m + 1)) + (1 / (~m + 1)) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m - 1) * Simp(~a * ~A * (~m + 1) + (~A * ~b * (~m + 1) + ~b * ~C * ~m) * csc(~e + ~f * ~x) + ~a * ~C * ~m * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(2 * ~m, 0))
    @apply_utils Antiderivative((~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2) / sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => Antiderivative((~A + (~B - ~C) * csc(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ~C * Antiderivative((csc(~e + ~f * ~x) * (1 + csc(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2) / sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => Antiderivative((~A - ~C * csc(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ~C * Antiderivative((csc(~e + ~f * ~x) * (1 + csc(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~f * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / (~a * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp((~A * ((~a) ^ 2 - (~b) ^ 2) * (~m + 1) - ~a * ((~A * ~b - ~a * ~B) + ~b * ~C) * (~m + 1) * csc(~e + ~f * ~x)) + ((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C) * (~m + 2) * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~A * (~b) ^ 2 + (~a) ^ 2 * ~C) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~f * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / (~a * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp((~A * ((~a) ^ 2 - (~b) ^ 2) * (~m + 1) - ~a * ~b * (~A + ~C) * (~m + 1) * csc(~e + ~f * ~x)) + (~A * (~b) ^ 2 + (~a) ^ 2 * ~C) * (~m + 2) * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(2 * ~m) && LtQ(~m, -1)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (1 / ~b) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~A * ~b + (~b * ~B - ~a * ~C) * csc(~e + ~f * ~x)), ~x) + (~C / ~b) * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (1 / ~b) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~A * ~b - ~a * ~C * csc(~e + ~f * ~x)), ~x) + (~C / ~b) * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~(b') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * sec(~(e') + ~(f') * ~x) + ~(C') * sec(~(e') + ~(f') * ~x) ^ 2), ~x) => (~b) ^ 2 * Antiderivative((~b * cos(~e + ~f * ~x)) ^ (~m - 2) * (~C + ~B * cos(~e + ~f * ~x) + ~A * cos(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (~b) ^ 2 * Antiderivative((~b * sin(~e + ~f * ~x)) ^ (~m - 2) * (~C + ~B * sin(~e + ~f * ~x) + ~A * sin(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(b') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(C') * sec(~(e') + ~(f') * ~x) ^ 2), ~x) => (~b) ^ 2 * Antiderivative((~b * cos(~e + ~f * ~x)) ^ (~m - 2) * (~C + ~A * cos(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (~b) ^ 2 * Antiderivative((~b * sin(~e + ~f * ~x)) ^ (~m - 2) * (~C + ~A * sin(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(a') * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~p) ^ ~m * (~(A') + ~(B') * sec(~(e') + ~(f') * ~x) + ~(C') * sec(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~a) ^ IntPart(~m) * (~a * (~b * sec(~e + ~f * ~x)) ^ ~p) ^ FracPart(~m)) / (~b * sec(~e + ~f * ~x)) ^ (~p * FracPart(~m))) * Antiderivative((~b * sec(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~B * sec(~e + ~f * ~x) + ~C * sec(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(a') * (~(b') * csc(~(e') + ~(f') * ~x)) ^ ~p) ^ ~m * (~(A') + ~(B') * csc(~(e') + ~(f') * ~x) + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~a) ^ IntPart(~m) * (~a * (~b * csc(~e + ~f * ~x)) ^ ~p) ^ FracPart(~m)) / (~b * csc(~e + ~f * ~x)) ^ (~p * FracPart(~m))) * Antiderivative((~b * csc(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~B * csc(~e + ~f * ~x) + ~C * csc(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(a') * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~p) ^ ~m * (~(A') + ~(C') * sec(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~a) ^ IntPart(~m) * (~a * (~b * sec(~e + ~f * ~x)) ^ ~p) ^ FracPart(~m)) / (~b * sec(~e + ~f * ~x)) ^ (~p * FracPart(~m))) * Antiderivative((~b * sec(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~C * sec(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(a') * (~(b') * csc(~(e') + ~(f') * ~x)) ^ ~p) ^ ~m * (~(A') + ~(C') * csc(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~a) ^ IntPart(~m) * (~a * (~b * csc(~e + ~f * ~x)) ^ ~p) ^ FracPart(~m)) / (~b * csc(~e + ~f * ~x)) ^ (~p * FracPart(~m))) * Antiderivative((~b * csc(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~C * csc(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
end

max_stars_repo_path: src/rules/4 Trig functions/4.5 Secant/4.5.4.1 (a+b sec)^m (A+B sec+C sec^2).jl


text:
# implementation/kernels.jl
#
# Implements the microkernels for solving the subproblems of the various problems.

@generated function add_micro!{N}(α, A::StridedData{N}, β, C::StridedData{N}, dims::NTuple{N, Int}, offsetA::Int, offsetC::Int)
    quote
        startA = A.start+offsetA
        stridesA = A.strides
        startC = C.start+offsetC
        stridesC = C.strides
        @stridedloops($N, dims, indA, startA, stridesA, indC, startC, stridesC, @inbounds C[indC]=axpby(α,A[indA],β,C[indC]))
        return C
    end
end

@generated function trace_micro!{N}(α, A::StridedData{N}, β, C::StridedData{N}, dims::NTuple{N, Int}, offsetA::Int, offsetC::Int)
    quote
        _scale!(C, β, dims, offsetC)
        startA = A.start+offsetA
        stridesA = A.strides
        startC = C.start+offsetC
        stridesC = C.strides
        @stridedloops($N, dims, indA, startA, stridesA, indC, startC, stridesC, @inbounds C[indC]=axpby(α,A[indA],_one,C[indC]))
        return C
    end
end

@generated function contract_micro!{N}(α, A::StridedData{N}, B::StridedData{N}, β, C::StridedData{N}, dims::NTuple{N, Int}, offsetA, offsetB, offsetC)
    quote
        _scale!(C, β, dims, offsetC)
        startA = A.start+offsetA
        stridesA = A.strides
        startB = B.start+offsetB
        stridesB = B.strides
        startC = C.start+offsetC
        stridesC = C.strides

        @stridedloops($N, dims, indA, startA, stridesA, indB, startB, stridesB, indC, startC, stridesC, @inbounds C[indC] = axpby(α, A[indA]*B[indB], _one, C[indC]))
        return C
    end
end


text:
function compute_temperature_gradient!(dT,dω,V,RTMMchol,T)
  b(v) = ∫(-1.0*DIV(v)*T)dω
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(dT), V)
  solve!(get_free_dof_values(dT),RTMMchol,get_free_dof_values(dT))
end

function compute_buoyancy_flux!(eF,dΩ,V,RTMMchol,e,F)
  b(v) = ∫((v⋅F)*e)dΩ
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(eF), V)
  solve!(get_free_dof_values(eF),RTMMchol,get_free_dof_values(eF))
end

# assume that H1chol factorization has already been performed
function compute_buoyancy!(e,dΩ,S,H1chol,E)
  b(s) = ∫(s*E)dΩ
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(e), S)
  solve!(get_free_dof_values(e),H1chol,get_free_dof_values(e))
end

function compute_velocity_tswe!(u1,dΩ,dω,V,RTMMchol,u2,qAPVM,eAPVM,F,ϕ,dT,n,dt1,dt2)
  b(v) = ∫(v⋅u2 - dt1*(qAPVM)*(v⋅⟂(F,n)) - dt1*(eAPVM)*v⋅dT)dΩ + ∫(dt2*DIV(v)*ϕ)dω
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(u1), V)
  solve!(get_free_dof_values(u1),RTMMchol,get_free_dof_values(u1))
end

function thermal_shallow_water_explicit_time_step!(
     h₂, u₂, E₂, hₚ, uₚ, Eₚ, ϕ, F, q₁, q₂, e₁, e₂, H1h, H1hchol, dT,  # in/out args
     model, dΩ, dω, V, Q, R, S, f, h₁, u₁, E₁, hₘ, uₘ, Eₘ,            # in args
     RTMMchol, L2MMchol, dt, τ, leap_frog,
     assem=SparseMatrixAssembler(SparseMatrixCSC{Float64,Int},Vector{Float64},R,S))                            # more in args

  # energetically balanced explicit second order thermal shallow water solver.
  # extends the explicit shallow water solver with the an additional buoyancy
  # modulated forcing term, and flux form equation for the density weighted
  # buoyancy

  n = get_normal_vector(Triangulation(model))
  # explicit step for provisional velocity, uₚ
  dt1 = dt
  if leap_frog
    dt1 = 2.0*dt
  end

  # 1.1: the mass flux
  compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*h₁)
  # 1.2: the bernoulli function
  compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,u₁⋅u₁,E₁,0.5)
  # 1.3: materially advected quantities (potential vorticity and buoyancy)
  compute_potential_vorticity!(q₁,H1h,H1hchol,dΩ,R,S,h₁,u₁,f,n,assem)
  compute_buoyancy!(e₁,dΩ,S,H1hchol,E₁)
  # 1.4: solve for the provisional velocity
  compute_temperature_gradient!(dT,dω,V,RTMMchol,0.5*h₁)
  compute_velocity_tswe!(uₚ,dΩ,dω,V,RTMMchol,uₘ,q₁-τ*u₁⋅∇(q₁),e₁-τ*u₁⋅∇(e₁),F,ϕ,dT,n,dt1,dt1)
  # 1.5: solve for the provisional depth
  compute_depth!(hₚ,dΩ,dω,Q,L2MMchol,hₘ,F,dt1)
  # 1.6: solve for the buoyancy weighted mass flux
  compute_buoyancy_flux!(dT,dΩ,V,RTMMchol,e₁-τ*u₁⋅∇(e₁),F)
  compute_depth!(Eₚ,dΩ,dω,Q,L2MMchol,Eₘ,dT,dt1)

  # 2.1: the mass flux
  compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*(2.0*h₁ + hₚ)/6.0+uₚ*(h₁ + 2.0*hₚ)/6.0)
  # 2.2: the bernoulli function
  compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,(u₁⋅u₁ + u₁⋅uₚ + uₚ⋅uₚ)/3.0,0.5*(E₁ + Eₚ),0.5)
  # 2.3: materially advected quantities (potential vorticity and buoyancy)
  compute_potential_vorticity!(q₂,H1h,H1hchol,dΩ,R,S,hₚ,uₚ,f,n,assem)
  compute_buoyancy!(e₂,dΩ,S,H1hchol,Eₚ)
  # 2.4: solve for the final velocity
  compute_temperature_gradient!(dT,dω,V,RTMMchol,0.25*(h₁+hₚ))
  compute_velocity_tswe!(u₂,dΩ,dω,V,RTMMchol,u₁,q₁-τ*u₁⋅∇(q₁)+q₂-τ*uₚ⋅∇(q₂),e₁-τ*u₁⋅∇(e₁)+e₂-τ*uₚ⋅∇(e₂),F,ϕ,dT,n,0.5*dt,dt)
  # 2.5: solve for the final depth
  compute_depth!(h₂,dΩ,dω,Q,L2MMchol,h₁,F,dt)
  # 2.6: solve for the buoyancy weighted mass flux
  compute_buoyancy_flux!(dT,dΩ,V,RTMMchol,0.5*(e₁-τ*u₁⋅∇(e₁)+e₂-τ*uₚ⋅∇(e₂)),F)
  compute_depth!(E₂,dΩ,dω,Q,L2MMchol,E₁,dT,dt)
end

function thermal_shallow_water_explicit_time_stepper(
    model, order, degree,
    h₀, u₀, E₀, f₀,
    dt, τ, N;
    mass_matrix_solver::Gridap.Algebra.LinearSolver=Gridap.Algebra.BackslashSolver(),
    write_diagnostics=true,
    write_diagnostics_freq=1,
    dump_diagnostics_on_screen=true,
    write_solution=false,
    write_solution_freq=N/10,
    output_dir="tswe_ncells_$(num_cells(model))_order_$(order)_explicit")

  # Forward integration of the shallow water equations
  Ω = Triangulation(model)
  dΩ = Measure(Ω, degree)
  dω = Measure(Ω, degree, ReferenceDomain())

  # Setup the trial and test spaces
  R, S, U, V, P, Q = setup_mixed_spaces(model, order)

  # assemble the mass matrices
  H1MM, _, L2MM, H1MMchol, RTMMchol, L2MMchol =
    setup_and_factorize_mass_matrices(dΩ, R, S, U, V, P, Q, mass_matrix_solver=mass_matrix_solver)

  # Project the initial conditions onto the trial spaces
  hn, un, f, hnv, unv, fv =  project_shallow_water_initial_conditions(dΩ, Q, V, S,
                               L2MMchol, RTMMchol, H1MMchol, h₀, u₀, f₀)

  b₄(q)   = ∫(q*E₀)dΩ
  rhs4    = assemble_vector(b₄, Q)
  En      = FEFunction(Q, copy(rhs4))
  solve!(get_free_dof_values(En), L2MMchol, get_free_dof_values(En))

  # work arrays
  h_tmp = copy(get_free_dof_values(hn))
  w_tmp = copy(get_free_dof_values(f))
  # build the potential vorticity lhs operator once just to initialise
  bmm(a,b) = ∫(a*hn*b)dΩ
  H1h      = assemble_matrix(bmm, R, S)
  H1hchol  = numerical_setup(symbolic_setup(mass_matrix_solver,H1h),H1h)

  function run_simulation(pvd=nothing)
    diagnostics_file = joinpath(output_dir,"tswe_diagnostics.csv")

    hm1    = clone_fe_function(Q,hn)
    hm2    = clone_fe_function(Q,hn)
    hp     = clone_fe_function(Q,hn)
    ϕ      = clone_fe_function(Q,hn)
    Em1    = clone_fe_function(Q,En)
    Em2    = clone_fe_function(Q,En)
    Ep     = clone_fe_function(Q,En)

    um1    = clone_fe_function(V,un)
    um2    = clone_fe_function(V,un)
    up     = clone_fe_function(V,un)
    F      = clone_fe_function(V,un)
    eF     = clone_fe_function(V,un)

    wn     = clone_fe_function(S,f)
    q1     = clone_fe_function(S,f)
    q2     = clone_fe_function(S,f)
    e1     = clone_fe_function(S,f)
    e2     = clone_fe_function(S,f)

    # first step, no leap frog integration
    thermal_shallow_water_explicit_time_step!(hn, un, En, hp, up, Ep, ϕ, F, q1, q2, e1, e2, H1h, H1hchol, eF,
                                              model, dΩ, dω, V, Q, R, S, f, hm1, um1, Em1, hm2, um2, Em2,
                                              RTMMchol, L2MMchol, dt, τ, false)

    if (write_diagnostics)
      initialize_csv(diagnostics_file,"time", "mass", "vorticity", "buoyancy", "kinetic", "internal", "power_k2p", "power_k2i")
    end

    if (write_diagnostics && write_diagnostics_freq==1)
      compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(Ω))
      dump_diagnostics_thermal_shallow_water!(h_tmp, w_tmp,
                                              model, dΩ, dω, S, L2MM, H1MM,
                                              hn, un, En, wn, ϕ, F, eF, 1, dt,
                                              diagnostics_file,
                                              dump_diagnostics_on_screen)
    end

    # subsequent steps, do leap frog integration
    # (now that we have the state at two previous time levels)
    for istep in 2:N
      h_aux = hm2
      hm2   = hm1
      hm1   = hn
      hn    = h_aux
      u_aux = um2
      um2   = um1
      um1   = un
      un    = u_aux
      E_aux = Em2
      Em2   = Em1
      Em1   = En
      En    = E_aux

      thermal_shallow_water_explicit_time_step!(hn, un, En, hp, up, Ep, ϕ, F, q1, q2, e1, e2, H1h, H1hchol, eF,
                                                model, dΩ, dω, V, Q, R, S, f, hm1, um1, Em1, hm2, um2, Em2,
                                                RTMMchol, L2MMchol, dt, τ, true)

      if (write_diagnostics && write_diagnostics_freq>0 && mod(istep, write_diagnostics_freq) == 0)
        compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(Ω))
        dump_diagnostics_thermal_shallow_water!(h_tmp, w_tmp,
                                                model, dΩ, dω, S, L2MM, H1MM,
                                                hn, un, En, wn, ϕ, F, eF, istep, dt,
                                                diagnostics_file,
                                                dump_diagnostics_on_screen)
      end
      if (write_solution && write_solution_freq>0 && mod(istep, write_solution_freq) == 0)
        compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(Ω))
	pvd[dt*Float64(istep)] = new_vtk_step(Ω,joinpath(output_dir,"n=$(istep)"),["hn"=>hn,"un"=>un,"wn"=>wn,"en"=>e2])
      end
    end
    hn, un, En
  end
  if (write_diagnostics || write_solution)
    rm(output_dir,force=true,recursive=true)
    mkdir(output_dir)
  end
  if (write_solution)
    pvdfile=joinpath(output_dir,"tswe_ncells_$(num_cells(model))_order_$(order)_explicit")
    paraview_collection(run_simulation,pvdfile)
  else
    run_simulation()
  end
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
const matrix_directory = Pkg.dir("CSparse", "Matrix")

function testmat(nm)
    fd = open(joinpath(matrix_directory, nm), "r")
    r1 = split(readline(fd))
    rows = [parse(Int32,r1[1])]
    cols = [parse(Int32,r1[2])]
    cmpx = length(r1) == 4
    if !cmpx && length(r1) != 3
        error("$(file_path(TestMats, nm)) is not a CXSparse test matrix")
    end
    nzs = cmpx ? [float(r1[3]) + im*float(r1[4])] : [float(r1[3])]
    for rr = EachLine(fd)
        vv = split(rr)
        push!(rows,parse(Int32,vv[1]))
        push!(cols,parse(Int32,vv[2]))
        push!(nzs, cmpx ? float(vv[3]) + im*float(vv[4]) : float(vv[3]))
    end
    close(fd)
    sparse(rows .+ one(Int32), cols .+ one(Int32), nzs)
end


text:
@testset "Data" begin

@test artificialInSARImage(2) == S1Point.(zeros(2,2))
@test artificialInSARImage(2,RnPoint) == RnPoint.(2*π*ones(2,2))
@test artificialInSARImage(2,SnPoint) == repeat([SnPoint([1.,0.])],2,2)

@test artificialS1SlopeSignal(20,0.) == repeat([S1Point(-π/2)],20)

@test ismissing( artificialS1Signal(-1.) )
@test ismissing( artificialS1Signal(2.) )
@test artificialS1Signal(2) == S1Point.([-3*π/4,-3*π/4])

# for the remainder check data types only
@test length( artificialS1Signal(20) ) == 20

@test size( artificialS2WhirlImage(64) ) == (64,64)
@test length(getValue(artificialS2WhirlImage(64)[1,1])) == 3
@test eltype( getValue(artificialS2WhirlImage(64)) ) == SnPoint{Float64}

@test size( artificialS2RotationsImage(64) ) == (64,64)
@test length(getValue(artificialS2RotationsImage(64)[1,1])) == 3
@test eltype( getValue(artificialS2RotationsImage(64)) ) == SnPoint{Float64}

@test size( artificialS2WhirlPatch(8) ) == (8,8)
@test length(getValue(artificialS2WhirlPatch(8)[1,1])) == 3
@test eltype( getValue(artificialS2WhirlPatch(8)) ) == SnPoint{Float64}

@test size( artificialSPDImage(8) ) == (8,8)
@test size(getValue(artificialSPDImage(8)[1,1])) == (3,3)
@test eltype( getValue(artificialSPDImage(8)) ) == SPDPoint{Float64}

@test size( artificialSPDImage2(8) ) == (8,8)
@test size(getValue(artificialSPDImage2(8)[1,1])) == (3,3)
@test eltype( getValue(artificialSPDImage2(8)) ) == SPDPoint{Float64}

@test length( artificialS2Lemniscate(SnPoint([0.,0.,1.]),20) ) == 20
@test length( getValue( artificialS2Lemniscate(SnPoint([0.,0.,1.]),20)[1])) == 3
@test eltype( artificialS2Lemniscate(SnPoint([0.,0.,1.]),20)) == SnPoint{Float64}

end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots, LaTeXStrings
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.19"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	alpha, beta = 8, 2
	prior(lam) = pdf(Gamma(alpha, 1/beta), lam)
	data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

	like(lam) = *([pdf(Poisson(lam),x) for x in data]...)
	posteriorUpToK(lam) = like(lam)*prior(lam)
end

# ╔═╡ 9ed0aafc-3afb-11eb-1f46-5bc8b9a49938
begin
	delta = 10^-4.
	lamRange = 0:delta:10
	K = sum([posteriorUpToK(lam)*delta for lam in lamRange])
	posterior(lam) = posteriorUpToK(lam)/K
end

# ╔═╡ 9ed0def2-3afb-11eb-2bd0-5b8a32cba3e1
bayesEstimate = sum([lam*posterior(lam)*delta for lam in lamRange])

# ╔═╡ 9ed16ca0-3afb-11eb-332d-3ba7b882417f
newAlpha, newBeta = alpha + sum(data), beta + length(data)

# ╔═╡ 9ed86364-3afb-11eb-0bb4-a586cc54fb14
closedFormBayesEstimate = mean(Gamma(newAlpha, 1/newBeta))

# ╔═╡ 9ee4be2c-3afb-11eb-02df-0f03e0f586fc
Text("Computational Bayes Estimate: $(bayesEstimate)")

# ╔═╡ 9ee54e44-3afb-11eb-19ef-511b9b6134a4
Text("Closed form Bayes Estimate: $(closedFormBayesEstimate)")

# ╔═╡ 9ef09b16-3afb-11eb-1cc0-a7d27312610f
begin
	plot(lamRange, prior.(lamRange),
		c=:blue, label="Prior distribution")
	plot!(lamRange, posterior.(lamRange),
		c=:red, label="Posterior distribution",
		xlims=(0, 10), ylims=(0, 1.2),
			xlabel=L"\lambda",ylabel="Density")
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.19"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═9ed0aafc-3afb-11eb-1f46-5bc8b9a49938
# ╠═9ed0def2-3afb-11eb-2bd0-5b8a32cba3e1
# ╠═9ed16ca0-3afb-11eb-332d-3ba7b882417f
# ╠═9ed86364-3afb-11eb-0bb4-a586cc54fb14
# ╠═9ee4be2c-3afb-11eb-02df-0f03e0f586fc
# ╠═9ee54e44-3afb-11eb-19ef-511b9b6134a4
# ╠═9ef09b16-3afb-11eb-1cc0-a7d27312610f
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using NBodyIPs
using JuLIP, Test, StaticArrays

# TODO: * monomials testset
#       * fast_polys testset

@testset "NBodyIPs" begin
   # @testset "Iterators" begin include("test_iterators.jl") end
   @testset "Invariants" begin include("test_invariants.jl") end
   @testset "Cutoffs" begin include("test_cutoffs.jl") end
   @testset "OneBody" begin include("test_onebody.jl") end
   @testset "Polys" begin include("test_polys.jl") end
   @testset "Statxic Polys" begin include("test_stpolys.jl") end
   @testset "NBodyIP IO" begin include("test_io.jl") end
   @testset "EnvironmentIPs" begin include("test_environ.jl") end
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# ------------------------------------------------------------------------------------------
# # Quantifying Chaos for a Dynamical System
#
# Topics:
# * Introduction to Lyapunov exponents
# * Maximum Lyapunov Exponent
# * Lyapunov Spectrum
#
# **WARNING** - Compilation of functions in this tutorial takes *a lot of time*.
#
# ---
#
# In the previous tutorial we saw that for example the Shinriki oscillator showed something
# that "could be chaotic behavior" for a specific parameter value. How can we quantify that?
#
# # Lyapunov exponents
# * Lyapunov exponents measure the exponential separation rate of trajectories that are
# (initially) close.
#     * Consider the following picture, where two nearby trajectories are evolved in time:
#
#
#
# <img src="lyapunov.png" alt="Sketch of the Lyapunov exponent" style="width: 500px;"/>
#
#
# * $\lambda$ denotes the "maximum Lyapunov exponent".
# * A $D$-dimensional system has $D$ exponents.
# * In general, a trajectory is called "chaotic" if
#     1. it follows nonlinear dynamics
#     2. it is *bounded* (does not escape to infinity)
#     2. it has a positive Lyapunov exponent
#
# *(please be aware that the above is an over-simplification! See the textbooks cited in our
# documentation for more)*
#
# ---
# ------------------------------------------------------------------------------------------

using DynamicalSystems, BenchmarkTools, PyPlot

# ------------------------------------------------------------------------------------------
# Before computing Lyapunov exponents, we'll demonstrate the concept of exponential
# separation using a simple map: the *towel map*
#
# $$
# \begin{aligned}
# x_{n+1} &= a x_n (1-x_n) -0.05 (y_n +0.35) (1-2z_n) \\
# y_{n+1} &= 0.1 \left( \left( y_n +0.35 \right)\left( 1+2z_n\right) -1 \right)
# \left( 1 -1.9 x_n \right) \\
# z_{n+1} &= 3.78 z_n (1-z_n) + b y_n
# \end{aligned}
# $$
# ------------------------------------------------------------------------------------------

towel = Systems.towel()

# ------------------------------------------------------------------------------------------
# First we'll generate a trajectory for the towel map, `tr1`, from the default initial
# condition,
# ------------------------------------------------------------------------------------------

tr1 = trajectory(towel, 100)
summary(tr1)

# ------------------------------------------------------------------------------------------
# and then we will generate a second trajectory, `tr2`, with a starting point slightly
# shifted from the initial condition of `tr1`.
# ------------------------------------------------------------------------------------------

u2 = get_state(towel) + (1e-9 * rand(3))
tr2 = trajectory(towel, 100, u2)
summary(tr2)

figure(figsize=(8,5))

# Plot the x-coordinate of the two trajectories:
ax1 = subplot(2,1,1)
plot(tr1[:, 1], alpha = 0.5)
plot(tr2[:, 1], alpha = 0.5)
ylabel("x")

# Plot their distance in a semilog plot:
ax2 = subplot(2,1,2, sharex = ax1)
d = [norm(tr1[i] - tr2[i]) for i in 1:length(tr2)]
ylabel("d")
xlabel("n")
semilogy(d);

# ------------------------------------------------------------------------------------------
# # Maximum Lyapunov Exponent
# `lyapunov` is a function that calculates the maximum Lyapunov exponent for a
# `DynamicalSystem` (for a given starting point).
#
# Since `lyapunov` is not a trivial function, it is best to read the documentation string
# first:
# ------------------------------------------------------------------------------------------

?lyapunov

# ------------------------------------------------------------------------------------------
# ---
#
# Let's apply this to the example of the previous section, the Shinriki oscillator!
#
# *Reminder:* we found something that "could" be chaotic behavior for the parameter `R1 =
# 21.0`
#
# ------------------------------------------------------------------------------------------

shi = Systems.shinriki(;R1 = 21.0)

lyapunov(shi, 1000.0, Ttr = 10.0)

# ------------------------------------------------------------------------------------------
# Positive Lyapunov exponent!?!? That is definitely chaotic behavior, right?
#
# *Right?*
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's increase accuracy of computation, as well as the transient time, to see whether an
# orbit actually enters a limit cycle.
# ------------------------------------------------------------------------------------------

lyapunov(shi, 2000.0, Ttr = 1000.0)

# ------------------------------------------------------------------------------------------
# **surprise**
#
# Let's find out whats going on!
#
# To see why the lyapunov exponent was positive in one case and negative in another, we can
# produce a more detailed orbit diagram, around the "critical" value of `R = 21.0`.
# ------------------------------------------------------------------------------------------

pvalues = linspace(20.9,21.1,101)
i = 1
plane = (2, 0.0)
tf = 1000.0
p_index = 1

# use extremely long transient time:
output = produce_orbitdiagram(shi, plane, i, p_index, pvalues; tfinal = tf,
                              Ttr = 2000.0, direction = -1, printparams = false);
length(output)

using PyPlot
figure(figsize=(6,4))
for (j, p) in enumerate(pvalues)
    plot(p .* ones(output[j]), output[j], lw = 0,
    marker = "o", ms = 0.2, color = "black")
end
plot([21, 21], [-2.1, 0.1], color = "red", alpha = 0.55)
xlabel("\$R_1\$"); ylabel("\$V_1\$");

# result of orbit diagram at R1 = 21.0
values = output[51]

# Amount of unique points
un = unique(round.(output[51], 8))
println("Total: $(length(output[51])), unique: $(length(un))")

# ------------------------------------------------------------------------------------------
# ## `lyapunov` for discrete system
#
# * All functions that accept a `DynamicalSystem` work with *any* instance of
# `DynamicalSystem`, regardless of whether it is continuous, discrete, in-place, out-of-
# place or whatever.
#
# ------------------------------------------------------------------------------------------

# Get the Henon map from the library of pre-defined systems:
hen = Systems.henon()
λ = lyapunov(hen, 10000)

# ------------------------------------------------------------------------------------------
# ---
#
# # Lyapunov Spectrum
#
# Besides the maximum Laypunov exponent, the function `lyapunovs` (with `s` at the end)
# returns the entire Lyapunov spectrum (or as many exponents the user wants).
#
#
#
# Once again, because the function `lyapunovs` is not trivial, we will view the
# documentation string first:
# ------------------------------------------------------------------------------------------

?lyapunovs

# ------------------------------------------------------------------------------------------
# ### Lyapunovs for discrete systems
#
# In our first example of calling `lyapunovs`, let's pass a discrete system.
# ------------------------------------------------------------------------------------------

towel = Systems.towel()
lyapunovs(towel, 2000; Ttr = 200)

# ------------------------------------------------------------------------------------------
# Here we're choosing to compute only the first two exponents.
# ------------------------------------------------------------------------------------------

lyapunovs(towel, 2000, 2; Ttr = 200)

# ------------------------------------------------------------------------------------------
# If you only want the first exponent (maximum), use the `lyapunov` function instead
# ------------------------------------------------------------------------------------------

lyapunov(towel, 2000; Ttr = 200)

# ------------------------------------------------------------------------------------------
# How much time does this take?
# ------------------------------------------------------------------------------------------

using BenchmarkTools
@btime lyapunovs($towel, 2000; Ttr = 200);

@btime lyapunov($towel, 2000; Ttr = 200);

# ------------------------------------------------------------------------------------------
# ### Lyapunov exponents for continuous systems
#
# Next, let's initialize the Lorenz system with random initial condition
# ------------------------------------------------------------------------------------------

lor = Systems.lorenz()

# ------------------------------------------------------------------------------------------
# We'll compute the Lyapunov spectrum with specified initial parallepiped matrix `Q0`:
# ------------------------------------------------------------------------------------------

Q0 = eye(3)
lyapunovs(lor, 2000, Q0; Ttr = 10.0)

# ------------------------------------------------------------------------------------------
# And we find that results "converge" already with 2000 iterations:
# ------------------------------------------------------------------------------------------

lyapunovs(lor, 3000, Q0; Ttr = 10.0)

# ------------------------------------------------------------------------------------------
# * Even the continuous systems are quite performant (note that compilation takes a **lot**
# of time):
# ------------------------------------------------------------------------------------------

@btime lyapunovs($lor, 2000, $Q0; Ttr = 10.0);

# ------------------------------------------------------------------------------------------
# * The above integration is done with a 9th order solver and tolerances of `1e-9`. But you
# can get away with lower tolerances.
#
# Let's load in `OrdinaryDiffEq` and specify keyword arguments for the integrators of
# DifferentialEquations.jl.
# ------------------------------------------------------------------------------------------

using OrdinaryDiffEq

dek = Dict(:solver => Tsit5(), :abstol => 1e-6, :reltol => 1e-6)

# ------------------------------------------------------------------------------------------
# And now we can call `lyapunovs` with the keyword `diff_eq_kwargs`
# ------------------------------------------------------------------------------------------

lyapunovs(lor, 2000, Q0; Ttr = 10.0, diff_eq_kwargs = dek)

@btime lyapunovs($lor, 2000.0, $Q0; Ttr = 10.0, diff_eq_kwargs = $dek);

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/introduction-to-dynamicalsystems.jl/3. Quantifying Chaos (Lyapunov Exponents).jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using VDPTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using PmapProgressMeter
using MCTS
using ProfileView

@everywhere begin
    using POMDPs
    using POMDPToolbox
    using MCTS
    using VDPTag

    mdp = VDPTagMDP()
    seed=27

    solver = DPWSolver(depth=20,
                       exploration_constant=40.0,
                       rng=MersenneTwister(seed),
                       estimate_value=RolloutEstimator(ToNextML(mdp)),
                       n_iterations=10000,
                       k_action=8.0,
                       alpha_action=1/20,
                       k_state=4.0,
                       alpha_state=1/20)
end

#=
N = 100
# s_rewards = SharedArray(Float64, N)
prog = Progress(N, desc="Simulating...")
rewards = pmap(prog, 1:N) do i
# @showprogress "Simulating..." for i in 1:N
    policy = solve(solver, mdp)

    hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(i))
    hist = simulate(hr, mdp, policy)
    discounted_reward(hist)
end
@show mean(rewards)
=#

policy = solve(solver, mdp)
a = action(policy, initial_state(mdp, MersenneTwister(1)))
clear_tree!(policy)

Profile.clear()
hr = HistoryRecorder(max_steps=20, rng=MersenneTwister(1))
hist = simulate(hr, mdp, policy)
clear_tree!(policy)
@profile hist = simulate(hr, mdp, policy)
# @profile a = action(policy, initial_state(mdp, MersenneTwister(1)))
ProfileView.view()

#=
hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(1))
hist = simulate(hr, mdp, policy)
@show discounted_reward(hist)

gr()
frames = Frames(MIME("image/png"), fps=2)
@showprogress "Creating gif..." for s in state_hist(hist)
    push!(frames, plot(mdp, s))
end

filename = string(tempname(), "_vdprun.gif")
write(filename, frames)
println(filename)
run(`setsid gifview $filename`)
=#


text:
# This test shows the interest to update the QR basis instead
# of computing from scratch at every time step

# function timing()
#     @btime begin
#         m = 1000
#         n = 50
#
#         A = randn(m, n)
#         qr(A)
#         for i=1:20
#             v = randn(m)
#             A = hcat(A, reshape(v,(m,1)))
#         end
#     end
#     @btime begin
#         m = 1000
#         n = 50
#
#         A = randn(m, n)
#         qr(A)
#         for i=1:20
#             v = randn(m)
#             A = hcat(A, reshape(v,(m,1)))
#             qr(A)
#         end
#     end
#     @btime begin
#                 m = 1000
#         n = 50
#
#         A = randn(m, n)
#         F = qr(A)
#         R = F.R
#         for i=1:20
#             v = randn(m)
#             R = qraddcol(A, R, v)
#
#             A = hcat(A, reshape(v,(m,1)))
#         end
#     end
# end
#
# 1.364 ms (86 allocations: 10.18 MiB)
# 17.929 ms (194 allocations: 20.08 MiB)
# 2.417 ms (368 allocations: 11.47 MiB)


text:
import PotentialUQ
using LinearAlgebra: I
using Random
Random.seed!(1234);

n = 100
β = ones(10)
A = 1e-1*randn(n, 10)
A[1:10, 1:10] += I(10)
b = randn(n)

Q = Matrix(1e-3*I(n))

snap = PotentialUQ.SNAP(A, b, β)
dsnap = PotentialUQ.SNAPDistribution(snap, Q)

@time fitsnap = PotentialUQ.MAP(dsnap)
dsnap = PotentialUQ.SNAPDistribution(snap, Q)
@time fitsnap = PotentialUQ.MAP(dsnap)

println("Solution ", A\b)
println("MAP Solution, Optimization ", fitsnap.x.β)


text:
# [1] Pratt, William K., Digital Image Processing, New York, John Wiley & Sons, Inc., 1991, p. 634.
function measure_feature(property::Measurement, t::IndexedTable, labels::AbstractArray, N::Int = 0)
    fill_properties(property, t)
end

function fill_properties(property::Measurement, t::IndexedTable)
    t = property.area ? compute_area(t) : t
    t = property.perimeter ? compute_perimeter(t) : t
end

function compute_area(t::IndexedTable)
    # Equation 18.2-8a
    # Pratt, William K., Digital Image Processing, New York, John Wiley & Sons, Inc., 1991, p. 629.
    @transform t {area = ((1/4)*:Q₁ + (1/2)*:Q₂ + (7/8)*:Q₃ + :Q₄ +(3/4)*:Qₓ) }
end

function compute_perimeter(t::IndexedTable)
    # perimiter₀ and perimter₁ are given by equations 18.2-8b and 18.2-7a in [1].
    # perimeter₂ is given by equation (32) in [2]
    # [1] Pratt, William K., Digital Image Processing, New York, John Wiley & Sons, Inc., 1991, p. 629.
    # [2] S. B. Gray, “Local Properties of Binary Images in Two Dimensions,” IEEE Transactions on Computers, vol. C–20, no. 5, pp. 551–561, May 1971. https://doi.org/10.1109/t-c.1971.223289
    @transform t {perimeter₀ = (:Q₂ + (1/sqrt(2))*(:Q₁ + :Q₃ + 2*:Qₓ)), perimeter₁ = (:Q₁ + :Q₂ + :Q₃ + 2*:Qₓ), perimeter₂ = (:Q₂ + (:Q₁ +:Q₃)/sqrt(2)) }
end


text:
using Flux, Flux.Tracker
using Distributions

struct VP
  b0
  b1

  VP() = new(param(randn(2)), param(randn(2)))
end

sigmoid(x) = 1.0 / (1.0 + exp(-x))

function loglike(y, x, params, N)
  b0 = params[:b0]
  b1 = params[:b1]
  prob = clamp.(sigmoid.(b0 .+ b1 .* x), .001, .999)
  # return sum(y .* log.(prob) .+ (1 .- y) .* log1p.(-prob)) * (N / length(y))
  return sum(logpdf.(Bernoulli.(prob), y)) * (N / length(y))
end

function log_p(params)
  # return 0.0
  b0 = params[:b0]
  b1 = params[:b1]
  return logpdf(Normal(0, 10), b0) + logpdf(Normal(0, 10), b1)
end

lpdf_vp(p, x) = logpdf(Normal(p[1], exp(p[2])), x)
function log_q(params, vp)
  return lpdf_vp(vp.b0, params[:b0]) + lpdf_vp(vp.b1, params[:b1])
end

rsample(p) = randn() * exp(p[2]) + p[1]
rsample(vp::VP) = Dict(:b0 => rsample(vp.b0), :b1 => rsample(vp.b1))

struct Metrics
  ll::Float64
  lp::Float64
  lq::Float64
end

function elbo(y, x, vp, N, metrics::Vector{Metrics})
  params = rsample(vp)

  ll = loglike(y, x, params, N)
  lp = log_p(params)
  lq = log_q(params, vp)

  m = Metrics(ll.data, lp.data, lq.data)
  push!(metrics, m)

  return ll + lp - lq
end


text:
struct Bottleneck
    layer
end

@functor Bottleneck

Bottleneck(in_planes, growth_rate) = Bottleneck(Chain(BatchNorm(in_planes, relu),
                                          Conv((1, 1), in_planes => 4growth_rate),
                                          BatchNorm(4growth_rate, relu),
                                          Conv((3, 3), 4growth_rate => growth_rate, pad = (1, 1))))

(b::Bottleneck)(x) = cat(b.layer(x), x, dims = 3)

Transition(chs::Pair{<:Int,<:Int}) = Chain(BatchNorm(chs[1], relu),
                                            Conv((1, 1), chs),
                                            MeanPool((2, 2)))

function _make_dense_layers(block, in_planes, growth_rate, nblock)
    local layers = []
    for i in 1:nblock
        push!(layers, block(in_planes, growth_rate))
        in_planes += growth_rate
    end
    Chain(layers...)
end

function _densenet(nblocks = [6, 12, 24, 16]; block = Bottleneck, growth_rate = 32, reduction = 0.5, num_classes = 1000)
    num_planes = 2growth_rate
    layers = []
    push!(layers, Conv((7, 7), 3 => num_planes, stride = (2, 2), pad = (3, 3)))
    push!(layers, BatchNorm(num_planes, relu))
    push!(layers, MaxPool((3, 3), stride = (2, 2), pad = (1, 1)))

    for i in 1:3
        push!(layers, _make_dense_layers(block, num_planes, growth_rate, nblocks[i]))
        num_planes += nblocks[i] * growth_rate
        out_planes = Int(floor(num_planes * reduction))
        push!(layers, Transition(num_planes => out_planes))
        num_planes = out_planes
    end

    push!(layers, _make_dense_layers(block, num_planes, growth_rate, nblocks[4]))
    num_planes += nblocks[4] * growth_rate
    push!(layers, BatchNorm(num_planes, relu))

    Chain(layers..., MeanPool((7, 7)),
        x->reshape(x, :, size(x, 4)),
        Dense(num_planes, num_classes), softmax)
end

function densenet_layers()
    weight = Metalhead.weights("densenet.bson")
    weights = Dict{Any,Any}()
    for ele in keys(weight)
        weights[string(ele)] = convert(Array{Float64,N} where N, weight[ele])
    end
    ls = _densenet()
    ls[1].weight .= weights["conv1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]
    ls[2].β .= weights["conv1/bn_b_0"]
    ls[2].γ .= weights["conv1/bn_w_0"]
    l = 4
    for (c, n) in enumerate([6, 12, 24, 16])
        for i in 1:n
            for j in [2, 4]
                ls[l][i].layer[j].weight .= weights["conv$(c + 1)_$i/x$(j ÷ 2)_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]
                ls[l][i].layer[j - 1].β .= weights["conv$(c + 1)_$i/x$(j ÷ 2)/bn_b_0"]
                ls[l][i].layer[j - 1].γ .= weights["conv$(c + 1)_$i/x$(j ÷ 2)/bn_w_0"]
            end
        end
        l += 2
    end
    for i in [5, 7, 9] # Transition Block Conv Layers
        ls[i][2].weight .= weights["conv$(i ÷ 2)_blk_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:]
        ls[i][1].β .= weights["conv$(i ÷ 2)_blk/bn_b_0"]
        ls[i][1].γ .= weights["conv$(i ÷ 2)_blk/bn_w_0"]
    end
    ls[end - 1].W .= transpose(dropdims(weights["fc6_w_0"], dims = (1, 2))) # Dense Layers
    ls[end - 1].b .= weights["fc6_b_0"]
    return ls
end

struct DenseNet <: ClassificationModel{ImageNet.ImageNet1k}
    layers::Chain
end

DenseNet() = DenseNet(densenet_layers())

Base.show(io::IO, ::DenseNet) = print(io, "DenseNet()")

@functor DenseNet

(m::DenseNet)(x) = m.layers(x)


text:
using Robust32s, LinearAlgebra, Test

streq(a, b) = string(a) === string(b)

include("primitive.jl")
include("rounding.jl")
include("arithmetic.jl")
include("vector.jl")
include("matrix.jl")
include("elementary.jl")


text:
using Luxor, Test, LaTeXStrings, MathTeXEngine

using MathTeXEngine

@testset "latex" begin
    d = Drawing(200,200,:svg)
    origin()
    fontsize(18)
    t = L"\sum^N_{n=1} \frac{x}{n}"

    for halign in [:left,:right,:center], valign in [:baseline, :bottom, :top, :middle]
        text(t, halign=halign, valign=valign)

        @test round.(latextextsize(t);digits=2) == (37.78, 45.98)

        if halign === :left
            @test round(latexboundingbox(t, halign=halign, valign=valign)[1][1];digits=2) == 0.0
            @test round(latexboundingbox(t, halign=halign, valign=valign)[2][1];digits=2) == 37.78
        elseif halign === :right
            @test round(latexboundingbox(t, halign=halign, valign=valign)[1][1];digits=2) == -37.78
            @test round(latexboundingbox(t, halign=halign, valign=valign)[2][1];digits=2) == 0.0
        elseif halign === :center
            @test round(latexboundingbox(t, halign=halign, valign=valign)[1][1];digits=2) == -37.78/2
            @test round(latexboundingbox(t, halign=halign, valign=valign)[2][1];digits=2) == 37.78/2
        end

        if valign === :baseline || valign === :bottom
            @test round(latexboundingbox(t, halign=halign, valign=valign)[1][2];digits=2) == 18.25
            @test round(latexboundingbox(t, halign=halign, valign=valign)[2][2];digits=2) == -27.72
        elseif valign === :top
            @test round(latexboundingbox(t, halign=halign, valign=valign)[1][2];digits=2) == 45.98
            @test round(latexboundingbox(t, halign=halign, valign=valign)[2][2];digits=2) == 0.0
        elseif valign === :middle
            @test round(latexboundingbox(t, halign=halign, valign=valign)[1][2];digits=2) == 45.98/2
            @test round(latexboundingbox(t, halign=halign, valign=valign)[2][2];digits=2) == -45.98/2
        end
    end

    @test round(rawlatexboundingbox(t)[1][1]; digits=2) == 0.0
    @test round(rawlatexboundingbox(t)[1][2]; digits=2) == round(18.25/get_fontsize(); digits=2)
    @test round(rawlatexboundingbox(t)[2][1]; digits=2) == round(37.78/get_fontsize(); digits=2)
    @test round(rawlatexboundingbox(t)[2][2]; digits=2) == round(-27.72/get_fontsize(); digits=2)

    finish()
end