text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using SuiteSparse.UMFPACK
using SuiteSparse
using SuiteSparse: increment!
using Serialization
using LinearAlgebra:
    I, det, issuccess, ldiv!, lu, lu!, Adjoint, Transpose, SingularException, Diagonal
using SparseArrays: nnz, sparse, sprand, sprandn, SparseMatrixCSC

@testset "UMFPACK wrappers" begin
    se33 = sparse(1.0I, 3, 3)
    do33 = fill(1., 3)
    @test isequal(se33 \ do33, do33)

    # based on deps/Suitesparse-4.0.2/UMFPACK/Demo/umfpack_di_demo.c

    A0 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
                increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
                [2.,1.,3.,4.,-1.,-3.,3.,6.,2.,1.,4.,2.], 5, 5)

    @testset "Core functionality for $Tv elements" for Tv in (Float64, ComplexF64)
        # We might be able to support two index sizes one day
        for Ti in Base.uniontypes(SuiteSparse.UMFPACK.UMFITypes)
            A = convert(SparseMatrixCSC{Tv,Ti}, A0)
            lua = lu(A)
            @test nnz(lua) == 18
            @test_throws ErrorException lua.Z
            L,U,p,q,Rs = lua.:(:)
            @test (Diagonal(Rs) * A)[p,q] ≈ L * U

            det(lua) ≈ det(Array(A))

            b = [8., 45., -3., 3., 19.]
            x = lua\b
            @test x ≈ float([1:5;])

            @test A*x ≈ b
            z = complex.(b)
            x = ldiv!(lua, z)
            @test x ≈ float([1:5;])
            @test z === x
            y = similar(z)
            ldiv!(y, lua, complex.(b))
            @test y ≈ x

            @test A*x ≈ b

            b = [8., 20., 13., 6., 17.]
            x = lua'\b
            @test x ≈ float([1:5;])

            @test A'*x ≈ b
            z = complex.(b)
            x = ldiv!(adjoint(lua), z)
            @test x ≈ float([1:5;])
            @test x === z
            y = similar(x)
            ldiv!(y, adjoint(lua), complex.(b))
            @test y ≈ x

            @test A'*x ≈ b
            x = transpose(lua) \ b
            @test x ≈ float([1:5;])

            @test transpose(A) * x ≈ b
            x = ldiv!(transpose(lua), complex.(b))
            @test x ≈ float([1:5;])
            y = similar(x)
            ldiv!(y, transpose(lua), complex.(b))
            @test y ≈ x

            @test transpose(A) * x ≈ b

            # Element promotion and type inference
            @inferred lua\fill(1, size(A, 2))
        end
    end

    @testset "More tests for complex cases" begin
        Ac0 = complex.(A0,A0)
        for Ti in Base.uniontypes(SuiteSparse.UMFPACK.UMFITypes)
            Ac = convert(SparseMatrixCSC{ComplexF64,Ti}, Ac0)
            x  = fill(1.0 + im, size(Ac,1))
            lua = lu(Ac)
            L,U,p,q,Rs = lua.:(:)
            @test (Diagonal(Rs) * Ac)[p,q] ≈ L * U
            b  = Ac*x
            @test Ac\b ≈ x
            b  = Ac'*x
            @test Ac'\b ≈ x
            b  = transpose(Ac)*x
            @test transpose(Ac)\b ≈ x
        end
    end

    @testset "Rectangular cases. elty=$elty, m=$m, n=$n" for
        elty in (Float64, ComplexF64),
            (m, n) in ((10,5), (5, 10))

        Random.seed!(30072018)
        A = sparse([1:min(m,n); rand(1:m, 10)], [1:min(m,n); rand(1:n, 10)], elty == Float64 ? randn(min(m, n) + 10) : complex.(randn(min(m, n) + 10), randn(min(m, n) + 10)))
        F = lu(A)
        L, U, p, q, Rs = F.:(:)
        @test (Diagonal(Rs) * A)[p,q] ≈ L * U
    end

    @testset "Issue #4523 - complex sparse \\" begin
        A, b = sparse((1.0 + im)I, 2, 2), fill(1., 2)
        @test A * (lu(A)\b) ≈ b

        @test det(sparse([1,3,3,1], [1,1,3,3], [1,1,1,1])) == 0
    end

    @testset "UMFPACK_ERROR_n_nonpositive" begin
        @test_throws ArgumentError lu(sparse(Int[], Int[], Float64[], 5, 0))
    end

    @testset "Issue #15099" for (Tin, Tout) in (
            (ComplexF16, ComplexF64),
            (ComplexF32, ComplexF64),
            (ComplexF64, ComplexF64),
            (Float16, Float64),
            (Float32, Float64),
            (Float64, Float64),
            (Int, Float64),
        )

        F = lu(sparse(fill(Tin(1), 1, 1)))
        L = sparse(fill(Tout(1), 1, 1))
        @test F.p == F.q == [1]
        @test F.Rs == [1.0]
        @test F.L == F.U == L
        @test F.:(:) == (L, L, [1], [1], [1.0])
    end

    @testset "BigFloat not supported" for T in (BigFloat, Complex{BigFloat})
        @test_throws ArgumentError lu(sparse(fill(T(1), 1, 1)))
    end

    @testset "size(::UmfpackLU)" begin
        m = n = 1
        F = lu(sparse(fill(1., m, n)))
        @test size(F) == (m, n)
        @test size(F, 1) == m
        @test size(F, 2) == n
        @test size(F, 3) == 1
        @test_throws ArgumentError size(F,-1)
    end

    @testset "Test aliasing" begin
        a = rand(5)
        @test_throws ArgumentError SuiteSparse.UMFPACK.solve!(a, lu(sparse(1.0I, 5, 5)), a, SuiteSparse.UMFPACK.UMFPACK_A)
        aa = complex(a)
        @test_throws ArgumentError SuiteSparse.UMFPACK.solve!(aa, lu(sparse((1.0im)I, 5, 5)), aa, SuiteSparse.UMFPACK.UMFPACK_A)
    end

    @testset "Issues #18246,18244 - lu sparse pivot" begin
        A = sparse(1.0I, 4, 4)
        A[1:2,1:2] = [-.01 -200; 200 .001]
        F = lu(A)
        @test F.p == [3 ; 4 ; 2 ; 1]
    end

    @testset "Test that A[c|t]_ldiv_B!{T<:Complex}(X::StridedMatrix{T}, lu::UmfpackLU{Float64}, B::StridedMatrix{T}) works as expected." begin
        N = 10
        p = 0.5
        A = N*I + sprand(N, N, p)
        X = zeros(ComplexF64, N, N)
        B = complex.(rand(N, N), rand(N, N))
        luA, lufA = lu(A), lu(Array(A))
        @test ldiv!(copy(X), luA, B) ≈ ldiv!(copy(X), lufA, B)
        @test ldiv!(copy(X), adjoint(luA), B) ≈ ldiv!(copy(X), adjoint(lufA), B)
        @test ldiv!(copy(X), transpose(luA), B) ≈ ldiv!(copy(X), transpose(lufA), B)
    end

    @testset "singular matrix" begin
        for A in sparse.((Float64[1 2; 0 0], ComplexF64[1 2; 0 0]))
            @test_throws SingularException lu(A)
            @test !issuccess(lu(A; check = false))
        end
    end

    @testset "deserialization" begin
        A  = 10*I + sprandn(10, 10, 0.4)
        F1 = lu(A)
        b  = IOBuffer()
        serialize(b, F1)
        seekstart(b)
        F2 = deserialize(b)
        for nm in (:colptr, :m, :n, :nzval, :rowval, :status)
            @test getfield(F1, nm) == getfield(F2, nm)
        end
    end

    @testset "Reuse symbolic LU factorization" begin
        A1 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
                    increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
                    [2.,1.,3.,4.,-1.,-3.,3.,9.,2.,1.,4.,2.], 5, 5)
        for Tv in (Float64, ComplexF64, Float16, Float32, ComplexF16, ComplexF32)
            for Ti in Base.uniontypes(SuiteSparse.UMFPACK.UMFITypes)
                A = convert(SparseMatrixCSC{Tv,Ti}, A0)
                B = convert(SparseMatrixCSC{Tv,Ti}, A1)
                b = Tv[8., 45., -3., 3., 19.]
                F = lu(A)
                lu!(F, B)
                @test F\b ≈ B\b ≈ Matrix(B)\b

                # singular matrix
                C = copy(B)
                C[4, 3] = Tv(0)
                F = lu(A)
                @test_throws SingularException lu!(F, C)

                # change of nonzero pattern
                D = copy(B)
                D[5, 1] = Tv(1.0)
                F = lu(A)
                @test_throws ArgumentError lu!(F, D)
            end
        end
    end

end

@testset "REPL printing of UmfpackLU" begin
    # regular matrix
    A = sparse([1, 2], [1, 2], Float64[1.0, 1.0])
    F = lu(A)
    facstring = sprint((t, s) -> show(t, "text/plain", s), F)
    lstring = sprint((t, s) -> show(t, "text/plain", s), F.L)
    ustring = sprint((t, s) -> show(t, "text/plain", s), F.U)
    @test facstring == "$(summary(F))\nL factor:\n$lstring\nU factor:\n$ustring"

    # singular matrix
    B = sparse(zeros(Float64, 2, 2))
    F = lu(B; check=false)
    facstring = sprint((t, s) -> show(t, "text/plain", s), F)
    @test facstring == "Failed factorization of type $(summary(F))"
end

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
# Examples presented in class - Lecture 1
# Example 2 - portfolio optimisation

using Ipopt, Plots
# Ipopt is for solving nonlinear problems - different technology
# Plots is for plotting (!)

# Read the daily prices data from a .csv file.
data = readcsv("prices.csv")
#First row has the names of the stocks
stock_names = data[1, 1:end-2]
#Last two columns has data and US$ rate
prices_data = data[2:end, 1:end-2]

#Having a peek at the data
cs = plot(prices_data, label = stock_names)

#Returns are calculated as (p(t+1) - p(t))/p(t-1)
returns_data = diff(prices_data) ./ prices_data[1:end-1,:]

#Number of days and stocks in data
T, n =  size(returns_data)
#Calculates expected return and covariance
μ, Σ =  mean(returns_data, 1), cov(returns_data)
# Input form the model: minimum average return required.
r_min = 0.2

port = Model(solver=IpoptSolver())

#allocation variables
@variable(port, 0 <= x[1:n] <= 1);

#notice the division by T to correct the average
@constraint(port, sum(μ[j]*x[j] for j=1:n) >= r_min/T);

@objective(port, Min , sum(x[i]*Σ[i,j]*x[j] for i=1:n,j=1:n));

solve(port);

alloc = getvalue(x)
# Should be the same as r_min
ret = μ*alloc*T
# Looking at risk as the st. deviation of the returns.
risk = sqrt(alloc'*Σ*alloc*T)

# Organising the data for plotting. In an array, they are
# considered different series (with different colours)
stock_names = convert(Array{String}, stock_names)
alloc = convert(Array{Float64}, alloc)
# plots a pie chart
pie(stock_names, alloc, legend = false)

ext: jl


text:
using HigherPrecision
using Base.Test

@testset "DoubleFloat64" begin
    @test DoubleFloat64(Float32(2.0)) == DoubleFloat64(2.0)
    @test DoubleFloat64(Float16(2.0)) == DoubleFloat64(2.0)
    @test DoubleFloat64(2) == DoubleFloat64(2.0)
    @test DoubleFloat64(big(π)) isa FastDouble
    @test DoubleFloat64(π) == DoubleFloat64(big(π))

    @test DoubleFloat64(2 // 3) isa FastDouble
    @test DoubleFloat64(2 // 3, ComputeAccurate) isa AccurateDouble
    @test DoubleFloat64(Float32(2.0), ComputeAccurate) isa AccurateDouble
    @test DoubleFloat64(Float16(2.0), ComputeAccurate) isa AccurateDouble
    @test DoubleFloat64(2, ComputeAccurate) isa AccurateDouble
    @test DoubleFloat64(big(π), ComputeAccurate) isa AccurateDouble
    @test DoubleFloat64(π, ComputeAccurate) isa AccurateDouble

    @test HigherPrecision.hi(DoubleFloat64(2.0)) == 2.0
    @test HigherPrecision.lo(DoubleFloat64(2.0)) == 0.0

    @test isbits(FastDouble) == true
    @test isbits(AccurateDouble) == true

    @test zero(DoubleFloat64(2.0)) == zero(FastDouble)

    @test convert(Integer, DoubleFloat64(2.0)) isa Int64
    @test 32 ≤ length(string(rand(FastDouble))) ≤ 40

    @test convert(DoubleFloat64, 2.0) isa FastDouble
    @test convert(DoubleFloat64, Float32(2.0)) isa FastDouble
    @test convert(DoubleFloat64, Float16(2.0)) isa FastDouble
    @test convert(DoubleFloat64, 2) isa FastDouble
    @test convert(DoubleFloat64, π) isa FastDouble

    @test promote_type(DoubleFloat64, Float64) == FastDouble
    @test promote_type(DoubleFloat64, Float32) == FastDouble
    @test promote_type(DoubleFloat64, Float16) == FastDouble

    @test double_sub(rand(), rand()) isa FastDouble
    @test double_add(rand(), rand()) isa FastDouble
    @test double_div(rand(), rand()) isa FastDouble
    @test double_mul(rand(), rand()) isa FastDouble

    x, y = rand(FastDouble), rand(AccurateDouble)

    @test x isa FastDouble
    @test y isa AccurateDouble
    @test x * y isa AccurateDouble
    @test y * x isa AccurateDouble
    @test x / y isa AccurateDouble
    @test y / x isa AccurateDouble
    @test x - y isa AccurateDouble
    @test y - x isa AccurateDouble
    @test x + y isa AccurateDouble
    @test y + x isa AccurateDouble
    @test x / 3.2 isa FastDouble
    @test y / 3.2 isa AccurateDouble

    @test x + 2 isa FastDouble
    @test 2 + x isa FastDouble

    @test im * rand(DoubleFloat64) isa Complex{FastDouble}
    @test x < 2 * x
    @test x ≤ x
    @test convert(Float64, x) ≤ x
    @test x ≤ 2 * convert(Float64, x)
    @test DoubleFloat64(2.0) == 2.0
    @test 2.0 == DoubleFloat64(2.0)

    @test isnan(DoubleFloat64(NaN))
    @test isinf(DoubleFloat64(Inf))

    @test floor(DoubleFloat64(3.2)) == 3.0
    @test ceil(DoubleFloat64(3.2)) == 4.0
    @test trunc(DoubleFloat64(3.2)) == 3.0
    @test isinteger(DoubleFloat64(3.2)) == false
    @test isinteger(DoubleFloat64(3.0)) == true

    @test DoubleFloat64(0.0)^0.14 == 1.0

    @test rand(Base.Random.GLOBAL_RNG, FastDouble, 4) isa Vector{FastDouble}
    @test rand(Base.Random.GLOBAL_RNG, Complex{FastDouble}, 4) isa Vector{Complex{FastDouble}}
    @test rand(Complex{DoubleFloat64}, 3) isa Vector{Complex{FastDouble}}
    @test rand(Base.Random.GLOBAL_RNG, Complex{DoubleFloat64}, 4) isa Vector{Complex{FastDouble}}
    @test length(rand(Base.Random.GLOBAL_RNG, FastDouble, 4)) == 4

    for k = 1:5000
        x = rand(DoubleFloat64) * 20 - 10
        y = rand(DoubleFloat64) * 20 - 10
        @test x * y ≈ big(x) * big(y) atol=1e-29
        @test x + y ≈ big(x) + big(y) atol=1e-30
        @test x - y ≈ big(x) - big(y) atol=1e-30
        @test x / y ≈ big(x) / big(y) atol=1e-26
        u = rand()
        @test x / u ≈ big(x) / big(u) atol=1e-26
        @test x^5 ≈ BigFloat(x)^5 atol=1e-25

        @test exp(x) ≈ exp(big(x)) atol=1e-26
        @test log(abs(x)) ≈ log(big(abs(x))) atol=1e-29
        @test sin(x) ≈ sin(big(x)) atol=1e-30
        @test cos(x) ≈ cos(big(x)) atol=1e-30
        @test atan(x) ≈ atan(big(x)) atol=1e-30
        @test atan2(y, x) ≈ atan2(big(y), big(x)) atol=1e-30
        @test tan(x) ≈ tan(big(x)) atol=1e-21

        sinh_x, cosh_x = sincosh(x)
        if abs(x) > 0.05
            @test sinh(x) ≈ sinh(big(x)) atol=1e-25
            @test tanh(x) ≈ tanh(big(x)) atol=1e-25
            @test sinh_x ≈ sinh(big(x)) atol=1e-25
            @test cosh_x ≈ cosh(big(x)) atol=1e-25
        else
            @test sinh(x) ≈ sinh(big(x)) atol=1e-16
            @test tanh(x) ≈ tanh(big(x)) atol=1e-16
            @test sinh_x ≈ sinh(big(x)) atol=1e-16
            @test cosh_x ≈ cosh(big(x)) atol=1e-16
        end
        @test cosh(x) ≈ cosh(big(x)) atol=1e-20

        sin_x, cos_x = HigherPrecision.sincos(x)
        @test sin_x ≈ sin(big(x)) atol=1e-30
        @test cos_x ≈ cos(big(x)) atol=1e-30

        z = DoubleFloat64(rand() * 2 - 1)
        @test asin(z) ≈ asin(big(z)) atol=1e-30
        @test acos(z) ≈ acos(big(z)) atol=1e-30

        n, p, d = Base.decompose(x)
        @test x ≈ BigFloat(n) * BigFloat(2)^p / BigFloat(d)
    end

    for k = 1:5000
        x = rand(DoubleFloat64{ComputeAccurate}) * 20 - 10
        y = rand(DoubleFloat64{ComputeAccurate}) * 20 - 10
        @test x * y ≈ big(x) * big(y) atol=1e-29
        @test x + y ≈ big(x) + big(y) atol=1e-31
        @test x - y ≈ big(x) - big(y) atol=1e-31
        @test x / y ≈ big(x) / big(y) atol=1e-28
        @test x^5 ≈ BigFloat(x)^5 atol=1e-26

        @test exp(x) ≈ exp(big(x)) atol=1e-26
        @test log(abs(x)) ≈ log(big(abs(x))) atol=1e-29
        @test sin(x) ≈ sin(big(x)) atol=1e-30
        @test cos(x) ≈ cos(big(x)) atol=1e-30
        @test atan(x) ≈ atan(big(x)) atol=1e-30
        @test atan2(y, x) ≈ atan2(big(y), big(x)) atol=1e-30
        @test tan(x) ≈ tan(big(x)) atol=1e-21

        sinh_x, cosh_x = sincosh(x)
        if abs(x) > 0.05
            @test sinh(x) ≈ sinh(big(x)) atol=1e-25
            @test tanh(x) ≈ tanh(big(x)) atol=1e-25
            @test sinh_x ≈ sinh(big(x)) atol=1e-25
            @test cosh_x ≈ cosh(big(x)) atol=1e-25
        else
            @test sinh(x) ≈ sinh(big(x)) atol=1e-16
            @test tanh(x) ≈ tanh(big(x)) atol=1e-16
            @test sinh_x ≈ sinh(big(x)) atol=1e-16
            @test cosh_x ≈ cosh(big(x)) atol=1e-16
        end
        @test cosh(x) ≈ cosh(big(x)) atol=1e-20

        sin_x, cos_x = HigherPrecision.sincos(x)
        @test sin_x ≈ sin(big(x)) atol=1e-30
        @test cos_x ≈ cos(big(x)) atol=1e-30

        z = DoubleFloat64(rand() * 2 - 1)
        @test asin(z) ≈ asin(big(z)) atol=1e-30
        @test acos(z) ≈ acos(big(z)) atol=1e-30
    end

    @test eps(DoubleFloat64(1.0)) == 4.93038065763132e-32 # 2^-104
    @test eps(DoubleFloat64) == 4.93038065763132e-32 # 2^-104
    @test eps(FastDouble) == 4.93038065763132e-32 # 2^-104
    @test eps(AccurateDouble) == 4.93038065763132e-32 # 2^-104
    @test realmin(DoubleFloat64) > realmin(Float64)
    @test realmax(DoubleFloat64) > realmax(Float64)
end


text:
struct Link
    href::URI
    rel::AbstractString
    name::AbstractString
    prompt::AbstractString
    render::AbstractString
end


text:
function _datareciperegistry(datasetregistry; name = "Dataset recipes")
    Registry(
        (;
            id = Field(String, name = "ID", formatfn = x -> sprint(show, x)),
            blocks = Field(
                Any,
                name = "Block types",
                description = "Types of blocks of the data container that this recipe loads.",
                filterfn = blocktypesmatch,
                formatfn = b -> FeatureRegistries.code_format(_formatblock(b))),
            description = Field(
                String,
                name = "Description",
                optional = true,
                description = "More information about the dataset recipe",
                formatfn=FeatureRegistries.md_format),
            downloaded = Field(
                Bool,
                name = "Is downloaded",
                description = """
                    Whether the dataset this recipe is based has been downloaded and is
                    available locally.
                """,
                computefn = (row, key) -> isavailable(datasetregistry[row.datasetid].loader)),
            datasetid = Field(
                String,
                name = "Dataset ID",
                description = "ID of the dataset this recipe is based on.",
                computefn = function (row, key)
                    val = row[key]
                    if !haskey(datasetregistry, val)
                        throw(ArgumentError("Could not find dataset with ID \"$val\"!"))
                    end
                    return val
                end,
            ),
            package = Field(
                Module,
                name = "Package"),
            recipe = Field(
                Datasets.DatasetRecipe,
                name = "Recipe",
                formatfn = FeatureRegistries.type_format
            )
        );
        name,
        loadfn = function loadrecipeentry(row)
            dataset = load(datasetregistry[row.datasetid])
            Datasets.loadrecipe(row.recipe, dataset)
        end,
        description = """
        A registry for dataset recipes. `load`ing an entry will download the
        corresponding dataset (see `datasets`) and return `data, blocks`, a
        data container and the `Block`s of the observations.

         ```julia
         data, blocks = load(datarecipes()[id])
         ```

        See `learningtasks` to find compatible learning tasks.
        """
    )
end

const DATARECIPES = _datareciperegistry(DATASETS)

"""
    datarecipes(; filters...)

Show a registry of available dataset recipes. A dataset recipe defines how
to load a dataset into a suitable format for use with a learning task.

Pass in filters as keyword arguments to look at a subset.

See also [finding functionality](/documents/docs/discovery.md), [`datasets`](#),
and [`learningtasks`](#). For more information on registries, see
[FeatureRegistries.jl](https://github.com/lorenzoh/FeatureRegistries.jl).

## Examples

Show all available dataset recipes:

{cell}
```julia
using FastAI
datarecipes()
```

Show all recipes for datasets that have "image" in their name:

{cell}
```julia
datarecipes(datasetid="image")
```

Show all data recipes usable for classification tasks, that is where the
target block is a [`Label`](#):

{cell}
```julia
datarecipes(blocks=(Any, Label))
```

Get an explanation of fields in the dataset recipe registry:

{cell}
```julia
info(datarecipes())
```
"""
datarecipes(; kwargs...) = isempty(kwargs) ? DATARECIPES : filter(DATARECIPES; kwargs...)


text:
# Debits & Credits

type UnbalancedException <: Exception end

function debit end
credit(args...) = -debit(args...)


text:
function count_histogram(nbr_outcomes, sequence)
    count = zeros(nbr_outcomes)
    for x in sequence
        count[x] += 1
    end
    count
end

using PyPlot

function simulate_die(ntrials)
    outcomes = collect(1:6)     # Simulate a fair die
    x = rand(outcomes, ntrials);
    bar(outcomes, count_histogram(6, x) / ntrials)
    xlabel("Die outcome")
    ylabel("Probability")
end

simulate_die(1000);

function simulate_sum_of_n_dice(ntrials, ndice)
    outcomes = collect(1:6)     # Simulate fair dice
    x = zeros(Int64, ntrials)
    for i = 1:ndice
        x .+= rand(outcomes, ntrials)
    end
    outcomesn = collect(1:6ndice)
    bar(outcomesn, count_histogram(6ndice, x) / ntrials)
    xlabel("Sum of n dice outcome")
    ylabel("Probability")
end

simulate_sum_of_n_dice(1000, 2);    # Two dice

simulate_sum_of_n_dice(10000, 50);

x = randn(10000)
nbins = 50

plt.hist(x, nbins)
xlabel("Random variable x")
ylabel("Count");


text:
#-----------------------------------------------------------------------
#   The functions presented here implement the relativistic location
#   algorithm of the authors
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   Levi-Civita
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   LEVI-CIVITA SYMBOL AND TENSOR
#-----------------------------------------------------------------------
"""
    ϵ( a , b , c , d )

The function `ϵ` computes the Levi-Civita symbol and implements the
corresponding tensor. The arguments `( a , b , c , d )` can either be
integers or vectors; in the case of integers, it returns the components
of the Levi-Civita symbol, and in the case of vectors, it returns the
contraction of the Levi-Civita symbol with the vectors.

This function also returns Hodge duals of vector tensor products. For
instance, the following

    ϵ( 0 , 0 , U , V )

yields the rank-2 quantity ``ϵ_{i,j,k,l} U^i V^i``. Make sure all `0`
valued arguments come first.

"""
function ϵ( a , b , c , d )              # Levi-Civita symbol and tensor
    if a == 0 && typeof(b) <: RealVec && typeof(c) <: RealVec &&
        typeof(d) <: RealVec
        v = zeros( typeof(d[1]) , 4 )
        for i=1:4, j=1:4, k=1:4, l=1:4
            v[i] += levicivita([i,j,k,l])*b[j]*c[k]*d[l]
        end
        return v
    elseif a == 0 && b == 0 && typeof(c) <: RealVec &&
        typeof(d) <: RealVec
        v = zeros( typeof(d[1]) , 4 , 4 )
        for i=1:4, j=1:4, k=1:4, l=1:4
            v[i,j] += levicivita([i,j,k,l])*c[k]*d[l]
        end
        return v
    elseif a == 0 && b == 0 && c == 0 && typeof(d) <: RealVec
        v = zeros( typeof(d[1]) , 4 , 4 , 4 )
        for i=1:4, j=1:4, k=1:4, l=1:4
            v[i,j,k] += levicivita([i,j,k,l])*d[l]
        end
        return v
    elseif typeof.((a,b,c,d)) == (Int,Int,Int,Int) && a != 0 || b != 0 ||
        c != 0 || d != 0
        return levicivita([a,b,c,d])
    elseif typeof(a) <: RealVec && typeof(b) <: RealVec &&
        typeof(c) <: RealVec && typeof(d) <: RealVec
        v = zero( typeof(a[1]) )
        for i=1:4, j=1:4, k=1:4, l=1:4
            v += levicivita([i,j,k,l])*a[i]*b[j]*c[k]*d[l]
        end
        return v
    else
        return 0
    end
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   HODGE VECTOR FUNCTION
#-----------------------------------------------------------------------
"""
    HodgeV( U::RealVec , V::RealVec , W::RealVec )

The function `HodgeV` computes `ϵ(0,U,V,W)`` and effectively raises the
index of the resulting (dual) vector.

"""
function HodgeV( U::RealVec , V::RealVec , W::RealVec )
    v = ϵ(0,U,V,W)
    v[1] = -v[1]        # Raise index
    return v
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   Transformation functions
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   LORENTZ TRANSFORMATION MATRIX CONSTRUCTOR
#-----------------------------------------------------------------------
"""
    LTM( NT::RealVec )

This function takes a timelike normal vector `NT` and constructs a
Lorentz transformation matrix so that the transformed vector `NT'` has
the form `NT'=[ N0 ; 0 ; 0 ; 0 ]`.

"""
function LTM( NT::RealVec )             # Lorentz transformation matrix
    tpfl=typeof(NT[1])
    normv = mnorm(NT)           # Normalization factor
    NV    = NT/normv            # Normalize Vector
    l     = one(tpfl)           # l has a value of one
    γ     = NV[1]               # γ is the Lorentz factor
    if γ<0
        NV = -NV
        γ  = -γ
    end
        δ    = (γ - l)             # A useful quantity
        norms = norm(NV[2:4])
    if γ != l || norms != zero(tpfl) || δ != zero(tpfl)
        nf = l/norms    # Normalization factor
        β = abs(norms/γ) # β is v/c

        vx = nf*NV[2]           # x-component of unit vector
        vy = nf*NV[3]           # y-component of unit vector
        vz = nf*NV[4]           # z-component of unit vector
        return [  γ        -γ*β*vx       -γ*β*vy       -γ*β*vz      ;
                 -γ*β*vx   l + δ*(vx^2)  δ*vx*vy       δ*vx*vz      ;
                 -γ*β*vy   δ*vy*vx       l + δ*(vy^2)  δ*vy*vz      ;
                 -γ*β*vz   δ*vz*vx       δ*vz*vy       l + δ*(vz^2) ]
    else
        return Matrix(l*I(4))
    end
end  # End LTM

#-----------------------------------------------------------------------
#   ROTATION MATRIX CONSTRUCTOR
#-----------------------------------------------------------------------
"""
    MRz( v::RealVec )

This function takes a vector `v` and constructs a rotation matrix so
that the transformed vector `v'` has the form `v'=[ 0 ; ... ; 0 ; vz ]`.
The three-dimensional and four-dimensional cases are considered.

"""
function MRz( v::RealVec )                   # Rotate to z-adapted frame
    tpfl=typeof(v[1])
    l = one(tpfl)
    o = zero(tpfl)
	if length(v)==3             # 3×3 rotation matrix
		rv=norm(v[1:3])
		Rv=norm(v[1:2])
        if Rv > 0 && rv > 0
		    Cosϑ=v[3]/rv
		    Sinϑ=Rv/rv
		    Cosφ=v[1]/Rv
		    Sinφ=v[2]/Rv
		    return [ Cosϑ*Cosφ  Cosϑ*Sinφ   -Sinϑ   ;
                     -Sinφ      Cosφ        o       ;
                     Sinϑ*Cosφ  Sinϑ*Sinφ   Cosϑ    ]
        else
            return Matrix(l*I(3))
        end
    elseif length(v)==4         # 4×4 rotation matrix
        rv=norm(v[2:4])
		Rv=norm(v[2:3])
        if Rv > 0 && rv > 0
            Cosϑ=v[4]/rv
            Sinϑ=Rv/rv
		    Cosφ=v[2]/Rv
		    Sinφ=v[3]/Rv
            return [ one(tpfl)  o           o           o       ;
                        o       Cosϑ*Cosφ   Cosϑ*Sinφ   -Sinϑ   ;
                        o       -Sinφ       Cosφ        o       ;
                        o       Sinϑ*Cosφ   Sinϑ*Sinφ   Cosϑ    ]
        else
            return Matrix(l*I(4))
        end
    else
        return Matrix(l*I(4))
	end
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   SPACETIME ROTATION OPERATOR
#-----------------------------------------------------------------------
"""
    Lrot( X::RealMtx , Λ::RealMtx )

This function applies the transformation transformation matrix `Λ` to
four emission points, which form the columns of a ``4×4`` matrix `X`.
The function returns a ``4×4`` matrix of transformed points `X'`.

"""
function Lrot( X::RealMtx , Λ::RealMtx )  # Spacetime rotation operator
    tpfl = typeof(X[1,1])

    Xprime = zeros(tpfl,4,4)

    for i=1:4
        Xprime[:,i] = Λ*X[:,i]
    end

    return Xprime
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   FLIP SPACELIKE VECTORS ACROSS A NULL DIRECTION
#-----------------------------------------------------------------------
"""
    NormflipS( Vsl::RealVec )

This function takes a spacelike vector `Vsl` and "flips" it across a
null direction to form a timelike vector. The null direction is chosen
to be the one tangent to the plane spanned by `Vsl` and the time axis
in the given coordinate system.

This function is used to construct a timelike vector from a spacelike
configuration vector (defined to be the vector normal to the plane
spanned by the emission points).

"""
function NormflipS( Vsl::RealVec )
    tpfl = typeof(Vsl[1])

    normVsl = mnorm(Vsl)
    normsq = ηdot(Vsl,Vsl)

    if Vsl[1] < 0                   # Ensures vector is future pointing
        Vsl = -Vsl
    end

    norms = norm(Vsl[2:4])          # Compute spatial norm

    if normsq > 0
        NVsl = Vsl/normVsl          # Normalization of spacetime vector
        nvsl = Vsl[2:4]/norms       # Normalized spatial part
        nlt  = NVsl[1]              # Temporal projection
        nls  = norm(NVsl[2:4])      # Spatial projection

        NVs  = zeros(tpfl,4)        # Construct and return flipped vec
        NVs[1]   = nls
        NVs[2:4] = nlt * nvsl
        return NVs
    else
        return Vsl/normVsl
    end
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   Adapted frame intersection point finders
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   INTERSECTION POINT FINDER (SPACELIKE CONFIGURATION HYPERPLANE)
#-----------------------------------------------------------------------
"""
    IPFinderS( Y::RealMtx )

Given a ``4×4`` matrix `Y` of emission points in an adapted frame
(defined such that the emission points have the same time coordinate),
this function computes the coordinates of the intersection point ``X_c``
for the future pointing light cones of the emission points.

The procedure here is to find the circumcenter of a sphere which
circumscribes the emission points. The spatial coordinates of the
circumcenter yields the spatial coordinates of the intersection point
``X_c`` and the time coordinate is given by the common time coordinate
of the emission points plus the radius of the circumsphere.

"""
function IPFinderS( Y::RealMtx )
    # Finds intersection of light cones in adapted frame
    tpfl=typeof(Y[1,1])

    x = zeros(tpfl,3,4)     # Create containers
    v = zeros(tpfl,3,3)

    x[1,:] = Y[2,:]         # Spatial points for corners of tetrahedron
    x[2,:] = Y[3,:]
    x[3,:] = Y[4,:]

    v[:,1] = x[:,2] - x[:,1]    # "Frame" vectors centered on x[:,1]
    v[:,2] = x[:,3] - x[:,1]
    v[:,3] = x[:,4] - x[:,1]

    B = [ x[:,2]⋅x[:,2] - x[:,1]⋅x[:,1] ; x[:,3]⋅x[:,3] - x[:,1]⋅x[:,1];
          x[:,4]⋅x[:,4] - x[:,1]⋅x[:,1] ] / 2    # Constructing B vector

    A = transpose([ v[:,1]  v[:,2]  v[:,3] ])    # Constructing A matrix

    xc = inv(A)*B           # Compute circumcenter
    rc = (norm(xc-x[:,1]) + norm(xc-x[:,2]) + # Distance to circumcenter
          norm(xc-x[:,3]) + norm(xc-x[:,4])) / 4
    tc = rc + Y[1,1]        # Compute time coordinate
    return [ tc ; xc[1] ; xc[2] ; xc[3] ]
end  # End IPfinder

#-----------------------------------------------------------------------
#   INTERSECTION POINT FINDER (TIMELIKE CONFIGURATION HYPERPLANE)
#-----------------------------------------------------------------------
"""
    IPFinderT( Y::RealMtx )

Given a ``4×4`` matrix `Y` of emission points in an adapted frame
(defined such that the emission points have the same z coordinate), this
function computes the coordinates of the two intersection points `` X_c
`` for the future pointing light cones of the emission points.

The procedure in this case is to find the vertex of the hyperboloid
which passes through the emission points. The coordinates of the vertex
on the `` (t,x,y) `` hyperplane yields the `` (t_c,x_c,y_c) ``
coordinates of the  intersection points, and the Minkowski distance `` R
`` between the vertex and the hyperboloid yields the `` z_c = z_e ± R ``
coordinate of the intersection points `` X_c ``, where `` z_e `` is the
common `` z `` coordinate for the emission points in the adapted frame.

"""
function IPFinderT( Y::RealMtx )
    # Finds intersection of light cones in adapted frame
    # Spacelike subconfiguration plane
    tpfl=typeof(Y[1,1])

    #   Defining variables
    z  = ( Y[4,1] + Y[4,2] + Y[4,3] + Y[4,4] )/4
    ( t1 , t2 , t3 ) = ( Y[1,1] , Y[1,2] , Y[1,3] )
    ( x1 , x2 , x3 ) = ( Y[2,1] , Y[2,2] , Y[2,3] )
    ( y1 , y2 , y3 ) = ( Y[3,1] , Y[3,2] , Y[3,3] )
    ( t4 , x4 , y4 ) = ( Y[1,4] , Y[2,4] , Y[3,4] )

    #   Denominator for time coordinate:
    dT = -(t3*x2*y1) + t4*x2*y1 + t2*x3*y1 - t4*x3*y1 - t2*x4*y1 +
        t3*x4*y1 + t3*x1*y2 - t4*x1*y2 - t1*x3*y2 + t4*x3*y2 +
        t1*x4*y2 - t3*x4*y2 - t2*x1*y3 + t4*x1*y3 + t1*x2*y3 -
        t4*x2*y3 - t1*x4*y3 + t2*x4*y3 + t2*x1*y4 - t3*x1*y4 -
        t1*x2*y4 + t3*x2*y4 + t1*x3*y4 - t2*x3*y4

    #   Denominator for spatial coordinates:
    dS = t3*x2*y1 - t4*x2*y1 - t2*x3*y1 + t4*x3*y1 + t2*x4*y1 -
         t3*x4*y1 - t3*x1*y2 + t4*x1*y2 + t1*x3*y2 - t4*x3*y2 -
         t1*x4*y2 + t3*x4*y2 + t2*x1*y3 - t4*x1*y3 - t1*x2*y3 +
         t4*x2*y3 + t1*x4*y3 - t2*x4*y3 - t2*x1*y4 + t3*x1*y4 +
         t1*x2*y4 - t3*x2*y4 - t1*x3*y4 + t2*x3*y4

    #   t coordinate:
    tc  = ( -(t3^2*x2*y1) + t4^2*x2*y1 + t2^2*x3*y1 - t4^2*x3*y1 -
            x2^2*x3*y1 + x2*x3^2*y1 - t2^2*x4*y1 + t3^2*x4*y1 +
            x2^2*x4*y1 - x3^2*x4*y1 - x2*x4^2*y1 + x3*x4^2*y1 +
            t3^2*x1*y2 - t4^2*x1*y2 - t1^2*x3*y2 + t4^2*x3*y2 +
            x1^2*x3*y2 - x1*x3^2*y2 + t1^2*x4*y2 - t3^2*x4*y2 -
            x1^2*x4*y2 + x3^2*x4*y2 + x1*x4^2*y2 - x3*x4^2*y2 +
            x3*y1^2*y2 - x4*y1^2*y2 - x3*y1*y2^2 + x4*y1*y2^2 -
            t2^2*x1*y3 + t4^2*x1*y3 + t1^2*x2*y3 - t4^2*x2*y3 -
            x1^2*x2*y3 + x1*x2^2*y3 - t1^2*x4*y3 + t2^2*x4*y3 +
            x1^2*x4*y3 - x2^2*x4*y3 - x1*x4^2*y3 + x2*x4^2*y3 -
            x2*y1^2*y3 + x4*y1^2*y3 + x1*y2^2*y3 - x4*y2^2*y3 +
            x2*y1*y3^2 - x4*y1*y3^2 - x1*y2*y3^2 + x4*y2*y3^2 +
            t2^2*x1*y4 - t3^2*x1*y4 - t1^2*x2*y4 + t3^2*x2*y4 +
            x1^2*x2*y4 - x1*x2^2*y4 + t1^2*x3*y4 - t2^2*x3*y4 -
            x1^2*x3*y4 + x2^2*x3*y4 + x1*x3^2*y4 - x2*x3^2*y4 +
            x2*y1^2*y4 - x3*y1^2*y4 - x1*y2^2*y4 + x3*y2^2*y4 +
            x1*y3^2*y4 - x2*y3^2*y4 - x2*y1*y4^2 + x3*y1*y4^2 +
            x1*y2*y4^2 - x3*y2*y4^2 - x1*y3*y4^2 + x2*y3*y4^2 ) /
          (2*dT)

    #   x coordinate:
    xc = ( -(t2^2*t3*y1) + t2*t3^2*y1 + t2^2*t4*y1 - t3^2*t4*y1 -
            t2*t4^2*y1 + t3*t4^2*y1 + t3*x2^2*y1 - t4*x2^2*y1 -
            t2*x3^2*y1 + t4*x3^2*y1 + t2*x4^2*y1 - t3*x4^2*y1 +
            t1^2*t3*y2 - t1*t3^2*y2 - t1^2*t4*y2 + t3^2*t4*y2 +
            t1*t4^2*y2 - t3*t4^2*y2 - t3*x1^2*y2 + t4*x1^2*y2 +
            t1*x3^2*y2 - t4*x3^2*y2 - t1*x4^2*y2 + t3*x4^2*y2 -
            t3*y1^2*y2 + t4*y1^2*y2 + t3*y1*y2^2 - t4*y1*y2^2 -
            t1^2*t2*y3 + t1*t2^2*y3 + t1^2*t4*y3 - t2^2*t4*y3 -
            t1*t4^2*y3 + t2*t4^2*y3 + t2*x1^2*y3 - t4*x1^2*y3 -
            t1*x2^2*y3 + t4*x2^2*y3 + t1*x4^2*y3 - t2*x4^2*y3 +
            t2*y1^2*y3 - t4*y1^2*y3 - t1*y2^2*y3 + t4*y2^2*y3 -
            t2*y1*y3^2 + t4*y1*y3^2 + t1*y2*y3^2 - t4*y2*y3^2 +
            t1^2*t2*y4 - t1*t2^2*y4 - t1^2*t3*y4 + t2^2*t3*y4 +
            t1*t3^2*y4 - t2*t3^2*y4 - t2*x1^2*y4 + t3*x1^2*y4 +
            t1*x2^2*y4 - t3*x2^2*y4 - t1*x3^2*y4 + t2*x3^2*y4 -
            t2*y1^2*y4 + t3*y1^2*y4 + t1*y2^2*y4 - t3*y2^2*y4 -
            t1*y3^2*y4 + t2*y3^2*y4 + t2*y1*y4^2 - t3*y1*y4^2 -
            t1*y2*y4^2 + t3*y2*y4^2 + t1*y3*y4^2 - t2*y3*y4^2) /
          (2*dS)

    #   y coordinate:
    yc = (  t2^2*t3*x1 - t2*t3^2*x1 - t2^2*t4*x1 + t3^2*t4*x1 +
            t2*t4^2*x1 - t3*t4^2*x1 - t1^2*t3*x2 + t1*t3^2*x2 +
            t1^2*t4*x2 - t3^2*t4*x2 - t1*t4^2*x2 + t3*t4^2*x2 +
            t3*x1^2*x2 - t4*x1^2*x2 - t3*x1*x2^2 + t4*x1*x2^2 +
            t1^2*t2*x3 - t1*t2^2*x3 - t1^2*t4*x3 + t2^2*t4*x3 +
            t1*t4^2*x3 - t2*t4^2*x3 - t2*x1^2*x3 + t4*x1^2*x3 +
            t1*x2^2*x3 - t4*x2^2*x3 + t2*x1*x3^2 - t4*x1*x3^2 -
            t1*x2*x3^2 + t4*x2*x3^2 - t1^2*t2*x4 + t1*t2^2*x4 +
            t1^2*t3*x4 - t2^2*t3*x4 - t1*t3^2*x4 + t2*t3^2*x4 +
            t2*x1^2*x4 - t3*x1^2*x4 - t1*x2^2*x4 + t3*x2^2*x4 +
            t1*x3^2*x4 - t2*x3^2*x4 - t2*x1*x4^2 + t3*x1*x4^2 +
            t1*x2*x4^2 - t3*x2*x4^2 - t1*x3*x4^2 + t2*x3*x4^2 +
            t3*x2*y1^2 - t4*x2*y1^2 - t2*x3*y1^2 + t4*x3*y1^2 +
            t2*x4*y1^2 - t3*x4*y1^2 - t3*x1*y2^2 + t4*x1*y2^2 +
            t1*x3*y2^2 - t4*x3*y2^2 - t1*x4*y2^2 + t3*x4*y2^2 +
            t2*x1*y3^2 - t4*x1*y3^2 - t1*x2*y3^2 + t4*x2*y3^2 +
            t1*x4*y3^2 - t2*x4*y3^2 - t2*x1*y4^2 + t3*x1*y4^2 +
            t1*x2*y4^2 - t3*x2*y4^2 - t1*x3*y4^2 + t2*x3*y4^2 ) /
          (2*dS)

    #   Hyperboloid distance calculation
    R  = (  mnorm( [t1-tc;x1-xc;y1-yc] ) +
            mnorm( [t2-tc;x2-xc;y2-yc] ) +
            mnorm( [t3-tc;x3-xc;y3-yc] ) +
            mnorm( [t4-tc;x4-xc;y4-yc] ) ) / 4

    #   z coordinate
    Δz = R

    return  (   [ tc ; xc ; yc ; z + Δz ]
              , [ tc ; xc ; yc ; z - Δz ] )
end  # End IPfinder

#-----------------------------------------------------------------------
#   Locator functions
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   FOUR POINT LOCATOR FUNCTION (FHC22)
#-----------------------------------------------------------------------
"""
    locator4FHC22( X::RealMtx )

This function implements the four point relativistic location algorithm
of the authors. It outputs a pair of location points.

"""
function locator4FHC22( X::RealMtx )
    #   Computes location from a single set of four emission points
    tpfl = typeof(X[1,1])
    XF = tpfl.(X)

    #   Containers
        nv  = zeros(tpfl,4)
        Xc  = zeros(tpfl,4)

        XA = zeros(tpfl,4,4)
        XB = zeros(tpfl,4,4)

        Mz  = zeros(tpfl,4,4)
        My  = zeros(tpfl,4,4)
        Λ   = zeros(tpfl,4,4)

        nv = HodgeV(X[:,1] - X[:,4],X[:,2] - X[:,4],X[:,3] - X[:,4])

    if ηdot(nv,nv) < 0         # Timelike normal vector n
        #---------------------------------------------------------------
        #   Computation for spacelike configuration plane
        #---------------------------------------------------------------
        Λ  = LTM(nv)
        Xc = inv(Λ)*IPFinderS( Lrot(XF,Λ) )
        return (Xc,Xc)
    elseif ηdot(nv,nv) > 0     # Spacelike normal vector n
        #---------------------------------------------------------------
        #   Computation for timelike configuration plane
        #---------------------------------------------------------------
        Λn  = LTM(NormflipS(nv))
        Mz  = MRz(nv)
        Λ   = Mz*Λn
        XT  = IPFinderT( Lrot(XF,Λ) )
        return (inv(Λ)*XT[1],inv(Λ)*XT[2])
    elseif ηdot(nv,nv) == 0     # Null normal vector n
		print("nv norm zero.")
        return (Xc,Xc)
    end
end     #---------------------------------------------------------------


text:
using Documenter, MLPreprocessing

makedocs(
    modules = [MLPreprocessing],
    clean = false,
    format = :html,
    assets = [
        joinpath("assets", "favicon.ico"),
    ],
    sitename = "MLPreprocessing.jl",
    authors = "Christof Stocker, Andre Bieler",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
        hide("Indices" => "indices.md"),
        "LICENSE.md",
    ],
    html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo = "github.com/JuliaML/MLPreprocessing.jl.git",
    target = "build",
    julia = "0.6",
    deps = nothing,
    make = nothing,
)


text:
"""
    FieldMaximum(mapping, field)

An object for calculating the maximum of a `mapping` function applied
element-wise to `field`.

Examples
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)));

julia> max_abs_u = FieldMaximum(abs, model.velocities.u);

julia> max_w² = FieldMaximum(x->x^2, model.velocities.w);
```
"""
struct FieldMaximum{F, M}
    mapping :: M
      field :: F
end

(m::FieldMaximum)(args...) = maximum(m.mapping, m.field.data.parent)

(m::FieldMaximum{<:NamedTuple})(args...) =
    NamedTuple{propertynames(m.field)}(maximum(m.mapping, f.data.parent) for f in m.field)


text:
export Receive, Pay, Buy, Sell, ZCB, Forward, Option, EuropeanCall, EuropeanPut, AmericanCall, AmericanPut, Scale

abstract type Contract end

"""
    Zero() <: Contract

A "null" contract.
"""
struct Zero <: Contract
end

"""
    Amount(o::Observable) <: Contract

Receive amount of `o`.
"""
struct Amount{O} <: Contract
    o::O
end

"""
    Scale(s::Observable, c::Contract)

Scale the contract `c` by `s`.
"""
struct Scale{O<:Observable, C<:Contract} <: Contract
    s::O
    c::C
end
Scale(x::Number, c::Contract) = Scale(ConstObs(x),c)

"""
    Both(c1::Contract, c2::Contract)

Acquire both contracts `c1` and `c2` (`and` in PJE2000).
"""
struct Both{C1<:Contract,C2<:Contract} <: Contract
    c1::C1
    c2::C2
end

"""
    Either(c1::Contract, c2::Contract)

Acquire either contract `c1` or `c2` (`or` in PJE2000).
"""
struct Either{C1<:Contract,C2<:Contract} <: Contract
    c1::C1
    c2::C2
end

"""
    Give(c::Contract)

Take the opposite side of contract `c`.
"""
struct Give{C<:Contract} <: Contract
    c::C
end

"""
    Cond(p::Observable{Bool}, c1::Contract, c2::Contract)

If `p` is true at the point of aquisition, acquire `c1` else acquire `c2`.
"""
struct Cond{P<:Observable{Bool}, T1<:Contract, T2<:Contract} <: Contract
    p::P
    c1::T1
    c2::T2
end

"""
    When(p::Observable{Bool}, c::Contract)

Acquire the contract `c` at the point when `p` becomes `true`.
"""
struct When{P<:Observable{Bool},C<:Contract} <: Contract
    p::P
    c::C
end

"""
    Anytime(p::Observable{Bool}, c::Contract)

May acquire `c` at any point when `p` is `true`.
"""
struct Anytime{P<:Observable{Bool}, C<:Contract} <: Contract
    p::P
    c::C
end

"""
    Until(p::Observable{Bool}, c::Contract)

Acts like contract `c` until `p` is `true`, at which point it is abandoned (and hence worthless).
"""
struct Until{P<:Observable{Bool}, C<:Contract} <: Contract
    p::P
    c::C
end

## Derived Contracts

const Receive{T} = Amount{ConstObs{T}}
Receive(x::Union{Real,CurrencyQuantity}) = Amount(ConstObs(x))

const Pay{T} = Give{Receive{T}}
Pay(x::Union{Real,CurrencyQuantity}) = Give(Receive(x))

const Buy{C,T} = Both{C, Pay{T}}
Buy(c::Contract, x::Union{Real,CurrencyQuantity}) = Both(c, Pay(x))

const Sell{C,T} = Both{Give{C}, Receive{T}}
Sell(c::Contract, x::Union{Real,CurrencyQuantity}) = Both(Give(c), Receive(x))

const WhenAt{C} = When{At, C}
WhenAt(date::Date, c::Contract) = When(At(date), c)

const ZCB{T} = WhenAt{Receive{T}}
ZCB(date::Date, x::Union{Real,CurrencyQuantity}) = WhenAt(date, Receive(x))

const Forward{C,T} = WhenAt{Buy{C,T}}
Forward(date::Date, c::Contract, strike::Union{Real,CurrencyQuantity}) = WhenAt(date, Buy(c, strike))

const Option{C} = Either{C, Zero}
Option(c::Contract) = Either(c, Zero())

const European{C} = WhenAt{Option{C}}
European(date::Date, c::Contract) = WhenAt(date, Option(c))

"""
    EuropeanCall(date, c, strike)

A European call contract, with maturity `date`, on underlying contract `c` at price `strike`.
"""
const EuropeanCall{C,T} = European{Buy{C,T}}
EuropeanCall(date::Date, c::Contract, strike::Union{Real,CurrencyQuantity}) = European(date, Buy(c, strike))

"""
    EuropeanPut(date, c, strike)

A European put contract, with maturity `date`, on underlying contract `c` at price `strike`.
"""
const EuropeanPut{C,T} = European{Sell{C,T}}
EuropeanPut(date::Date, c::Contract, strike::Union{Real,CurrencyQuantity}) = European(date, Sell(c, strike))

const AnytimeBefore{C} = Anytime{Before, C}
AnytimeBefore(date::Date, c::Contract) = Anytime(Before(date), c)

const American{C} = AnytimeBefore{Option{C}}
American(date::Date, c::Contract) = AnytimeBefore(date, Option(c))

"""
    AmericanCall(date, c, strike)

An American call contract, with maturity `date`, on underlying contract `c` at price `strike`.
"""
const AmericanCall{C,T} = American{Buy{C,T}}
AmericanCall(date::Date, c::Contract, strike::Union{Real,CurrencyQuantity}) = American(date, Buy(c, strike))

"""
    AmericanPut(date, c, strike)

An American put contract, with maturity `date`, on underlying contract `c` at price `strike`.
"""
const AmericanPut{C,T} = American{Sell{C,T}}
AmericanPut(date::Date, c::Contract, strike::Union{Real,CurrencyQuantity}) = American(date, Sell(c, strike))

## helper functions

"""
    maturitydate(c::Contract)

The maturity date of a contract `c` is the date at which contract is completed or expires.
"""
maturitydate(c::WhenAt{C}) where {C} = c.p.a[2].val
maturitydate(c::AnytimeBefore{C}) where {C} = c.p.a[2].val

"""
    strikeprice(c::Contract)

The price at which an option or forward can be exercised.
"""
strikeprice(c::Option{C}) where {C} = price(c.c1)
strikeprice(c::American{C}) where {C} = strikeprice(c.c)
strikeprice(c::European{C}) where {C} = strikeprice(c.c)
strikeprice(c::WhenAt{C}) where {C} = price(c.c)

"""
    price(c::Contract)

The price of a `Buy` or `Sell` contract.
"""
price(c::Buy{C,T}) where {C,T} = c.c2.c.o.val
price(c::Sell{C,T}) where {C,T} = c.c2.o.val


text:
# ---
# title: 434. Number of Segments in a String
# id: problem434
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/number-of-segments-in-a-string/description/>
# hidden: true
# ---
#
# You are given a string `s`, return _the number of segments in the string_.
#
# A **segment** is defined to be a contiguous sequence of **non-space
# characters**.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "Hello, my name is John"
#     Output: 5
#     Explanation: The five segments are ["Hello,", "my", "name", "is", "John"]
#
#
# **Example 2:**
#
#
#
#     Input: s = "Hello"
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: s = "love live! mu'sic forever"
#     Output: 4
#
#
# **Example 4:**
#
#
#
#     Input: s = ""
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 300`
#   * `s` consists of lower-case and upper-case English letters, digits or one of the following characters `"!@#/$/%^&*()_+-=',.:"`.
#   * The only space character in `s` is `' '`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using BioSequences.FASTQ
using CodecZlib

# Zlib could possibly be made more efficient with transcodingstreams
# https://github.com/bicycle1885/TranscodingStreams.jl
# NOTE: if you assume that each FASTQ files are sorted in same order then
# you do not have to hash
# NOTE: Implementation is for single set of files
# Assumes file names will be:
# $input_1.fq.gz
# $input_2.fq.gz
# $input_1.tagged_filter.fastq.gz
# $input_2.tagged_filter.fastq.gz
# the search function is used to build the union and map_search is the intersection
function match_screens(input, out; merge_type = "intersection", raw_fq=".fq.gz", filt_fq = ".tagged_filter.fastq.gz", out_fq=".fs.fq.gz")
    r1_map = make_map(string(input,"_1", filt_fq))
    println("Map for read1 completed")
    r2_map = make_map(string(input,"_2", filt_fq))
    println("Map for read2 completed")
    writer_r1 = FASTQ.Writer(GzipCompressorStream(open(string(out,"_1",out_fq),"w")))
    writer_r2 = FASTQ.Writer(GzipCompressorStream(open(string(out,"_2", out_fq),"w")))
    map_search!(r1_map,r2_map,writer_r1,writer_r2)
    if merge_type == "union"
        search!(string(input,"_2", raw_fq),r1_map,writer_r2)
        search!(string(input,"_1". raw_fq),r2_map,writer_r1)
    end
    close(writer_r1)
    close(writer_r2)
end

# search thru r1_map for records in r2_map. If matched, write to file and
# delete record from both r1_map and r2_map
function map_search!(r1_map, r2_map, w1, w2)
    #i = 0 # TESTING THAT ALL FUNCTIONS ARE WORKING
    for record in values(r2_map)
        header = FASTQ.identifier(record)
        #p = parse(header)
        paired_rec = pop!(r1_map,header,0)
        if paired_rec != 0
            write(w1, record)
            write(w2, paired_rec)
            delete!(r2_map,header)
        end
    end
end
# search thru reader for records with identifiers in rmap
# rmap will have records deleted as we go based on if they are found or not
# For primer on mutable functions in Julia: https://docs.julialang.org/en/latest/manual/faq/?highlight=Mutable#functions
function search!(file, rmap, w1, w2)
    #i = 0
    reader = FASTQ.Reader(GzipDecompressorStream(open(file,"r")))
    for record in reader
        header = FASTQ.identifier(record)
        #p = parse(header)
        paired_rec = pop!(rmap,header,0)
        if paired_rec != 0
            write(w1, record)
            write(w2, paired_rec)
        end
    end
    close(reader)
end
# make read map
# reader is FASTQ Reader object
function make_map(file)
    reader = FASTQ.Reader(GzipDecompressorStream(open(file,"r")))
    rmap = Dict()
    for record in reader
        header = FASTQ.identifier(record)
        #p = parse(header)
        rmap[header] = record
    end
    return(rmap)
end
# parse header in some way that is mappable to the other read
# appears to not be needed for now
#function parse(header)
#end
#@time #match_screens("/Users/aragaven/scratch/test_julia/test_fq_screen/D1137A_1m", "/Users/aragaven/scratch/test_julia/test_fq_screen/D1137A", filt_fq="_tr.tagged_filter.fastq.gz")
@time match_screens(ARGS[1], ARGS[2], filt_fq="_tr.tagged_filter.fastq.gz")

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

# See https://github.com/JuliaLang/Pkg.jl/issues/2942
# Once this Pkg issue is resolved, this must be removed
uuid = Base.UUID("a83860b7-747b-57cf-bf1f-3e79990d037f")
delete!(Pkg.Types.get_last_stdlibs(v"1.6.3"), uuid)

julia_versions = [v"1.6.3", v"1.7.0"]

name = "OpenSpiel"
version = v"1.1.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/deepmind/open_spiel/archive/v1.1.1.tar.gz", "73d46e2d9ce7a86a420aad141b789e98237463b5c98608401faf86b15ec5eaf9"),
    ArchiveSource("https://github.com/findmyway/dds/archive/v0.1.1.tar.gz", "fd13ee77feb5b5c3dfcc3333a0523266beb2a3d27715703cf508313af25306e5"),
    ArchiveSource("https://github.com/abseil/abseil-cpp/archive/20211102.0.tar.gz", "dcf71b9cba8dc0ca9940c4b316a0c796be8fab42b070bb6b7cab62b48f0e66c4"),
    ArchiveSource("https://github.com/findmyway/hanabi-learning-environment/archive/v0.1.0.tar.gz", "6126936fd13a95f8cadeacaa69dfb38a960eaf3bd588aacc8893a6e07e4791a3"),
    ArchiveSource("https://github.com/findmyway/project_acpc_server/archive/v0.1.0.tar.gz", "e29f969dd62ba354b7019cae3f7f1dbfbd9a744687ea4a8f7494c2bb1ee87382"),
]

# Bash recipe for building across all platforms
script = raw"""

mv open_spiel-* open_spiel
mv abseil-cpp-* open_spiel/open_spiel/abseil-cpp
mv dds-* open_spiel/open_spiel/games/bridge/double_dummy_solver
mv hanabi-learning-environment-* open_spiel/open_spiel/games/hanabi/hanabi-learning-environment
mv project_acpc_server-* open_spiel/open_spiel/games/universal_poker/acpc

mkdir open_spiel/build
cd open_spiel/build
export OPEN_SPIEL_BUILD_WITH_JULIA=ON OPEN_SPIEL_BUILD_WITH_PYTHON=OFF OPEN_SPIEL_BUILD_WITH_HANABI=ON OPEN_SPIEL_BUILD_WITH_ACPC=OFF
cmake \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DJulia_PREFIX=${prefix} \
    ../open_spiel/
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/open_spiel/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc", julia_version=v)
    for v in julia_versions
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libspieljl", :libspieljl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"; platforms=platforms),
    BuildDependency("libjulia_jll"; platforms=platforms)
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"9",
    julia_compat="1.6"
)

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
@recipe function f(p::PowseekerProblem)
    n = 100
    x = linspace(-4000, 4000, n)
    y = linspace(-4000, 4000, n)
    ratio --> :equal
    seriestype --> :contour
    x, y, mdp(p).topology
end

@recipe function f(s::SkierState)
    len = 0.2
    u = len*cos(s.psi)
    v = len*sin(s.psi)
    seriestype --> :quiver
    quiver --> ([u], [v])
    [s.pos[1]], [s.pos[2]]
end

@recipe function f(v::AbstractVector{SkierState})
    len = 200
    psis = collect(s.psi for s in v)
    xs = collect(s.pos[1] for s in v)
    ys = collect(s.pos[2] for s in v)
    us = len*cos.(psis)
    vs = len*sin.(psis)
    @series begin
        label --> "True Path"
        xs, ys
    end
    @series begin
        label --> "direction"
        seriestype --> :quiver
        quiver --> (us, vs)
        xs, ys
    end
end

@recipe function f(hist::AbstractMDPHistory{SkierState, Float64})
    @series begin
        state_hist(hist)
    end
    len = 200
    sh = state_hist(hist)[1:end-1]
    psis = collect(s.psi for s in sh)
    xs = collect(s.pos[1] for s in sh)
    ys = collect(s.pos[2] for s in sh)
    aus = len*cos(psis + collect(action_hist(hist)))
    avs = len*sin(psis + collect(action_hist(hist)))
    @series begin
        label --> "action"
        seriestype --> :quiver
        quiver --> (aus, avs)
        xs, ys
    end
end

@recipe function f(b::ParticleCollection{SkierState})
    v = particles(b)
    len = 50
    psis = collect(s.psi for s in v)
    xs = collect(s.pos[1] for s in v)
    ys = collect(s.pos[2] for s in v)
    us = len*cos.(psis)
    vs = len*sin.(psis)
    @series begin
        color --> "black"
        seriestype --> :quiver
        quiver := (us, vs)
        xs, ys
    end
end

@recipe function f(h::AbstractPOMDPHistory{SkierState, GPSOrAngle})
    @series begin
        state_hist(h)[1:end-1]
    end
    @series begin
        belief_hist(h)[end]
    end
    oh = observation_hist(h)[1:end-1]
    gps_x = [get(o.pos)[1] for o in oh if !isnull(o.pos)]
    gps_y = [get(o.pos)[2] for o in oh if !isnull(o.pos)]
    if !isempty(gps_x)
        @series begin
            seriestype := :scatter
            marker --> :circle
            label := "GPS Measurements"
            gps_x, gps_y
        end
    end
end


text:
@testset "show" begin
    Mag = Arblib.Mag
    @testset "string" begin
        @test Arblib._string(Mag()) isa String
        @test !isempty(Arblib._string(Mag(2.3)))
        @test Arblib.string_nice(Arb()) isa String
        @test Arblib.string_nice(Acb()) isa String
    end

    @testset "dump" begin
        x = Mag(1.1)
        for x in (Mag(1.1), Arf(1.1), Arb(1.1))
            y = zero(x)
            str = Arblib.dump_string(x)
            Arblib.load_string!(y, str)
            @test isequal(x, y)
        end
    end

    @testset "show" begin
        prec = 32

        P = ArbPoly(Arb[1, 2, 0, π], prec = prec)
        @test "$P" == "1.00000000 + 2.00000000⋅x + [3.14159265 +/- 3.59e-9]⋅x^3"
        P = AcbPoly([Acb[1, 2, 0, π]; Acb(1, 1)], prec = prec)
        @test "$P" ==
              "1.00000000 + 2.00000000⋅x + [3.14159265 +/- 3.59e-9]⋅x^3 + (1.00000000 + 1.00000000im)⋅x^4"
        P = ArbSeries(Arb[1, 2, 0, π], degree = 4, prec = prec)
        @test "$P" == "1.00000000 + 2.00000000⋅x + [3.14159265 +/- 3.59e-9]⋅x^3 + 𝒪(x^5)"
        P = AcbSeries([Acb[1, 2, 0, π]; Acb(1, 1)], degree = 5, prec = prec)
        @test "$P" ==
              "1.00000000 + 2.00000000⋅x + [3.14159265 +/- 3.59e-9]⋅x^3 + (1.00000000 + 1.00000000im)⋅x^4 + 𝒪(x^6)"

        @test "$(ArbPoly())" == "$(AcbPoly())" == "0"
        @test "$(ArbSeries())" == "$(AcbSeries())" == "𝒪(x)"
    end
end


text:
struct CKAN <: DataRepo
end

description(repo::CKAN, mainpage) = mainpage["notes"]

author(repo::CKAN, mainpage) = [mainpage["author"]]

maintainer(repo::CKAN, mainpage) = mainpage["maintainer"]

license(repo::CKAN, mainpage) = mainpage["license_title"]

function create_date(repo::CKAN, mainpage)
    rawdate = Dates.Date(mainpage["metadata_created"][1:10], "yyyy-mm-dd")
    Dates.format(rawdate, "U d, yyyy")
end

function get_urls(repo::CKAN, page)
    urls = []
    for i = 1:page["num_resources"]
        push!(urls, page["resources"][i]["url"])
    end
    urls
end

data_fullname(::CKAN, mainpage) = mainpage["title"]

function website(repo::CKAN, mainpage_url, mainpage)
    replace(mainpage_url, "/api/3/action/package_show?id=" =>  "/dataset/")
end

function mainpage_url(repo::CKAN, dataname)
    if startswith(dataname, "http")
        url = replace(dataname, "/dataset/" =>  "/api/3/action/package_show?id=")
        return JSON.parse(text_only(getpage(url).root))["result"], dataname
    else
        throw(GeneratorError(repo))
    end
end


text:
# ps3000_open_unit.jl

export ps3000_open_unit

function ps3000_open_unit()
  ps_status = ccall((:ps3000_open_unit, ps3000driver),Int16,())
  if ps_status == -1
    error("ps3000 Error: Failed To Open")
  end
  return ps_status
end


text:
############################# abstract QR type #################################
abstract type AbstractQR{T} <: Factorization{T} end
# AbstractQR assumes existence of n, m, Q, R fields
Base.size(F::AbstractQR) = (F.n, F.m)
Base.size(F::AbstractQR, i::Int) = i > 2 ? 1 : size(F)[i]
Base.eltype(F::AbstractQR{T}) where T = T

Base.:\(F::AbstractQR, x::AbstractVecOrMat) = ldiv!(F, copy(x))

# some helpers
"""
```
    number_of_rotations(n::Int, m::Int)
```
Computes the number of Givens rotations that are necessary to compute the QR
factorization of a general matrix of size n by m.
"""
number_of_rotations(n::Int, m::Int) = n*m - (m*(m+1)) ÷ 2
number_of_rotations(A::AbstractMatOrFac) = number_of_rotations(size(A)...)
# min_nm = min(m, m) # IDEA: for overdetermined systems

# to append k columns
number_of_rotations_to_append_column(n::Int, m::Int, k::Int = 1) = k*(n-m) - (k*(k+1)) ÷ 2
# to remove a single column at index k in R
number_of_rotations_to_remove_column(m::Int, k::Int) = m-k

"""
```
    allocate_rotations(T::DataType, n::Int, m::Int)
```
Allocates a vector of Givens rotations of a length that is necessary to compute
the QR factorization of a general matrix of size n by m.
"""
function allocate_rotations(T::DataType, n::Int, m::Int)
    Vector{Givens{T}}(undef, number_of_rotations(n, m))
end
function allocate_rotations(A::AbstractMatOrFac)
    allocate_rotations(eltype(A), size(A)...)
end