text:
module PLINK

using StatsFuns
### idea from Julia Hoffimann Mendes to conditionally load OpenCL module
# if no OpenCL library is available, then issue a warning
# set "cl" variable to Void,
# then conditionally load GPU code based on value of "cl"
try
    using OpenCL
catch e
    warn("PLINK.jl cannot find an OpenCL library and will not load GPU functions correctly.")
    global cl = nothing
end
using DataFrames

import Base.size
import Base.==
import Base.isequal
import Base.mean
import Base.copy
import Base.getindex
import Base.length
import Base.ndims
import Base.display
import Base.convert
import Base.A_mul_B!
import Base.At_mul_B!
import Base.setindex!
#import Base.*

export BEDFile
export decompress_genotypes!, decompress_genotypes
export A_mul_B!, A_mul_B
export At_mul_B!, At_mul_B
#export *
export sumsq!, sumsq
export mean!, prec!
export maf
export getindex
export compress
export read_plink_data
export prednames

# constants used for decompression purposes
const ZERO8  = convert(Int8,0)
const ONE8   = convert(Int8,1)
const TWO8   = convert(Int8,2)
const THREE8 = convert(Int8,3)

# PLINK magic numbers
const MNUM1  = convert(Int8,108)
const MNUM2  = convert(Int8,27)

# old PLINK magic numbers
#const MNUM1  = convert(Int8, -17)
#const MNUM2  = convert(Int8,-85)

# typealias floating point operations
const Float = Union{Float32,Float64}
"""
This lookup table encodes the following PLINK format for genotypes:

- 00 is homozygous for allele 1
- 01 is heterozygous
- 10 is missing
- 11 is homozygous for allele 2

We will represent missing (NA) with NaN.
Further note that the bytes are read from right to left.
That is, if we label each of the 8 position as A to H, we would label backwards:

    01101100
    HGFEDCBA

and so the first four genotypes are read as follows:

    01101100
    HGFEDCBA

          AB   00  -> homozygote (first)        -> 0
        CD     11  -> other homozygote (second) -> 2
      EF       01  -> heterozygote (third)      -> 1
    GH         10  -> missing genotype (fourth) -> NaN

Finally, when we reach the end of a SNP (or if in individual-mode, the end of an individual),
then we skip to the start of a new byte (i.e. skip any remaining bits in that byte).
For a precise desceiption of PLINK BED files, see the file type reference in the [PLINK documentation](http://pngu.mgh.harvard.edu/~purcell/plink/binary.shtml).

The implementation here uses bitshifting and a bit threshold against `Int8` value 3 to interpret the compressed data.
The bitshifting trick works left-to-right, in contrast to the PLINK convention of reading right-to-left.
Thus, the map in PLINK.jl is slightly different:

- 00 -> 0 -> 0.0
- 11 -> 3 -> 2.0
- **10 -> 2 -> 1.0**
- **01 -> 1 -> NaN**
"""
const genofloat = [0.0, NaN, 1.0, 2.0]
const genoint   = [0 -99 1 2]

# dictionaries used in compression to convert floating point numbers to Int8 numbers
const bin32  = Dict{Float32, Int8}(0.0f0 => ZERO8, NaN32 => ONE8, 1.0f0 => TWO8, 2.0f0 => THREE8)
const bin64  = Dict{Float64, Int8}(0.0 => ZERO8, NaN => ONE8, 1.0 => TWO8, 2.0 => THREE8)

# preloaded GPU kernels
# these are merely long strongs that contain the code in the file
const gpucode64 = readstring(open(Pkg.dir() * "/PLINK/src/kernels/iht_kernels64.cl"))
const gpucode32 = readstring(open(Pkg.dir() * "/PLINK/src/kernels/iht_kernels32.cl"))

include("covariate.jl")
include("data.jl")
include("genomatrix.jl")
include("bedfile.jl")
include("compression.jl")
include("decompression.jl")
if cl != nothing
    include("gpu.jl")
end
include("linalg.jl")

end # end module PLINK

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
"""
    Reluplex(optimizer, eager::Bool)

Reluplex uses binary tree search to find an activation pattern that maps a feasible input to an infeasible output.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle
3. Output: halfspace

# Return
`CounterExampleResult`

# Method
Binary search of activations (0/1) and pruning by optimization.

# Property
Sound and complete.

# Reference
[G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer,
"Reluplex: An Efficient SMT Solver for Verifying Deep Neural Networks," in
*International Conference on Computer Aided Verification*, 2017.](https://arxiv.org/abs/1702.01135)
"""
struct Reluplex end

function solve(solver::Reluplex, problem::Problem)
    basic_model = new_model(solver)
    bs, fs = encode(solver, basic_model, problem)
    layers = problem.network.layers
    initial_status = [zeros(Int, n) for n in n_nodes.(layers)]

    return reluplex_step(solver, problem, basic_model, bs, fs, initial_status)
end

function find_relu_to_fix(b_vars, f_vars)
    for i in 1:length(f_vars), j in 1:length(f_vars[i])
        b = getvalue(b_vars[i+1][j])
        f = getvalue(f_vars[i][j])

        if type_one_broken(b, f) ||
           type_two_broken(b, f)
            return (i, j)
        end
    end
    return (0, 0)
end

type_one_broken(b::Real, f::Real) = (f >= 0.0) && (f != b)  # TODO consider renaming to `inactive_broken` and `active_broken`
type_two_broken(b::Real, f::Real) = (f == 0.0) && (b > 0.0)

# NOTE that the b that is passed in should be: bs[i+1] relative to fs[i]
function type_one_repair!(m::Model, b::Variable, f::Variable)
    @constraint(m, b == f)
    @constraint(m, b >= 0.0)
    return nothing
end
function type_two_repair!(m::Model, b::Variable, f::Variable)
    @constraint(m, b <= 0.0)
    @constraint(m, f == 0.0)
    return nothing
end

function encode(solver::Reluplex, model::Model,  problem::Problem)
    layers = problem.network.layers
    bs = init_neurons(model, layers)
    fs = init_forward_facing_vars(model, layers)

    # Positivity contraint for forward-facing variables
    for i in 1:length(fs)
        @constraint(model, fs[i] .>= 0.0)
    end

    # All layers (input and hidden) get an "input" constraint
    # In addition, each set of back-facing and forward-facing
    # variables are related to each other by a constraint.
    bounds = get_bounds(problem)
    for (i, L) in enumerate(layers)
        ## layerwise input constraint
        add_set_constraint!(model, bounds[i], bs[i])

        ## b<——>f[next] constraint
        # first layer technically has only vars which are forward facing,
        # but they behave like b-vars, so they are treated as such.
        W, bias = L.weights, L.bias
        if (i == 1)
            @constraint(model, -bs[i+1] .+ W*bs[i] .== -bias)
        else
            @constraint(model, -bs[i+1] .+ W*fs[i-1] .== -bias)
        end

        ## f >= b always, by definition
        if i <= length(fs)
            @constraint(model, fs[i] .>= bs[i+1])
        end
    end
    add_complementary_set_constraint!(model, problem.output, last(bs))

    zero_objective!(model)

    return bs, fs
end

function enforce_repairs!(model::Model, bs, fs, relu_status)
    # Need to decide what to do with last layer, this assumes there is no ReLU.
    for i in 1:(length(relu_status)-1), j in 1:length(relu_status[i])
        b = bs[i+1][j]
        f = fs[i][j]
        if relu_status[i][j] == 1
            type_one_repair!(model, b, f)
        elseif relu_status[i][j] == 2
            type_two_repair!(model, b, f)
        end
    end
end

function reluplex_step(solver::Reluplex,
                       problem::Problem,
                       model::Model,
                       b_vars::Vector{Vector{Variable}},
                       f_vars::Vector{Vector{Variable}},
                       relu_status::Vector{Vector{Int}})

    status = solve(model, suppress_warnings = true)
    if status == :Infeasible
        return CounterExampleResult(:SAT)

    elseif status == :Optimal
        i, j = find_relu_to_fix(b_vars, f_vars)

        # in case no broken relus could be found, return the "input" as a countereexample
        i == 0 && return CounterExampleResult(:UNSAT, getvalue.(first(b_vars)))

        for repair_type in 1:2
            relu_status[i][j] = repair_type

            new_m  = new_model(solver)
            bs, fs = encode(solver, new_m, problem)
            enforce_repairs!(new_m, bs, fs, relu_status)

            result = reluplex_step(solver, problem, new_m, bs, fs, relu_status)

            relu_status[i][j] = 0

            return result
        end
    else
        error("unexpected status $status") # are there alternatives to the if and elseif?
    end
end

# for convenience:
new_model(::Reluplex) = Model(solver = GLPKSolverLP(method = :Exact))

# doesn't do what it should! TODO open feature request issue on JuMP
# function extract_bs_fs(m::Model)
#     vars = collect(keys(m.varData))
#     L1, L2 = map(length, vars)
#     # "b_vars" is the longer of the two, and should be returned first
#     if L1 < L2
#         reverse!(vars)
#     end
#     return vars
# end

ext: jl


text:
# just some Julia for plotting the benchmarks and examples

using Plots
using DelimitedFiles
using Statistics
A = readdlm("areaData.txt",',');
T = readdlm("runtimeData.txt",',');

D = collect(0:size(A,1))
inds = [0,8*4 .^D...].+1
Areas = [Float64.(A[i,inds[i]:inds[i+1]-1]) for i in 1:size(D,1)-1];

histogram(Areas[11],normed=true,label="",title="Histogram of Trixel areas at depth 10")
savefig("Depth10Histogram.pdf")
p1 = scatter(Int.(collect(0:size(Areas,1)-1)),mean.(Areas),ribbon=std.(Areas),label="",title="Mean and std Trixel Area at Depth d",xlabel="Depth d", ylabel="Area")
p2 = scatter(Int.(collect(0:size(Areas,1)-1)),mean.(Areas),ribbon=std.(Areas),label="",title="Mean and std Trixel Area at Depth d",xlabel="Depth d",yaxis=:log10,ylabel="Log10 Area")
plot(p1,p2,size=(900,600))
savefig("MeanStdTrixelArea.pdf")
plot(collect(0:(size(T',1)-2)).+1,T[1,1:end-1],label="",xaxis=:log10,yaxis=:log10,title = "Runtime loglog",xlabel="depth+1",ylabel="time, s")
savefig("TimeComplexity.pdf")

mutable struct Trixel{T<:AbstractArray}
    vertices::T
    id::String
end
function Trixel(;vertices=zeros(2,3))
    if size(vertices,1) != 3
        println("Must have exactly 3 vertices")
        return Triangle(zeros(2,3),"")
    end
    return Trixel(vertices,"")
end
function DrawTrixel!(T::Trixel;subplot=nothing)
    if subplot == nothing
        plot!([T.vertices[1,1],T.vertices[1,2]],
              [T.vertices[2,1],T.vertices[2,2]],
              [T.vertices[3,1],T.vertices[3,2]],label="")
        plot!([T.vertices[1,1],T.vertices[1,3]],
          [T.vertices[2,1],T.vertices[2,3]],
          [T.vertices[3,1],T.vertices[3,3]],label="")
        plot!([T.vertices[1,2],T.vertices[1,3]],
          [T.vertices[2,2],T.vertices[2,3]],
          [T.vertices[3,2],T.vertices[3,3]],label="")
    else
        plot!([T.vertices[1,1],T.vertices[1,2]],
              [T.vertices[2,1],T.vertices[2,2]],
              [T.vertices[3,1],T.vertices[3,2]],label="",subplot=subplot)
        plot!([T.vertices[1,1],T.vertices[1,3]],
          [T.vertices[2,1],T.vertices[2,3]],
          [T.vertices[3,1],T.vertices[3,3]],label="",subplot=subplot)
        plot!([T.vertices[1,2],T.vertices[1,3]],
          [T.vertices[2,2],T.vertices[2,3]],
          [T.vertices[3,2],T.vertices[3,3]],label="",subplot=subplot)
    end
end

p = plot(layout=5,aspect_ratio=:equal)
for i in 0:4
    HTM = readdlm("HTM-$i.txt",',')
    T = []
    for j in 1:size(HTM,1)
        v = Float64.(HTM[j,1:9])
        trix = Trixel(vertices=reshape(v,(3,3)))
        push!(T,trix)
        DrawTrixel!(trix,subplot=i+1)
    end
    plot!(title="HTM: depth = $i",subplot=i+1)
end
plot!(xaxis=nothing,yaxis=nothing,zaxis=nothing)
savefig(p,"ExampleHTMs.pdf")
savefig(p,"ExampleHTMs.png")

max_stars_repo_name: harveydevereux/HTM

max_issues_repo_name: harveydevereux/HTM

max_forks_repo_name: harveydevereux/HTM


text:
export AbstractMarkovDecisionProcess, MarkovDecisionProcess,
        reward, transition_model, actions,
        GridMarkovDecisionProcess, go_to, show_grid, to_arrows,
        value_iteration, expected_utility, optimal_policy,
        policy_evaluation, policy_iteration;

abstract type AbstractMarkovDecisionProcess end;

#=

    MarkovDecisionProcess is a MDP implementation of AbstractMarkovDecisionProcess.

    A Markov decision process is a sequential decision problem with fully observable

    and stochastic environment with a transition model and rewards function.

    The discount factor (gamma variable) describes the preference for current rewards

    over future rewards.

=#
struct MarkovDecisionProcess{T} <: AbstractMarkovDecisionProcess
    initial::T
    states::Set{T}
    actions::Set{T}
    terminal_states::Set{T}
    transitions::Dict
    gamma::Float64
    reward::Dict

    function MarkovDecisionProcess{T}(initial::T, actions_list::Set{T}, terminal_states::Set{T}, transitions::Dict, states::Union{Nothing, Set{T}}, gamma::Float64) where T
        if (!(0 < gamma <= 1))
            error("MarkovDecisionProcess(): The gamma variable of an MDP must be between 0 and 1, the constructor was given ", gamma, "!");
        end
        local new_states::Set{typeof(initial)};
        if (typeof(states) <: Set)
            new_states = states;
        else
            new_states = Set{typeof(initial)}();
        end
        return new(initial, new_states, actions_list, terminal_states, transitions, gamma, Dict());
    end
end

MarkovDecisionProcess(initial, actions_list::Set, terminal_states::Set, transitions::Dict; states::Union{Nothing, Set}=nothing, gamma::Float64=0.9) = MarkovDecisionProcess{typeof(initial)}(initial, actions_list, terminal_states, transitions, states, gamma);

"""
    reward(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}

Return a reward based on the given 'state'.
"""
function reward(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}
    return mdp.reward[state];
end

"""
    transition_model(mdp::T, state, action) where {T <: AbstractMarkovDecisionProcess}

Return a list of (P(s'|s, a), s') pairs given the state 's' and action 'a'.
"""
function transition_model(mdp::T, state, action) where {T <: AbstractMarkovDecisionProcess}
    if (length(mdp.transitions) == 0)
        error("transition_model(): The transition model for the given 'mdp' could not be found!");
    else
        return mdp.transitions[state][action];
    end
end

"""
    actions(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}

Return a set of actions that are possible in the given state.
"""
function actions(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}
    if (state in mdp.terminal_states)
        return Set{Nothing}([nothing]);
    else
        return mdp.actions;
    end
end

#=

    GridMarkovDecisionProcess is a two-dimensional environment MDP implementation

    of AbstractMarkovDecisionProcess. Obstacles in the environment are represented

    by a null.

=#
struct GridMarkovDecisionProcess <: AbstractMarkovDecisionProcess
    initial::Tuple{Int64, Int64}
    states::Set{Tuple{Int64, Int64}}
    actions::Set{Tuple{Int64, Int64}}
    terminal_states::Set{Tuple{Int64, Int64}}
    grid::Array{Union{Nothing, Float64}, 2}
    gamma::Float64
    reward::Dict

    function GridMarkovDecisionProcess(initial::Tuple{Int64, Int64}, terminal_states::Set{Tuple{Int64, Int64}}, grid::Array{Union{Nothing, Float64}, 2}; states::Union{Nothing, Set{Tuple{Int64, Int64}}}=nothing, gamma::Float64=0.9)
        if (!(0 < gamma <= 1))
            error("GridMarkovDecisionProcess(): The gamma variable of an MDP must be between 0 and 1, the constructor was given ", gamma, "!");
        end
        local new_states::Set{Tuple{Int64, Int64}};
        if (typeof(states) <: Set)
            new_states = states;
        else
            new_states = Set{Tuple{Int64, Int64}}();
        end
        local orientations::Set = Set{Tuple{Int64, Int64}}([(1, 0), (0, 1), (-1, 0), (0, -1)]);
        local reward::Dict = Dict();
        for i in 1:getindex(size(grid), 1)
            for j in 1:getindex(size(grid, 2))
                reward[(i, j)] = grid[i, j]
                if (!(grid[i, j] === nothing))
                    push!(new_states, (i, j));
                end
            end
        end
        return new(initial, new_states, orientations, terminal_states, grid, gamma, reward);
    end
end

"""
    go_to(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, direction::Tuple{Int64, Int64})

Return the next state given the current state and direction.
"""
function go_to(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, direction::Tuple{Int64, Int64})
    local next_state::Tuple{Int64, Int64} = map(+, state, direction);
    if (next_state in gmdp.states)
        return next_state;
    else
        return state;
    end
end

function transition_model(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, action::Nothing)
    return [(0.0, state)];
end

function transition_model(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, action::Tuple{Int64, Int64})
    return [(0.8, go_to(gmdp, state, action)),
            (0.1, go_to(gmdp, state, utils.turn_heading(action, -1))),
            (0.1, go_to(gmdp, state, utils.turn_heading(action, 1)))];
end

function show_grid(gmdp::GridMarkovDecisionProcess, mapping::Dict)
    local grid::Array{Union{Nothing, String}, 2};
    local rows::AbstractVector = [];
    for i in 1:getindex(size(gmdp.grid), 1)
        local row::Array{Union{Nothing, String}, 1} = Array{Union{Nothing, String}, 1}();
        for j in 1:getindex(size(gmdp.grid), 2)
            push!(row, get(mapping, (i, j), nothing));
        end
        push!(rows, reshape(row, (1, length(row))));
    end
    grid = reduce(vcat, rows);
    return grid;
end

# (0, 1) will move the agent rightward.
# (-1, 0) will move the agent upward.
# (0, -1) will move the agent leftward.
# (1, 0) will move the agent downward.
function to_arrows(gmdp::GridMarkovDecisionProcess, policy::Dict)
    local arrow_characters::Dict = Dict([Pair((0, 1), ">"),
                                        Pair((-1, 0), "^"),
                                        Pair((0, -1), "<"),
                                        Pair((1, 0), "v"),
                                        Pair(nothing, ".")]);
    return show_grid(gmdp, Dict(collect(Pair(state, arrow_characters[action])
                                    for (state, action) in policy)));
end

# An example sequential decision problem (Fig. 17.1a) where an agent does not
# terminate until it reaches a terminal state in the 4x3 environment (Fig. 17.1a).
#
# Matrices in Julia start from the upper-left corner and index (1, 1).
sequential_decision_environment = GridMarkovDecisionProcess((1, 1),
                                            Set([(2, 4), (3, 4)]),
                                            [-0.04 -0.04 -0.04 -0.04;
                                            -0.04 nothing -0.04 -1;
                                            -0.04 -0.04 -0.04 +1]);

"""
    value_iteration(mdp::T; epsilon::Float64=0.001) where {T <: AbstractMarkovDecisionProcess}

Return the utilities of the MDP's states as a Dict by applying the value iteration algorithm (Fig. 17.4)
on the given Markov decision process 'mdp' and a arbitarily small positive number 'epsilon'.
"""
function value_iteration(mdp::T; epsilon::Float64=0.001) where {T <: AbstractMarkovDecisionProcess}
    local U_prime::Dict = Dict(collect(Pair(state, 0.0) for state in mdp.states));
    while (true)
        local U::Dict = copy(U_prime);
        local delta::Float64 = 0.0;
        for state in mdp.states
            U_prime[state] = (reward(mdp, state)
                                + (mdp.gamma
                                * max((sum(collect(p * U[state_prime]
                                                    for (p, state_prime) in transition_model(mdp, state, action)))
                                        for action in actions(mdp, state))...)));
            delta = max(delta, abs(U_prime[state] - U[state]));
        end
        if (delta < ((epsilon * (1 - mdp.gamma))/mdp.gamma))
            return U;
        end
    end
end

function value_iteration(gmdp::GridMarkovDecisionProcess; epsilon::Float64=0.001)
    local U_prime::Dict = Dict(collect(Pair(state, 0.0) for state in gmdp.states));
    while (true)
        local U::Dict = copy(U_prime);
        local delta::Float64 = 0.0;
        for state in gmdp.states
            U_prime[state] = (reward(gmdp, state)
                            + (gmdp.gamma
                            * max((sum(collect(p * U[state_prime]
                                                for (p, state_prime) in transition_model(gmdp, state, action)))
                                        for action in actions(gmdp, state))...)));
            delta = max(delta, abs(U_prime[state] - U[state]));
        end
        if (delta < ((epsilon * (1 - gmdp.gamma))/gmdp.gamma))
            return U;
        end
    end
end

function expected_utility(mdp::T, U::Dict, state::Tuple{Int64, Int64}, action::Tuple{Int64, Int64}) where {T <: AbstractMarkovDecisionProcess}
    return sum((p * U[state_prime] for (p, state_prime) in transition_model(mdp, state, action)));
end

function expected_utility(mdp::T, U::Dict, state::Tuple{Int64, Int64}, action::Nothing) where {T <: AbstractMarkovDecisionProcess}
    return sum((p * U[state_prime] for (p, state_prime) in transition_model(mdp, state, action)));
end

"""
    optimal_policy(mdp::T, U::Dict) where {T <: AbstractMarkovDecisionProcess}

Return the optimal_policy 'π*(s)' (Equation 17.4) given the Markov decision process 'mdp'
and the utility function 'U'.
"""
function optimal_policy(mdp::T, U::Dict) where {T <: AbstractMarkovDecisionProcess}
    local pi::Dict = Dict();
    for state in mdp.states
        pi[state] = argmax(collect(actions(mdp, state)), (function(action::Union{Nothing, Tuple{Int64, Int64}})
                                                                return expected_utility(mdp, U, state, action);
                                                            end));
    end
    return pi;
end

"""
    policy_evaluation(pi::Dict, U::Dict, mdp::T; k::Int64=20) where {T <: AbstractMarkovDecisionProcess}

Return the updated utilities of the MDP's states by applying the modified policy iteration
algorithm on the given Markov decision process 'mdp', utility function 'U', policy 'pi',
and number of Bellman updates to use 'k'.
"""
function policy_evaluation(pi::Dict, U::Dict, mdp::T; k::Int64=20) where {T <: AbstractMarkovDecisionProcess}
    for i in 1:k
        for state in mdp.states
            U[state] = (reward(mdp, state)
                        + (mdp.gamma
                        * sum((p * U[state_prime] for (p, state_prime) in transition_model(mdp, state, pi[state])))));
        end
    end
    return U;
end

function policy_evaluation(pi::Dict, U::Dict, gmdp::GridMarkovDecisionProcess; k::Int64=20)
    for i in 1:k
        for state in gmdp.states
            U[state] = (reward(gmdp, state)
                        + (gmdp.gamma
                        * sum((p * U[state_prime] for (p, state_prime) in transition_model(gmdp, state, pi[state])))));
        end
    end
    return U;
end

"""
    policy_iteration(mdp::T) where {T <: AbstractMarkovDecisionProcess}

Return a policy using the policy iteration algorithm (Fig. 17.7) given the Markov decision process 'mdp'.
"""
function policy_iteration(mdp::T) where {T <: AbstractMarkovDecisionProcess}
    local U::Dict = Dict(collect(Pair(state, 0.0) for state in mdp.states));
    local pi::Dict = Dict(collect(Pair(state, rand(RandomDeviceInstance, collect(actions(mdp, state))))
                                    for state in mdp.states));
    while (true)
        U = policy_evaluation(pi, U, mdp);
        local unchanged::Bool = true;
        for state in mdp.states
            local action = argmax(collect(actions(mdp, state)), (function(action::Union{Nothing, Tuple{Int64, Int64}})
                                                                    return expected_utility(mdp, U, state, action);
                                                                end));
            if (action != pi[state])
                pi[state] = action;
                unchanged = false;
            end
        end
        if (unchanged)
            return pi;
        end
    end
end

max_stars_repo_name: mikhail-j/aima-julia

max_issues_repo_name: mikhail-j/aima-julia

max_forks_repo_name: mikhail-j/aima-julia


text:
# Performance of global constant vs. global variable
#
# Concept:
# Global variables hurts performance because compiler cannot generate
# efficient code.  Make it constant helps the compiler tremendously.

using BenchmarkTools

# using global variable
variable = 10
function add_using_global_variable(x)
    return x + variable
end

@btime add_using_global_variable(10);
#=
julia> @btime add_using_global_variable(10);
  25.572 ns (0 allocations: 0 bytes)
=#

function add_using_function_arg(x, y)
    return x + y
end

@btime add_using_function_arg(10, 10);
#=
julia> @btime add_using_function_arg(10, 10)
  0.030 ns (0 allocations: 0 bytes)
=#

@code_llvm add_using_function_arg(10, 10)
#=
julia> @code_llvm add_using_function_arg(10, 10)

;  @ REPL[158]:6 within `add_using_function_arg'
define i64 @julia_add_using_function_arg_15028(i64, i64) {
top:
; ┌ @ int.jl:53 within `+'
   %2 = add i64 %1, %0
; └
  ret i64 %2
}
=#

@code_llvm add_using_global_variable(10)
#=
julia> @code_llvm add_using_global_variable(10)

;  @ REPL[156]:2 within `add_using_global_variable'
define nonnull %jl_value_t addrspace(10)* @julia_add_using_global_variable_14990(i64) {
top:
  %1 = alloca %jl_value_t addrspace(10)*, i32 3
  %gcframe = alloca %jl_value_t addrspace(10)*, i32 4
  %2 = bitcast %jl_value_t addrspace(10)** %gcframe to i8*
  call void @llvm.memset.p0i8.i32(i8* %2, i8 0, i32 32, i32 0, i1 false)
  %3 = call %jl_value_t*** inttoptr (i64 4493620896 to %jl_value_t*** ()*)() #2
  %4 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %gcframe, i32 0
  %5 = bitcast %jl_value_t addrspace(10)** %4 to i64*
  store i64 4, i64* %5
  %6 = getelementptr %jl_value_t**, %jl_value_t*** %3, i32 0
  %7 = load %jl_value_t**, %jl_value_t*** %6
  %8 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %gcframe, i32 1
  %9 = bitcast %jl_value_t addrspace(10)** %8 to %jl_value_t***
  store %jl_value_t** %7, %jl_value_t*** %9
  %10 = bitcast %jl_value_t*** %6 to %jl_value_t addrspace(10)***
  store %jl_value_t addrspace(10)** %gcframe, %jl_value_t addrspace(10)*** %10
  %11 = load %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** inttoptr (i64 4757987176 to %jl_value_t addrspace(10)**), align 8
  %12 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %gcframe, i32 2
  store %jl_value_t addrspace(10)* %11, %jl_value_t addrspace(10)** %12
  %13 = call %jl_value_t addrspace(10)* @jl_box_int64(i64 signext %0)
  %14 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %gcframe, i32 3
  store %jl_value_t addrspace(10)* %13, %jl_value_t addrspace(10)** %14
  %15 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %1, i32 0
  store %jl_value_t addrspace(10)* addrspacecast (%jl_value_t* inttoptr (i64 4582122192 to %jl_value_t*) to %jl_value_t addrspace(10)*), %jl_value_t addrspace(10)** %15
  %16 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %1, i32 1
  store %jl_value_t addrspace(10)* %13, %jl_value_t addrspace(10)** %16
  %17 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %1, i32 2
  store %jl_value_t addrspace(10)* %11, %jl_value_t addrspace(10)** %17
  %18 = call nonnull %jl_value_t addrspace(10)* @jl_apply_generic(%jl_value_t addrspace(10)** %1, i32 3)
  %19 = getelementptr %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %gcframe, i32 1
  %20 = load %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)** %19
  %21 = getelementptr %jl_value_t**, %jl_value_t*** %3, i32 0
  %22 = bitcast %jl_value_t*** %21 to %jl_value_t addrspace(10)**
  store %jl_value_t addrspace(10)* %20, %jl_value_t addrspace(10)** %22
  ret %jl_value_t addrspace(10)* %18
}
=#

const constant = 10
function add_using_global_constant(x)
    return constant + x
end

@btime add_using_global_constant(10);
#=
julia> @btime add_using_global_constant(10);
  0.030 ns (0 allocations: 0 bytes)
=#

@code_llvm add_using_global_constant(10)
#=
julia> @code_llvm add_using_global_constant(10)

;  @ REPL[166]:2 within `add_using_global_constant'
define i64 @julia_add_using_global_constant_15035(i64) {
top:
; ┌ @ int.jl:53 within `+'
   %1 = add i64 %0, 10
; └
  ret i64 %1
}
=#

# ----------------------------------------------------------------------
# what if variable is annotated with type?
# ----------------------------------------------------------------------
variable::Int = 10

#=
julia> variable::Int = 10
ERROR: syntax: type declarations on global variables are not yet supported
=#

function add_using_global_variable_typed(x)
    return x + variable::Int
end

variable = 10
@btime add_using_global_variable_typed(10);
#=
julia> @btime add_using_global_variable_typed(10)
  5.213 ns (0 allocations: 0 bytes)
=#

variable = 10.0
add_using_global_variable_typed(10)
#= what if the variable got assigned to a value with different type?
julia> variable = 10.0
10.0

julia> add_using_global_variable_typed(10)
ERROR: TypeError: in typeassert, expected Int64, got Float64
=#


text:
using TimeParallel
using Test

@testset "ErrorControl" begin
    # Write your tests here.
end

@testset "Coarse" begin
    # Write your tests here.
end

@testset "Parareal" begin
    # Write your tests here.
end


text:
# Methods for extending Distributions.jl interface to InfExtendedReal and Intervals

_pdf(dbn, x) = pdf(dbn, x)
_logpdf(dbn, x) = logpdf(dbn, x)
_cdf(dbn, x) = cdf(dbn, x)
_ccdf(dbn, x) = ccdf(dbn, x)

# Extend to InfExtendedReal

function _pdf(dbn, x::InfExtendedReal)
    isinf(x) ? 0.0 : pdf(dbn, x.val)
end
function _logpdf(dbn, x::InfExtendedReal)
    isinf(x) ? -Inf : logpdf(dbn, x.val)
end
function _cdf(dbn, x::InfExtendedReal)
    isposinf(x) ? 1.0 : isneginf(x) ? 0.0 : cdf(dbn, x.val)
end
function _ccdf(dbn, x::InfExtendedReal)
    isneginf(x) ? 1.0 : isposinf(x) ? 0.0 : ccdf(dbn, x.val)
end

function cdf(dbn, int::Union{<:ClosedClosed, <:OpenClosed})
    _cdf(dbn, rightendpoint(int))
end
function cdf(dbn::ContinuousUnivariateDistribution, int::Union{<:ClosedOpen, <:OpenOpen})
    _cdf(dbn, rightendpoint(int))
end
function cdf(dbn::DiscreteDistribution, int::Union{<:ClosedOpen, <:OpenOpen})
    b = rightendpoint(int)
    _cdf(dbn, b) - _pdf(dbn, b)
end

function pdf(dbn::ContinuousUnivariateDistribution, int::OpenClosed)
    a, b = endpoints(int)
    isposinf(b) ? _ccdf(dbn, a) : isneginf(a) ? _cdf(dbn, b) : _cdf(dbn, b) - _cdf(dbn, a)
end

function pdf(dbn::ContinuousUnivariateDistribution, int::AbstractInterval)
    a, b = endpoints(int)
    pdf(dbn, Interval{:open,:closed}(a,b))
end

function pdf(dbn::DiscreteDistribution, int::OpenOpen)
    a, b = endpoints(int)
    pdf(dbn, Interval{:open,:closed}(a,b)) - _pdf(dbn, b)
end

function pdf(dbn::DiscreteDistribution, int::ClosedOpen)
    a, b = endpoints(int)
    pdf(dbn, Interval{:open,:closed}(a,b)) - _pdf(dbn, b) + _pdf(dbn, a)
end

#logcdf, logccdf, logdiffcdf


text:
# recipe.jl

using RecipesBase

@recipe function g(f::TaylorModel1)
    ffp = fp_rpa(f)
    fT = polynomial(ffp)
    Δ = remainder(ffp)
    ξ0 = ffp.x0

    alpha --> 0.3
    seriestype := :shape

    xs = range(f.dom.lo, stop=f.dom.hi, length=100)
    evals = fT.(xs .- ξ0) .+ Δ

    # make polygon:
    xs = [xs; reverse(xs); xs[1]]
    ys = [inf.(evals); reverse(sup.(evals)); inf(evals[1])]

    xs, ys
end

@recipe function g(f::RTaylorModel1)
    ffp = fp_rpa(f)
    fT = polynomial(ffp)
    Δ = remainder(ffp)
    ξ0 = ffp.x0
    order = get_order(f)+1

    alpha --> 0.5
    seriestype := :shape

    xs = range(f.dom.lo, stop=f.dom.hi, length=100)
    evals = fT.(xs .- ξ0)

    corrs = (xs .- ξ0) .^ order
    Δrel = Δ .* corrs
    evalslo = inf.(evals + Δrel)
    evalshi = sup.(evals + Δrel)

    xs = [xs; reverse(xs); xs[1]]
    ys = [evalslo; reverse(evalshi); evalslo[1]]

    xs, ys
end


text:
import Base: hcat, vcat, \
export hcat, vcat, \, diagOp, A_mul_B

include("NFFTOp.jl")
include("ExplicitOp.jl")
include("SensitivityOp.jl")
include("SamplingOp.jl")
include("MapSliceOp.jl")
include("FieldmapNFFTOp.jl")
include("EncodingOp.jl")
include("SparseOp.jl")

"""
    hcat(A::AbstractLinearOperator, n::Int)

horizontally concatenates alinear operator `A` n times
"""
function hcat(A::AbstractLinearOperator, n::Int)
  op = A
  for i = 2:n
    op = LinearOperators.hcat(op, A)
  end
  return op
end

"""
    hcat(A::AbstractLinearOperator, n::Int)

vertically concatenates alinear operator `A` n times
"""
function vcat(A::AbstractLinearOperator, n::Int)
  op = A
  for i = 2:n
    op = LinearOperators.vcat(op, A)
  end
  return op
end

function diagOpProd(x::Vector{T}, nrow::Int, ops :: AbstractLinearOperator...) where T
  y = Vector{T}(undef, nrow)
  xIdx=1
  yIdx=1
  for i=1:length(ops)
    y[yIdx:yIdx+ops[i].nrow-1] = ops[i]*x[xIdx:xIdx+ops[i].ncol-1]
    xIdx += ops[i].ncol
    yIdx += ops[i].nrow
  end
  return y
end

function diagOpTProd(x::Vector{T}, ncol::Int, ops :: AbstractLinearOperator...) where T
  y = Vector{T}(undef, ncol)
  xIdx=1
  yIdx=1
  for i=1:length(ops)
    y[yIdx:yIdx+ops[i].ncol-1] = transpose(ops[i])*x[xIdx:xIdx+ops[i].nrow-1]
    xIdx += ops[i].nrow
    yIdx += ops[i].ncol
  end
  return y
end

function diagOpCTProd(x::Vector{T}, ncol::Int, ops :: AbstractLinearOperator...) where T
  y = Vector{T}(undef, ncol)
  xIdx=1
  yIdx=1
  for i=1:length(ops)
    y[yIdx:yIdx+ops[i].ncol-1] = adjoint(ops[i])*x[xIdx:xIdx+ops[i].nrow-1]
    xIdx += ops[i].nrow
    yIdx += ops[i].ncol
  end
  return y
end

"""
    diagOp(ops :: AbstractLinearOperator...)

create a bloc-diagonal operator out of the `LinearOperator`s contained in ops
"""
function diagOp(ops :: AbstractLinearOperator...)
  nrow=0
  ncol=0
  S = eltype(ops[1])
  for i = 1:length(ops)
    nrow += ops[i].nrow
    ncol += ops[i].ncol
    S = promote_type(S, eltype(ops[i]))
  end

  return LinearOperator{S}( nrow, ncol, false, false
                            , x->diagOpProd(x,nrow,ops...)
                            , y->diagOpTProd(y,ncol,ops...)
                            , y->diagOpCTProd(y,ncol,ops...) )
end

# implement A_mul_B for the product
A_mul_B(A::AbstractLinearOperator, x::Vector) = A*x

#
# use hermitian conjugate as an estimate for the inverse (fallback)
#
\(A::AbstractLinearOperator, x::Vector) = adjoint(A)*x


text:
#unum-hlayer.jl - human layer things in the unum library.

#modify the show() directive so that a text display of both unum types outputs
#as a "Unum" object and hides the underlying "UnumSmall"/"UnumLarge" distinction.
#
#N.B. typeof() will correctly identify the Unums.UnumSmall and Unums.UnumLarge
#types.

@universal function Base.show(io::IO, x::Unum)

  @typenames

  #we want to be able to represent having g-layer flags as part of this.
  gflags = (x.flags & (~UNUM_FLAG_MASK))
  #for nan, let's also show the noisy nan bit.
  isnan(x) && (gflags |= (x.flags & UNUM_SIGN_MASK))
  gflagstring = (gflags == 0) ? "" : @sprintf ", 0x%04X" gflags
  is_pos_inf(x) && (print(io, "inf($uname{$ESS,$FSS}$gflagstring)"); return)
  is_pos_mmr(x) && (print(io, "mmr($uname{$ESS,$FSS}$gflagstring)"); return)
  is_pos_sss(x) && (print(io, "sss($uname{$ESS,$FSS}$gflagstring)"); return)
  is_neg_inf(x) && (print(io, "-inf($uname{$ESS,$FSS}$gflagstring)"); return)
  is_neg_mmr(x) && (print(io, "-mmr($uname{$ESS,$FSS}$gflagstring)"); return)
  is_neg_sss(x) && (print(io, "-sss($uname{$ESS,$FSS}$gflagstring)"); return)
  is_zero(x)    && (print(io, "zero($uname{$ESS,$FSS})"); return)
  isnan(x)      && (print(io,"nan($uname{$ESS,$FSS}$gflagstring)");return)

  fsize_string = @sprintf "0x%04X" x.fsize
  esize_string = @sprintf "0x%04X" x.esize
  flags_string = @sprintf "0x%04X" x.flags

  (FSS < 7) ? (fraction_string = @sprintf "0x%016X" x.fraction) : fraction_string = string(x.fraction.a)
  exponent_string = @sprintf "0x%016X" x.exponent
  print(io, "$uname{$ESS,$FSS}($exponent_string, $fraction_string, $flags_string, $esize_string, $fsize_string)")
end

@universal function Base.show(io::IO, T::Type{Unum})
  #strips the Large or Small suffix when displaying this type.
  @typenames
  print(io, "$uname{$ESS,$FSS}")
end
function Base.show{ESS,FSS}(io::IO, T::Type{Unum{ESS,FSS}})
  @typenames
  print(io, "$uname{$ESS,$FSS}")
end
function Base.show(io::IO, T::Type{Unum})
  print(io, "Unum")
end

@universal function Base.bits(x::Unum, space::String = "")
  res = string((x.flags & 0b10) >> 1,       space,
         bits(x.exponent)[64 - x.esize:64], space,
         bits(x.fraction)[1:(x.fsize + 1)], space,
         x.flags & 0b1)
  ESS > 0 && (res = string(res, space, bits(x.esize)[17-ESS:16]))
  FSS > 0 && (res = string(res, space, bits(x.fsize)[17-FSS:16]))
  res
end

___textual{ESS,FSS}(v::UnumSmall{ESS,FSS}) = (ESS > 3) ? calculate(v) : Float64(v)
___textual{ESS}(v::UnumSmall{ESS,6}) = calculate(v)
___textual{ESS,FSS}(v::UnumLarge{ESS,FSS}) = calculate(v)

@universal function __textual(v::Unum)
  is_neg_inf(v) && return "-∞"
  is_pos_inf(v) && return "∞"
  ___textual(v)
end

@universal function describe(v::Unum)
  if isnan(v)
    print("NaN")
  elseif is_exact(v)
    print(__textual(v))
    print(" ex")
  else
    print(__textual(glb(v)))
    print(" op → ")
    print(__textual(lub(v)))
    print(" op")
  end
  println()
end

@universal function describe(v::Ubound)
  print("Ubound{$ESS,$FSS}(")
  if is_exact(v.lower)
    print(__textual(v.lower))
    print(" ex")
  else
    print(__textual(glb(v.lower)))
    print(" op")
  end
  print(" → ")
  if is_exact(v.upper)
    print(__textual(v.upper))
    print(" ex")
  else
    print(__textual(lub(v.upper)))
    print(" op")
  end
  print(")")
  println()
end
export describe

doc"""
  `ℜ` is used to generate special sets of real numbers in prettyprint directives.

  * `Unum{4,6}(ℜ) == Ubound(neg_mmr(Unum{4,6}), pos_mmr(Unum{4,6}))`
  * `Unum{4,6}(ℜ(∘)) == Ubound(neg_inf(Unum{4,6}), pos_inf(Unum{4,6}))`
"""
type ℜ; end

doc"""
  `ℜ⁺` is used to generate special sets of real numbers in prettyprint directives.

  * `Unum{4,6}(ℜ⁺) == Ubound(zero(Unum{4,6}), pos_mmr(Unum{4,6}))`
  * `Unum{4,6}(ℜ⁺(*)) == Ubound(pos_sss(Unum{4,6}), pos_mmr(Unum{4,6}))`
  * `Unum{4,6}(ℜ⁺(∘)) == Ubound(zero(Unum{4,6}), pos_inf(Unum{4,6}))`
"""
type ℜ⁺; end  #positive real numbers

doc"""
  `ℜ⁻` is used to generate special sets of real numbers in prettyprint directives.

  * `Unum{4,6}(ℜ⁻) == Ubound(neg_mmr(Unum{4,6}), zero(Unum{4,6}))`
  * `Unum{4,6}(ℜ⁻(*)) == Ubound(neg_mmr(Unum{4,6}), neg_sss(Unum{4,6}))`
  * `Unum{4,6}(ℜ⁻(∘)) == Ubound(neg_inf(Unum{4,6}), zero(Unum{4,6}))`
"""
type ℜ⁻; end
type ∘; end

(::Type{ℜ})(::Type{∘}) = :_rextended
(::Type{ℜ⁺})(f::Function) = (f == (*)) ? :_rposstar : nothing
(::Type{ℜ⁻})(f::Function) = (f == (*)) ? :_rnegstar : nothing
(::Type{ℜ⁺})(::Type{∘}) = :_rposext
(::Type{ℜ⁻})(::Type{∘}) = :_rnegext

export ℜ, ℜ⁺, ℜ⁻, ∘

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
module AssigmentProblem

using JuMP, GLPK

export solve_assigm_prob

include("solving_AP.jl")

end # module


text:
#EXERCISE MADE BY ERIK MARTÍN GARZÓN & DAVID BOCES OBIS
# the packages we need
using Gtk, Graphics, Logging, Printf

include("affin_transformation.jl")
#File opening and reading matrix and vor together

# the main window
win = GtkWindow("SO(3)")

function init_canvas(h,w)
    # create the drawing canvas
    c = @GtkCanvas(h,w)

    # create the initial drawing inside a try/catch loop
    @guarded draw(c) do widget
        #  draw the background with the canvas drawing context
        # the code for this comes later
        draw_the_canvas(c)
    end

    show(c)
    return c
end

# make the canvas
the_canvas = init_canvas(1500,750)

# an array to store the entry boxes
entry_list = []

# an array of labels that we use to display normalized inputs,
# and which also gets modified from the callback
normalized_labels = []

function find_by_name(list, name)
    for item in list
        if get_gtk_property(item, :name, String) == name
            return item
        end
    end
    @warn name, "not found in list"
    @warn "available names are"
    for item in list
        @warn get_gtk_property(item, :name, String)
    end
end

function bold_label(label_string)
    label = GtkLabel("")
    GAccessor.markup(label, """<b>""" * label_string * """</b>""")
    return label
end

# Now put everything into the window,
# including the canvas
function init_window(win, canvas)
    # make a box for the drawing canvas
    canvas_box = GtkBox(:v)
    push!(canvas_box, bold_label("EJERCICIO 4.4:                                                                                EJERCICIO 4.5:"))
    push!(canvas_box, bold_label("                  Dibuja en 3d la escena, en las coordenadas globales.                                                                              Dibuja en 3d la misma escena, pero esta vez en las coordenadas de la camara."))
    push!(canvas_box, canvas)

    # make a containing box that will stack the widgets and the canvas side by side
    global_box = GtkBox(:h)
    push!(global_box, GtkLabel("   ")) # a very basic separator
    push!(global_box, canvas_box)

    # put it all inside the window
    push!(win, global_box)
end

# the background drawing
function draw_the_canvas(canvas)
    h   = height(canvas)
    w   =  width(canvas)
    ctx =  getgc(canvas)

    # clear the canvas
    rectangle(ctx, 750, 0, w, h)
    set_source_rgb(ctx, 1, 1, 1)
    fill(ctx)

    # Paint red rectangle
    # rectangle(ctx, 3w/4, 0, w/4, h/4)
    # set_source_rgb(ctx, 1, 0, 0)
    # fill(ctx)

    # Paint blue rectangle
    # rectangle(ctx, 0, 3h/4, w/4, h/4)
    # set_source_rgb(ctx, 0, 0, 1)
    # fill(ctx)

    # Paint an example line
    # set_line_width(ctx, 5)
    # set_source_rgb(ctx, 1, 0, 1)
    # move_to(ctx, 10, 10)
    # line_to(ctx, 100, 50)
    # stroke(ctx)

    # ------------------------------------------------------------------ DATA
    originX = 200
    originY = 500

    #AXIS R3
    X = [1;0;0]
    Y = [0;1;0]
    Z = [0;0;1]

    #WE PROJECT THE ROTATED AXIS IN R3 TO R2
    X2 = to_2d(X)
    Y2 = to_2d(Y)
    Z2 = to_2d(Z)

    cam = [2.2837;5.0798;2.2328]

    angle = -150 #degree

    u = [0.01;-0.2;1.0] #camera's rotation AXIS

    #Setting camera coords
    camCoord2d = to_2d(cam)*100

    #Camera axis rotations
    camXR = rotate_phi_z(angle, u, X)
    camYR = rotate_phi_z(angle, u, Y)
    camZR = rotate_phi_z(angle, u, Z)

    camX2d = to_2d(camXR)
    camY2d = to_2d(camYR)
    camZ2d = to_2d(camZR)

    # n = hcat(camXR,camYR,camZR)
    # println(det(n))

    #Setting up variables we will need for calculations
    C = hcat(camXR, camYR, camZR)
    D = hcat(C, cam)
    B = vcat(D, [0 0 0 1])
    Binv = B^-1

    focus = 1

    for i=1:4
        println(B[i,:])
    end

    a = [0.9115;1.9397;3.3304]
    b = [3.7207;2.8794;4.4372]
    c = [1.9659;1.0000;3.2588]
    d = [2.6663;3.8191;4.5087]

    a2d = to_2d(a)
    b2d = to_2d(b)
    c2d = to_2d(c)
    d2d = to_2d(d)

    # --------------------------------------------------------------- EXERCISE 4.3 / 4.5
    focal_camera(Binv, focus, a)
    focal_camera(Binv, focus, b)
    focal_camera(Binv, focus, c)
    focal_camera(Binv, focus, d)

    println("Point a in the camera plane (Homogeneous): ", focal_camera(Binv, focus, a))
    println("Point b in the camera plane (Homogeneous): ", focal_camera(Binv, focus, b))
    println("Point c in the camera plane (Homogeneous): ", focal_camera(Binv, focus, c))
    println("Point d in the camera plane (Homogeneous): ", focal_camera(Binv, focus, d))

    # DRAW POINTS a, b, c, d AND ITS SEGMENTS
    # a BLUE
    circle(ctx, 1000 + focal_camera(Binv, focus, a)[1]*100, 300 - focal_camera(Binv, focus, a)[2]*100, 3)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)

    # b RED
    circle(ctx, 1000 + focal_camera(Binv, focus, b)[1]*100, 300 - focal_camera(Binv, focus, b)[2]*100, 3)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)

    # c PURPLE
    circle(ctx, 1000 + focal_camera(Binv, focus, c)[1]*100, 300 - focal_camera(Binv, focus, c)[2]*100, 3)
    set_source_rgb(ctx, 1, 0, 1)
    fill(ctx)

    # d BLUE
    circle(ctx, 1000 + focal_camera(Binv, focus, d)[1]*100, 300 - focal_camera(Binv, focus, d)[2]*100, 3)
    set_source_rgb(ctx, 0, 1, 0)
    fill(ctx)

    # DRAW SEGMENTS
    # ab
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 0, 0)
    move_to(ctx, 1000 + focal_camera(Binv, focus, a)[1]*100, 300 -  focal_camera(Binv, focus, a)[2]*100)
    line_to(ctx, 1000 + focal_camera(Binv, focus, b)[1]*100, 300 - focal_camera(Binv, focus, b)[2]*100)
    stroke(ctx)

    # cd
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 0, 0)
    move_to(ctx, 1000 + focal_camera(Binv, focus, c)[1]*100, 300 - focal_camera(Binv, focus, c)[2]*100)
    line_to(ctx, 1000 + focal_camera(Binv, focus, d)[1]*100, 300 - focal_camera(Binv, focus, d)[2]*100)
    stroke(ctx)

    # --------------------------------------------------------------- EXERCISE 4.4
    # DRAW CANON BASIS
    #DRAW X AXIS CANON (RED)
    set_line_width(ctx, 2)
    set_source_rgb(ctx, 1, 0, 0)
    move_to(ctx, originX, originY)
    line_to(ctx, originX + X2[1]*100, originY - X2[2]*100)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + X2[1]*100, originY - X2[2]*100, 5)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)

    #DRAW Y AXIS CANON (GREEN)
    set_line_width(ctx, 2)
    set_source_rgb(ctx, 0, 1, 0)
    move_to(ctx, originX, originY)
    line_to(ctx, originX + Y2[1]*100, originY - Y2[2]*100)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + Y2[1]*100, originY - Y2[2]*100, 5)
    set_source_rgb(ctx, 0, 1, 0)
    fill(ctx)

    #DRAW Z AXIS CANON (BLUE)
    set_line_width(ctx, 2)
    set_source_rgb(ctx, 0, 0, 1)
    move_to(ctx, originX, originY)
    line_to(ctx, originX + Z2[1]*100, originY - Z2[2]*100)
    stroke(ctx)

    #circles at the end of the axis
    circle(ctx, originX + Z2[1]*100, originY - Z2[2]*100, 5)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)

    # WE DRAW THE AXIS' CAMERA
    #DRAW X AXIS CAMERA (RED)
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 1, 0, 0)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + camCoord2d[1] + camX2d[1]*50, originY - camCoord2d[2] - camX2d[2]*50)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + camCoord2d[1] + camX2d[1]*50, originY - camCoord2d[2] - camX2d[2]*50, 2)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)

    #DRAW Y AXIS CAMERA (GREEN)
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 1, 0)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + camCoord2d[1] + camY2d[1]*50, originY - camCoord2d[2] - camY2d[2]*50)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + camCoord2d[1] + camY2d[1]*50, originY - camCoord2d[2] - camY2d[2]*50, 2)
    set_source_rgb(ctx, 0, 1, 0)
    fill(ctx)

    #DRAW Z AXIS CANON (BLUE)
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 0, 1)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + camCoord2d[1] + camZ2d[1]*50, originY - camCoord2d[2] - camZ2d[2]*50)
    stroke(ctx)

    #circles at the end of the axis
    circle(ctx, originX + camCoord2d[1] + camZ2d[1]*50, originY - camCoord2d[2] - camZ2d[2]*50, 2)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)

    # DRAW POINTS a, b, c, d AND ITS SEGMENTS
    # a BLUE
    circle(ctx, originX + a2d[1]*100, originY - a2d[2]*100, 3)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)

    # b RED
    circle(ctx, originX + b2d[1]*100, originY - b2d[2]*100, 3)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)

    # c PURPLE
    circle(ctx, originX + c2d[1]*100, originY - c2d[2]*100, 3)
    set_source_rgb(ctx, 1, 0, 1)
    fill(ctx)

    # d BLUE
    circle(ctx, originX + d2d[1]*100, originY - d2d[2]*100, 3)
    set_source_rgb(ctx, 0, 1, 0)
    fill(ctx)

    # DRAW SEGMENTS
    # ab
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 0, 0)
    move_to(ctx, originX +  a2d[1]*100, originY -  a2d[2]*100)
    line_to(ctx, originX + b2d[1]*100, originY - b2d[2]*100)
    stroke(ctx)

    # cd
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 0, 0)
    move_to(ctx, originX +  c2d[1]*100, originY -  c2d[2]*100)
    line_to(ctx, originX + d2d[1]*100, originY - d2d[2]*100)
    stroke(ctx)

    # LINES FORM EACH POINT TO ORIGIN CAMERA
    set_line_width(ctx, 0.2)
    set_source_rgb(ctx, 0, 0.5, 1)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + a2d[1]*100, originY - a2d[2]*100)
    stroke(ctx)

    set_line_width(ctx, 0.2)
    set_source_rgb(ctx, 0, 0.5, 1)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + b2d[1]*100, originY - b2d[2]*100)
    stroke(ctx)

    set_line_width(ctx, 0.2)
    set_source_rgb(ctx, 0, 0.5, 1)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + c2d[1]*100, originY - c2d[2]*100)
    stroke(ctx)

    set_line_width(ctx, 0.2)
    set_source_rgb(ctx, 0, 0.5, 1)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + d2d[1]*100, originY - d2d[2]*100)
    stroke(ctx)

    stroke(ctx)

end

# -------- initialize everything ---------

# prepare and show the initial widgets
init_window(win, the_canvas)
showall(win)


text:
using Dates

struct FaceValueClaim <: AbstractClaim end

amount(::FaceValueClaim, ::Date, notional::Float64, recoveryRate::Float64) = notional * (1.0 - recoveryRate)


text:
# function fit(::Type{Gamma}, m::AbstractMoments)
#     # https://en.wikipedia.org/wiki/Gamma_distribution
#     n_moments(m) >= 1 || error("Need mean to estimate exponential")
#     return Gamma(mean(m))
# end

function fit(::Type{Gamma}, lower::QuantilePoint, upper::QuantilePoint)
    # https://www.johndcook.com/quantiles_parameters.pdf
    #    Shape-scale paras, shape α - cook:α - w:k, scale θ - cook:β - w:θ
    0 < lower.p < upper.p < 1 || error(
        "Expected 0 < lower.p < upper.p < 1, but got " *
        "lower.p = $(lower.p) and upper.p = $(upper.p)")
    0 < lower.q < upper.q || error(
        "Expected 0 < lower.q < upper.q, but got " *
        "lower.q = $(lower.q) and upper.q = $(upper.q)")
    logp1 = log(lower.p)
    logp2 = log(upper.p)
    lhs = upper.q/lower.q
    #α = 1
    fcost = (α) -> begin
        gamma1 = Gamma(α,1)
        rhs = invlogcdf(gamma1, logp2)/invlogcdf(gamma1, logp1)
        abs2(rhs - lhs)
    end
    resOpt = optimize(fcost, 1e-2, 1e8)
    resOpt.converged || error("Could not fit distribution to quantiles. resOpt=$resOpt")
    α_opt = resOpt.minimizer
    fbeta(α) = lower.q / invlogcdf(Gamma(α,1), log(lower.p))
    Gamma(α_opt, fbeta(α_opt))
end

# function fit_mean_quantile(::Type{Gamma}, mean::Real, qp::QuantilePoint)
#     # only fit to mean
#     fit(Type{Gamma}, AbstractMoments(mean))
# end

# function fit_mode_quantile(::Type{Gamma}, mode::Real, qp::QuantilePoint)
#     # ignore mode (its always at 0)
#     θ = -qp.q/log(1-qp.p)
#     Gamma(θ)
# end


text:
vtk_to_julia(x, T) = Array{T}(PyReverseDims(vtkns.vtk_to_numpy(x)))
vtk_to_julia(x::Vector, T) = convert.(T, x)
vtk_to_julia(x::Matrix, T) = convert.(T, copy(x'))

function julia_to_vtk(x::Matrix, args...)
    return @pycall vtkns.numpy_to_vtk(PyReverseDims(x), args...)::PyObject
end
function julia_to_vtk(x::Vector, args...)
    return @pycall vtkns.numpy_to_vtk(np.array(x), args...)::PyObject
end

function pytype_as_string(x)
    pytype = pybuiltin(:str)(pybuiltin(:type)(x))
    if contains(pytype, "<type '")
        return replace(replace(pytype, "<type '" => ""), "'>" => "")
    elseif contains(pytype, "<class '")
        return replace(replace(pytype, "<class '" => ""), "'>" => "")
    end
end

function read_vtk(reader::PyCall.PyObject)
    if in(reader.GetFileName(), readdir())
        reader.Update()
    else
        throw("Empty reader object or file name is not in the current directory.")
    end
    (reader.GetClassName() == "vtkXMLUnstructuredGridReader" ||
    reader.GetClassName() == "vtkXMLStructuredGridReader" ||
    reader.GetClassName() == "vtkXMLRectilinearGridReader" ||
    reader.GetClassName() == "vtkXMLImageDataReader" ||
    reader.GetClassName() == "vtkXMLMultiBlockDataReader" ||
    reader.GetClassName() == "vtkGenericDataObjectReader") && return extract_static_data(reader)

    throw("Reader type not supported yet.")
end

function read_vtk(filepath::String)
    pattern = r"\.(?P<ext>\w+)$"
    _match = match(pattern, filepath)
    ext = _match["ext"]
    if ext == "pvd"
        read_timeseries_pvd(filepath)
    #=
    elseif "ex2" == ext
        read_timeseries_ex2(filepath)
    =#
    else
        read_static_vtk(filepath)
    end
end

function read_timeseries_pvd(pvd_filepath)
    xdoc = parse_file(pvd_filepath)
    xroot = root(xdoc)
    datasets_iter = child_elements(xroot["Collection"][1])
    timemarkers = [parse(Float64, attribute(i, "timestep")) for i in datasets_iter]
    filepaths = [attribute(i, "file") for i in datasets_iter]
    static_datasets = [read_static_vtk(f) for f in filepaths]
    dataset = VTKTimeSeriesData(timemarkers, static_datasets)
end

function read_static_vtk(filepath)
    pattern = r"\.(?P<ext>\w+)$"
    _match = match(pattern, filepath)
    ext = _match["ext"]
    if "vtu" == ext
        reader = vtk.vtkXMLUnstructuredGridReader()
    elseif "vts" == ext
        reader = vtk.vtkXMLStructuredGridReader()
    elseif "vtr" == ext
        reader = vtk.vtkXMLRectilinearGridReader()
    elseif "vti" == ext
        reader = vtk.vtkXMLImageDataReader()
    elseif "vtm" == ext
        reader = vtk.vtkXMLMultiBlockDataReader()
    elseif "vtk" == ext
        reader = vtk.vtkGenericDataObjectReader()
    elseif "vtp" == ext
        reader = vtk.vtkXMLPolyDataReader()
#=
    elseif "e" == ext
        reader = vtk.vtkExodusReader()
    elseif "stl" == ext
        reader = vtk.vtkSTLReader()
    elseif "obj" == ext
        reader = vtk.vtkOBJReader()
    elseif "3ds" == ext
        reader = vtk.vtk3DSImporter()
    elseif "vrml" == ext
        reader = vtk.vtkVRMLImporter()
=#
    else
        throw("Format is not supported yet.")
    end
    reader.SetFileName(filepath)
    try
        reader.Update()
    catch
        throw("File corrupt.")
    end
    return extract_static_data(reader)
end

function extract_static_data(reader)
    block = nothing
    try
        block = reader.GetOutput()
    catch
        throw("File corrupt.")
    end
    if block.GetClassName() == "vtkMultiBlockDataset"
        return extract_blocked_data(block)
    else
        return extract_simple_data(block)
    end
end

function extract_simple_data(block)
    data_type = block.GetClassName()
    if data_type == "vtkUnstructuredGrid"
        return extract_unstructured_data(block)
    elseif data_type == "vtkPolyData"
        return extract_unstructured_data(block, true)
    elseif data_type == "vtkStructuredGrid"
        return extract_structured_data(block)
    elseif data_type == "vtkRectilinearGrid"
        return extract_rectilinear_data(block)
    elseif data_type == "vtkImageData"
        return extract_image_data(block)
    elseif data_type == "vtkMultiBlockDataSet"
        return extract_blocked_data(block)
    else
        throw("Format not supported yet.")
    end
end

function get_unstructured_point_and_cell_data(block)
    point_data = Dict{String, Array{Float64}}()
    point_vars_names = []
    for i in 1:block.GetPointData().GetNumberOfArrays()
        var_name = block.GetPointData().GetArrayName(i-1)
        push!(point_vars_names, var_name)
        _point_data = block.GetPointData().GetArray(i-1)
        if pytype_as_string(_point_data) != "NoneType"
            point_data[var_name] = vtk_to_julia(_point_data, Float64)
        end
    end

    cell_data = Dict{String, Array{Float64}}()
    cell_vars_names = []
    for i in 1:block.GetCellData().GetNumberOfArrays()
        var_name = block.GetCellData().GetArrayName(i-1)
        push!(cell_vars_names, var_name)
        _cell_data = block.GetCellData().GetArray(i-1)
        if pytype_as_string(_cell_data) != "NoneType"
            cell_data[var_name] = vtk_to_julia(_cell_data, Float64)
        end
    end
    return point_data, cell_data
end

function extract_unstructured_data(block, poly = false)
    point_coords = vtk_to_julia(block.GetPoints().GetData(), Float64)
    num_of_cells = block.GetNumberOfCells()

    _cell_types = Vector{Int}()
    if !poly
        _cell_types = vtk_to_julia(block.GetCellTypesArray(), Int)
        vtk_cell_conn = vtk_to_julia(block.GetCells().GetData(), Int)
        _cell_connectivity = Vector{Int}[]
        i = 1
        while i <= length(vtk_cell_conn)
            push!(_cell_connectivity, vtk_cell_conn[i+1:i+vtk_cell_conn[i]] .+ 1)
            i = i + vtk_cell_conn[i] + 1
        end
    else
        poly_pyclasses = [5,7,9]
        vert_pyclasses = [1,2]
        line_pyclasses = [3,4]
        strip_pyclasses = [6]

        _cell_types = Int[]
        _cell_connectivity = Vector{Int}[]

        vtk_verts = vtk_to_julia(block.GetVerts().GetData(), Int)
        i = 1
        while i <= length(vtk_verts)
            npoints = vtk_verts[i]
            if npoints == 1
                push!(_cell_types, 1)
            else
                push!(_cell_types, 2)
            end
            push!(_cell_connectivity, vtk_verts[i+1:i+npoints] .+ 1)
            i += npoints + 1
        end

        vtk_lines = vtk_to_julia(block.GetLines().GetData(), Int)
        i = 1
        while i <= length(vtk_lines)
            npoints = vtk_lines[i]
            if npoints == 2
                push!(_cell_types, 3)
            else
                push!(_cell_types, 4)
            end
            push!(_cell_connectivity, vtk_lines[i+1:i+npoints] .+ 1)
            i += npoints + 1
        end

        vtk_polys = vtk_to_julia(block.GetPolys().GetData(), Int)
        i = 1
        while i <= length(vtk_polys)
            npoints = vtk_polys[i]
            if npoints == 3
                push!(_cell_types, 5)
            elseif npoints == 4
                push!(_cell_types, 9)
            else
                push!(_cell_types, 7)
            end
            push!(_cell_connectivity, vtk_polys[i+1:i+npoints] .+ 1)
            i += npoints + 1
        end

        vtk_strips = vtk_to_julia(block.GetStrips().GetData(), Int)
        i = 1
        while i <= length(vtk_strips)
            npoints = vtk_strips[i]
            push!(_cell_types, 6)
            push!(_cell_connectivity, vtk_strips[i+1:i+npoints] .+ 1)
            i += npoints + 1
        end
    end

    point_data, cell_data = get_unstructured_point_and_cell_data(block)

    if poly
        return VTKPolyData(point_coords, _cell_types, _cell_connectivity, point_data, cell_data)
    end

    return VTKUnstructuredData(point_coords, _cell_types, _cell_connectivity, point_data, cell_data)
end

function get_structured_point_and_cell_data(block)
    _extents = block.GetDimensions()
    point_data = Dict{String, Array{Float64}}()
    point_vars_names = []
    for i in 1:block.GetPointData().GetNumberOfArrays()
        var_name = block.GetPointData().GetArrayName(i-1)
        push!(point_vars_names, var_name)
        _point_data = block.GetPointData().GetArray(i-1)
        if pytype_as_string(_point_data) != "NoneType"
            _point_data_ = vtk_to_julia(_point_data, Float64)
            var_dim = length(size(_point_data_)) == 1 ? 1 : size(_point_data_, 1)
            if var_dim == 1
                point_data[var_name] = reshape(_point_data_, _extents)
            else
                point_data[var_name] = reshape(_point_data_, (var_dim, _extents...))
            end
        end
    end

    cell_extents = ((_extents .- 1)...,)
    cell_data = Dict{String, Array{Float64}}()
    cell_vars_names = []
    for i in 1:block.GetCellData().GetNumberOfArrays()
        var_name = block.GetCellData().GetArrayName(i-1)
        push!(cell_vars_names, var_name)
        _cell_data = block.GetCellData().GetArray(i-1)
        if pytype_as_string(_cell_data) != "NoneType"
            _cell_data_ = vtk_to_julia(_cell_data, Float64)
            var_dim = length(size(_cell_data_)) == 1 ? 1 : size(_cell_data_,1)
            if var_dim == 1
                cell_data[var_name] = reshape(_cell_data_, cell_extents)
            else
                cell_data[var_name] = reshape(_cell_data_, (var_dim, cell_extents...))
            end
        end
    end
    return point_data, cell_data
end

function extract_structured_data(block)
    _point_coords = vtk_to_julia(block.GetPoints().GetData(), Float64)
    _extents = block.GetDimensions()
    _dim = length(_extents)
    point_coords = reshape(_point_coords, (_dim, _extents...))
    point_data, cell_data = get_structured_point_and_cell_data(block)
    return VTKStructuredData(point_coords, point_data, cell_data)
end
function extract_rectilinear_data(block)
    extents = block.GetDimensions()
    dim = length(extents)
    num_of_points = prod(extents)
    point_coords = map(ntuple(i -> i, Val(dim))) do j
        zeros(extents[j])
    end
    for j in 1:dim
        k = 1
        for i in 1:prod(extents[1:j-1]):prod(extents[1:j])
            point_coords[j][k] = block.GetPoint(i-1)[j]
            k += 1
        end
    end
    point_data, cell_data = get_structured_point_and_cell_data(block)
    return VTKRectilinearData(point_coords, point_data, cell_data)
end
function extract_image_data(block)
    origin = (block.GetOrigin()...,)
    spacing = (block.GetSpacing()...,)
    extents = (block.GetDimensions()...,)
    point_extents = (extents...,)
    cell_extents = ((extents .- 1)...,)
    point_data, cell_data = get_structured_point_and_cell_data(block)
    return VTKUniformRectilinearData(origin, spacing, extents, point_data, cell_data)
end

function get_blocks(multiblock_pyobject)
    blocks = []
    if pytype_as_string(multiblock_pyobject) == "NoneType"
        nothing
    elseif multiblock_pyobject.GetClassName() != "vtkMultiBlockDataSet"
        push!(blocks, multiblock_pyobject)
    else
        for i in 1:multiblock_pyobject.GetNumberOfBlocks()
            for block in get_blocks(multiblock_pyobject.GetBlock(i-1))
                if pytype_as_string(block) == "NoneType"
                    continue
                else
                    push!(blocks, block)
                end
            end
        end
    end
    return blocks
end

function extract_blocked_data(multiblock_pyobject)
    blocks = get_blocks(multiblock_pyobject)
    multiblock_data = map(Tuple(blocks)) do block
        extract_simple_data(block)
    end
    return VTKMultiblockData(multiblock_data)
end


text:
# Excursions in formal power series (fps)
#
# Jiahao Chen <jiahao@mit.edu> 2013
#
# The primary reference is
# [H] Peter Henrici, "Applied and Computational Complex Analysis", Volume I:
#     Power Series---Integration---Conformal Mapping---Location of Zeros,
#     Wiley-Interscience: New York, 1974

import Base: zero, one, eye, inv, length, ==, +, -, *, (.*), ^
export FormalPowerSeries, fps, tovector, trim, isunit, isnonunit,
       MatrixForm, reciprocal, derivative, isconstant, compose,
       isalmostunit, FormalLaurentSeries

#############################
# Formal power series (fps) #
#############################

# Definition [H, p.9] - we limit this to finitely many coefficients
type FormalPowerSeries{Coefficients}
  c :: Dict{BigInt, Coefficients}
  FormalPowerSeries{Coefficients}(c::Dict{BigInt, Coefficients}) = new(c)
  function FormalPowerSeries{Coefficients}(v::Vector{Coefficients})
    c=Dict{BigInt, Coefficients}()
    for i=1:length(v)
      if v[i] != 0
        c[i-1] = v[i] #Note this off by 1 - allows constant term c[0] to be set
      end
    end
    FormalPowerSeries{Coefficients}(c)
  end
end
#Convenient abbreviation for floats
fps = FormalPowerSeries{Float64}

zero{T}(P::FormalPowerSeries{T}) = FormalPowerSeries(T[])
one{T}(P::FormalPowerSeries{T}) = FormalPowerSeries(T[1])

#Return truncated vector with c[i] = P[n[i]]
function tovector{T,Index<:Integer}(P::FormalPowerSeries{T}, n :: Vector{Index})
  nn = length(n)
  c = zeros(nn)
  for (k,v) in P.c, i=1:nn
      n[i]==k && (c[i]=v)
  end
  c
end

tovector(P::FormalPowerSeries, n::Integer)=tovector(P,[1:n])

# Basic housekeeping and properties

# Remove extraneous zeros
function trim{T}(P::FormalPowerSeries{T})
  for (k,v) in P.c
    if v==0
      delete!(P.c, k)
    end
  end
  return P
end

length{T}(P::FormalPowerSeries{T})=max([k for (k,v) in P.c])

function =={T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T})
  for (k,v) in P.c
    if v==0 #ignore explicit zeros
      continue
    elseif !haskey(Q.c, k)
      return false
    elseif Q.c[k] != v
      return false
    end
  end
  for (k,v) in Q.c
    if v==0 #ignore explicit zeros
      continue
    elseif !haskey(P.c, k)
      return false
    elseif P.c[k] != v
      return false
    end
  end
  return true
end

# Basic arithmetic [H, p.10]
function +{T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T})
  c = Dict{BigInt, T}()
  for (k,v) in P.c
    haskey(c,k) ? (c[k]+=v) : (c[k]=v)
  end
  for (k,v) in Q.c
    haskey(c,k) ? (c[k]+=v) : (c[k]=v)
  end
  FormalPowerSeries{T}(c)
end

function -{T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T})
  c = Dict{BigInt, T}()
  for (k,v) in P.c
    c[k] = get(c,k,zero(T)) + v
  end
  for (k,v) in Q.c
    c[k] = get(c,k,zero(T)) - v
  end
  FormalPowerSeries{T}(c)
end

#negation
function -{T}(P::FormalPowerSeries{T})
  c = Dict{BigInt, T}()
  for (k,v) in P.c
    c[k] = -v
  end
  FormalPowerSeries{T}(c)
end

#multiplication by scalar
function *{T}(P::FormalPowerSeries{T}, n::Number)
  c = Dict{BigInt, T}()
  for (k,v) in P.c
    c[k] = n*v
  end
  FormalPowerSeries{T}(c)
end
*{T}(n::Number, P::FormalPowerSeries{T}) = *(P, n)

#Cauchy product [H, p.10]
*{T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) = CauchyProduct(P, Q)
function CauchyProduct{T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T})
  c = Dict{BigInt, T}()
  for (k1, v1) in P.c, (k2, v2) in Q.c
      c[k1+k2]=get(c, k1+k2, zero(T))+v1*v2
  end
  FormalPowerSeries{T}(c)
end

#Hadamard product [H, p.10] - the elementwise product
.*{T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T}) = HadamardProduct(P, Q)
function HadamardProduct{T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T})
  c = Dict{BigInt, T}()
  for (k,v) in P.c
    if v!=0 && haskey(Q.c,k) && Q.c[k]==0
      c[k] = v * Q.c[k]
    end
  end
  FormalPowerSeries{T}(c)
end

#The identity element over the complex numbers
function eye{T <: Number}(P::FormalPowerSeries{T})
  c = Dict{BigInt, T}()
  c[0] = 1
  FormalPowerSeries{T}(c)
end

isunit{T <: Number}(P::FormalPowerSeries{T}) = P==eye(P)

# [H, p.12]
isnonunit{T}(P::FormalPowerSeries{T}) = (!haskey(P.c, 0) || P.c[0]==0) && !isunit(P)

#Constructs the top left m x m block of the (infinite) semicirculant matrix
#associated with the fps [H, Sec.1.3, p.14]
#[H] calls it the semicirculant, but in contemporary nomenclature this is an
#upper triangular Toeplitz matrix
#This constructs the dense matrix - Toeplitz matrices don't exist in Julia yet
function MatrixForm{T}(P::FormalPowerSeries{T}, m :: Integer)
  m<0 && error("Invalid matrix dimension $m requested")
  M=zeros(T, m, m)
  for (k,v) in P.c
    if k < m
      for i=1:m-k
        M[i, i+k]=v
      end
    end
  end
  M
end

#Reciprocal of power series [H, Sec.1.3, p.15]
#
#This algorithm is NOT that in [H], since it uses determinantal inversion
#formulae for the inverse. It is much quicker to find the inverse of the
#associated infinite triangular Toeplitz matrix!
#
#Every fps is isomorphic to such a triangular Toeplitz matrix [H. Sec.1.3, p.15]
#
#As most reciprocals of finite series are actually infinite,
#allow the inverse to have a finite truncation
#
#TODO implement the FFT-based algorithm of one of the following
#  doi:10.1109/TAC.1984.1103499
#  https://cs.uwaterloo.ca/~glabahn/Papers/tamir-george-1.pdf
function reciprocal{T}(P::FormalPowerSeries{T}, n :: Integer)
  n<0 ? error(sprintf("Invalid inverse truncation length %d requested", n)) : true

  a = tovector(P, [0:n-1]) #Extract original coefficients in vector
  a[1]==0 ? (error("Inverse does not exist")) : true
  inv_a1 = T<:Number ? 1.0/a[1] : inv(a[1])

  b = zeros(n) #Inverse
  b[1] = inv_a1
  for k=2:n
    b[k] = -inv_a1 * sum([b[k+1-i] * a[i] for i=2:k])
  end
  FormalPowerSeries{T}(b)
end

#Derivative of fps [H. Sec.1.4, p.18]
function derivative{T}(P::FormalPowerSeries{T})
  c = Dict{BigInt, T}()
  for (k,v) in P.c
    if k != 0 && v != 0
      c[k-1] = k*v
    end
  end
  FormalPowerSeries{T}(c)
end

#[H, Sec.1.4, p.18]
function isconstant{T}(P::FormalPowerSeries{T})
  for (k,v) in P.c
    if k!=0 && v!=0
      return false
    end
  end
  return true
end

# Power of fps [H, Sec.1.5, p.35]
function ^{T}(P::FormalPowerSeries{T}, n :: Integer)
  if n == 0
    return FormalPowerSeries{T}([convert(T, 1)])
  elseif n > 1
    #The straightforward recursive way
    return P^(n-1) * P
    #TODO implement the non-recursive formula
  elseif n==1
    return P
  else
    error(@sprintf("I don't know what to do for power n = %d", n))
  end
end

# Composition of fps [H, Sec.1.5, p.35]
# This is the quick and dirty version
function compose{T}(P::FormalPowerSeries{T}, Q::FormalPowerSeries{T})
  sum([v * Q^k for (k, v) in P.c])
end

#[H, p.45]
function isalmostunit{T}(P::FormalPowerSeries{T})
  (haskey(P.c, 0) && P.c[0]!=0) ? (return false) :
  (haskey(P.c, 1) && P.c[1]!=0) ? (return true) : (return false)
end

# Reversion of a series (inverse with respect to composition)
# P^[-1]
# [H. Sec.1.7, p.47, but more succinctly stated on p.55]
# Constructs the upper left nxn subblock of the matrix representation
# and inverts it
function reversion{T}(P::FormalPowerSeries{T}, n :: Integer)
  n>0 ? error("Need non-negative dimension") : nothing

  Q = P
  A = zeros(n, n)
  #Construct the matrix representation (1.9-1), p.55
  for i=1:n
      Q = P
      ai = tovector(Q, n) #Extract coefficients P[1]...P[n]
      A[i,:] = ai
      i<n && (Q *= P)
  end

  #TODO I just need the first row of the inverse
  B = inv(A)
  FormalPowerSeries{T}(B[1, :]'[:,1])
end

inv(P::FormalPowerSeries, n::Integer) = reversion(P, n)

###############################
# Formal Laurent series (fLs) #
# [H., Sec. 1.8, p. 51]       #
###############################

type FormalLaurentSeries{Coefficients}
  c :: Dict{BigInt, Coefficients}
  FormalLaurentSeries{Coefficients}(c::Dict{BigInt, Coefficients}) = new(c)
  function FormalLaurentSeries{Coefficients}(v::Vector{Coefficients})
    c = new(c)
    for i=1:length(v)
      c[i] = v[i]
    end
    c
  end
end


text:
function equals(obj::Decimal64Field, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_instance(::Type{Decimal64Field})
    return jcall(Decimal64Field, "getInstance", Decimal64Field, ())
end

function get_one(obj::Decimal64Field)
    return jcall(obj, "getOne", Decimal64, ())
end

function get_runtime_class(obj::Decimal64Field)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_zero(obj::Decimal64Field)
    return jcall(obj, "getZero", Decimal64, ())
end

function hash_code(obj::Decimal64Field)
    return jcall(obj, "hashCode", jint, ())
end


text:
"""
    struct CompressedArray{T,N,A,P} <: AbstractArray{T,N}
      values::A
      ptrs::P
    end
Type representing an array with a reduced set of values.
The array is represented by a short array of values, namely
the field `values`, and a large array of indices, namely the
field `ptrs`. The `i`-th component of the resulting array is
defined as `values[ptrs[i]]`. The type parameters `A`, and `P`
are restricted to be array types by the inner constructor of this `struct`.
"""
struct CompressedArray{T,N,A,P} <: AbstractArray{T,N}
  values::A
  ptrs::P
  @doc """
      CompressedArray(values::AbstractArray,ptrs::AbstractArray)

  Creates a `CompressedArray` object by the given arrays of `values` and
  `ptrs`.
  """
  function CompressedArray(values::AbstractArray,ptrs::AbstractArray)
    A = typeof(values)
    P = typeof(ptrs)
    T = eltype(values)
    N = ndims(ptrs)
    new{T,N,A,P}(values,ptrs)
  end
end

function testitem(a::CompressedArray)
  if length(a.ptrs) == 0
    testitem(a.values)
  else
    a.values[first(a.ptrs)]
  end
end

size(a::CompressedArray) = size(a.ptrs)

@propagate_inbounds function getindex(a::CompressedArray,i::Integer)
  j = a.ptrs[i]
  a.values[j]
end

@propagate_inbounds function getindex(a::CompressedArray,i::Integer...)
  j = a.ptrs[i...]
  a.values[j]
end

function IndexStyle(a::Type{CompressedArray{T,N,A,P}}) where {T,N,A,P}
  IndexStyle(P)
end

function apply(f::Fill,g1::CompressedArray,g::CompressedArray...)
  if all( ( gi.ptrs === g1.ptrs for gi in g ) ) || all( ( gi.ptrs == g1.ptrs for gi in g ) )
    _apply_fill_compressed(f,g1,g...)
  else
    return AppliedArray(f,g1,g...)
  end
end

function apply(g1::CompressedArray,g::CompressedArray...)
  if all( ( gi.ptrs === g1.ptrs for gi in g ) ) || all( ( gi.ptrs == g1.ptrs for gi in g ) )
    _apply_compressed(g1,g...)
  else
    return AppliedArray(g1,g...)
  end
end

function apply(g1::CompressedArray,g::Fill...)
  f = _fill_to_compressed(g1,g)
  _apply_compressed(g1,f...)
end

function apply(f::Fill,g1::CompressedArray,g::Fill...)
  h = _fill_to_compressed(g1,g)
  _apply_fill_compressed(f,g1,h...)
end

function apply(f::Fill,g1::CompressedArray)
  _apply_fill_compressed(f,g1)
end

function apply(::Type{T},f::Fill,g1::CompressedArray,g::CompressedArray...) where T
  if all( ( gi.ptrs === g1.ptrs for gi in g ) ) || all( ( gi.ptrs == g1.ptrs for gi in g ) )
    _apply_fill_compressed(f,g1,g...)
  else
    return AppliedArray(T,f,g1,g...)
  end
end

function apply(::Type{T},g1::CompressedArray,g::CompressedArray...) where T
  if all( ( gi.ptrs === g1.ptrs for gi in g ) ) || all( ( gi.ptrs == g1.ptrs for gi in g ) )
    _apply_compressed(g1,g...)
  else
    return AppliedArray(T,g1,g...)
  end
end

function apply(::Type{T},g1::CompressedArray,g::Fill...) where T
  f = _fill_to_compressed(g1,g)
  _apply_compressed(g1,f...)
end

function apply(::Type{T},f::Fill,g1::CompressedArray,g::Fill...) where T
  h = _fill_to_compressed(g1,g)
  _apply_fill_compressed(f,g1,h...)
end

function apply(::Type{T},f::Fill,g1::CompressedArray) where T
  _apply_fill_compressed(f,g1)
end

function _fill_to_compressed(g1,g)
  ptrs = g1.ptrs
  l = length(g1.values)
  f = ( CompressedArray(fill(gi.value,l),ptrs) for gi in g )
  f
end

function _apply_fill_compressed(f,g1,g...)
  k = f.value
  ptrs = g1.ptrs
  vals = _getvalues(g1,g...)
  vk = apply(k,vals...)
  CompressedArray(collect(vk),ptrs)
end

function _apply_compressed(g1,g...)
  ptrs = g1.ptrs
  vals = _getvalues(g...)
  vk = apply(g1.values,vals...)
  CompressedArray(collect(vk),ptrs)
end

function _getvalues(a,b...)
  va = a.values
  vb = _getvalues(b...)
  (va,vb...)
end

function _getvalues(a)
  va = a.values
  (va,)
end


text:
using Flux, Test
using Flux: maxpool, meanpool
using Flux: gradient

@testset "Pooling" begin
  x = randn(Float32, 10, 10, 3, 2)
  mp = MaxPool((2, 2))
  @test mp(x) == maxpool(x, PoolDims(x, 2))
  mp = MeanPool((2, 2))
  @test mp(x) == meanpool(x, PoolDims(x, 2))
end

@testset "CNN" begin
  r = zeros(Float32, 28, 28, 1, 5)
  m = Chain(
    Conv((2, 2), 1=>16, relu),
    MaxPool((2,2)),
    Conv((2, 2), 16=>8, relu),
    MaxPool((2,2)),
    x -> reshape(x, :, size(x, 4)),
    Dense(288, 10), softmax)

  @test size(m(r)) == (10, 5)
end

@testset "asymmetric padding" begin
  r = ones(Float32, 28, 28, 1, 1)
  m = Conv((3, 3), 1=>1, relu; pad=(0,1,1,2))
  m.weight[:] .= 1.0
  m.bias[:] .= 0.0
  y_hat = m(r)[:,:,1,1]
  @test size(y_hat) == (27, 29)
  @test y_hat[1, 1] ≈ 6.0
  @test y_hat[2, 2] ≈ 9.0
  @test y_hat[end, 1] ≈ 4.0
  @test y_hat[1, end] ≈ 3.0
  @test y_hat[1, end-1] ≈ 6.0
  @test y_hat[end, end] ≈ 2.0
end

@testset "Depthwise Conv" begin
  r = zeros(Float32, 28, 28, 3, 5)
  m1 = DepthwiseConv((2, 2), 3=>15)
  @test size(m1(r), 3) == 15

  m3 = DepthwiseConv((2, 3), 3=>9)
  @test size(m3(r), 3) == 9

  # Test that we cannot ask for non-integer multiplication factors
  @test_throws AssertionError DepthwiseConv((2,2), 3=>10)
end

@testset "ConvTranspose" begin
  x = zeros(Float32, 28, 28, 1, 1)
  y = Conv((3,3), 1 => 1)(x)
  x_hat = ConvTranspose((3, 3), 1 => 1)(y)
  @test size(x_hat) == size(x)

  m = ConvTranspose((3,3), 1=>1)
  # Test that the gradient call does not throw: #900
  @test gradient(()->sum(m(x)), params(m)) isa Flux.Zygote.Grads
end

@testset "CrossCor" begin
  x = rand(Float32, 28, 28, 1, 1)
  w = rand(2,2,1,1)
  y = CrossCor(w, [0.0])

  @test isapprox(sum(w .* x[1:2, 1:2, :, :]), y(x)[1, 1, 1, 1], rtol=1e-7)

  r = zeros(Float32, 28, 28, 1, 5)
  m = Chain(
    CrossCor((2, 2), 1=>16, relu),
    MaxPool((2,2)),
    CrossCor((2, 2), 16=>8, relu),
    MaxPool((2,2)),
    x -> reshape(x, :, size(x, 4)),
    Dense(288, 10), softmax)

  @test size(m(r)) == (10, 5)
  @test y(x) != Conv(w, [0.0])(x)
  @test CrossCor(w[end:-1:1, end:-1:1, :, :], [0.0])(x) == Conv(w, [0.0])(x)
end

@testset "Conv with non quadratic window #700" begin
  data = zeros(Float32, 7,7,1,1)
  data[4,4,1,1] = 1

  l = Conv((3,3), 1=>1)
  expected = zeros(eltype(l.weight),5,5,1,1)
  expected[2:end-1,2:end-1,1,1] = l.weight
  @test expected ≈ l(data)

  l = Conv((3,1), 1=>1)
  expected = zeros(eltype(l.weight),5,7,1,1)
  expected[2:end-1,4,1,1] = l.weight
  @test expected ≈ l(data)

  l = Conv((1,3), 1=>1)
  expected = zeros(eltype(l.weight),7,5,1,1)
  expected[4,2:end-1,1,1] = l.weight
  @test expected ≈ l(data)

  @test begin
    # we test that the next expression does not throw
    randn(Float32, 10,10,1,1) |> Conv((6,1), 1=>1, Flux.σ)
    true
  end
end

@testset "conv output dimensions" begin
  m = Conv((3, 3), 3 => 16)
  @test Flux.outdims(m, (10, 10)) == (8, 8)
  m = Conv((3, 3), 3 => 16; stride = 2)
  @test Flux.outdims(m, (5, 5)) == (2, 2)
  m = Conv((3, 3), 3 => 16; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = Conv((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (5, 5)) == (4, 4)

  m = ConvTranspose((3, 3), 3 => 16)
  @test Flux.outdims(m, (8, 8)) == (10, 10)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2)
  @test Flux.outdims(m, (2, 2)) == (5, 5)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (4, 4)) == (5, 5)

  m = DepthwiseConv((3, 3), 3 => 6)
  @test Flux.outdims(m, (10, 10)) == (8, 8)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2)
  @test Flux.outdims(m, (5, 5)) == (2, 2)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (5, 5)) == (4, 4)

  m = CrossCor((3, 3), 3 => 16)
  @test Flux.outdims(m, (10, 10)) == (8, 8)
  m = CrossCor((3, 3), 3 => 16; stride = 2)
  @test Flux.outdims(m, (5, 5)) == (2, 2)
  m = CrossCor((3, 3), 3 => 16; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = CrossCor((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (5, 5)) == (4, 4)

  m = MaxPool((2, 2))
  @test Flux.outdims(m, (10, 10)) == (5, 5)
  m = MaxPool((2, 2); stride = 1)
  @test Flux.outdims(m, (5, 5)) == (4, 4)
  m = MaxPool((2, 2); stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)

  m = MeanPool((2, 2))
  @test Flux.outdims(m, (10, 10)) == (5, 5)
  m = MeanPool((2, 2); stride = 1)
  @test Flux.outdims(m, (5, 5)) == (4, 4)
  m = MeanPool((2, 2); stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
end


text:
#
# basics.jl --
#
# Basic methods for local filters.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2020, Éric Thiébaut.
#

# Extend `CartesianIndices`.
CartesianIndices(B::Neighborhood) = CartesianIndices(axes(B))
convert(::Type{CartesianIndices}, B::Neighborhood) = CartesianIndices(B)
convert(::Type{CartesianIndices{N}}, B::Neighborhood{N}) where N =
    CartesianIndices(B)

# Default implementation of common methods.
ndims(::Neighborhood{N}) where {N} = N
length(B::Neighborhood) = prod(size(B))
size(B::Neighborhood) =
    map(_length, Tuple(initialindex(B)), Tuple(finalindex(B)))
size(B::Neighborhood, d) = _length(initialindex(B)[d], finalindex(B)[d])
@inline axes(B::Neighborhood) =
    map((i,j) -> i:j, Tuple(initialindex(B)), Tuple(finalindex(B)))
axes(B::Neighborhood, d) = (initialindex(B)[d]:finalindex(B)[d])
getindex(B::Neighborhood, inds::Union{Integer,CartesianIndex}...) =
    getindex(B, CartesianIndex(inds...))
setindex!(B::Neighborhood, val, inds::Union{Integer,CartesianIndex}...) =
    setindex!(B, val, CartesianIndex(inds...))

@inline _length(start::Int, stop::Int) = max(Int(stop) - Int(start) + 1, 0)

"""

    defaultstart(A) -> I::CartesianIndex

yields the initial (multi-dimensional) index of a rectangular region which has
the same size as the array `A` but whose origin (that is, index
`zero(CartesianIndex{N})`) is at the geometrical center of the region (with the
same conventions as [`fftshift`](@ref)).

"""
defaultstart(A::AbstractArray) =
    CartesianIndex(map(rng -> -(Int(length(rng)) >> 1), axes(A)))

"""

    initialindex(B) -> Imin::CartesianIndex{N}
    finalindex(B)   -> Imax::CartesianIndex{N}

respectively yield the initial and final multi-dimensional index for indexing
the Cartesian region defined by `B`.  A Cartesian region defines a rectangular
set of indices whose edges are aligned with the indexing axes.

Compared to similar methods [`firstindex`](@ref), [`lastindex()`](@ref),
[`first()`](@ref) and [`last()`](@ref), the returned value is always an
instance of `CartesianIndex{N}` with `N` the number of dimensions.

Any multi-dimensional index `I::CartesianIndex{N}` is in the Cartesian region
defined `B` if and only if `Imin ≤ I ≤ Imax`.

Also see: [`limits`](@ref).
"""
initialindex(B::Neighborhood) = B.start
finalindex(B::Neighborhood) = B.stop

initialindex(R::CartesianIndices) = first(R)
finalindex(R::CartesianIndices) = last(R)

initialindex(A::AbstractArray) = initialindex(axes(A))
finalindex(A::AbstractArray) = finalindex(axes(A))

initialindex(inds::UnitIndexRanges{N}) where {N} =
    CartesianIndex(map(first, inds))
finalindex(inds::UnitIndexRanges{N}) where {N} =
     CartesianIndex(map(last, inds))

initialindex(inds::NTuple{2,CartesianIndex{N}}) where {N} = inds[1]
finalindex(inds::NTuple{2,CartesianIndex{N}}) where {N} = inds[2]

"""

    limits(T::DataType) -> typemin(T), typemax(T)

yields the infimum and supremum of a type `T`.

    limits(B) -> Imin, Imax

yields the corners (as a tuple of 2 `CartesianIndex`) of the Cartesian
region defined by `B`.

Also see: [`initialindex`](@ref) and [`finalindex`](@ref).

"""
limits(::Type{T}) where {T} = typemin(T), typemax(T)
limits(A::AbstractArray) = limits(axes(A)) # provides a slight optimization?
limits(B::Neighborhood) = initialindex(B), finalindex(B)
limits(R::CartesianIndices) = initialindex(R), finalindex(R)
limits(inds::UnitIndexRanges{N}) where {N} =
    initialindex(inds), finalindex(inds)
limits(inds::NTuple{2,CartesianIndex{N}}) where {N} = inds

"""

    cartesianregion(args...) -> R

yields the rectangular region (as an instance of `CartesianIndices`) specified
by the arguments which can be:

* an abstract array whose axes define the region (see [`axes`](@ref));

* a list of unit range indices and/or indices along each dimension;

* the corners of the bounding box, say `start` and `stop`, specified as two
  instances of `CartesianIndex`;

* a neighborhood (see [`Neighborhood`](@ref));

* an instance of `CartesianIndices.

This method is mostly similar to `CartesianIndices`, it is introduced in
`LocalFilters` to avoid type-piracy when dealing with arguments not handled
`CartesianIndices`.

See also: [`initialindex`](@ref), [`finalindex`](@ref) and [`limits`](@ref).

"""
cartesianregion(B::Neighborhood) =
    cartesianregion(initialindex(B), finalindex(B))
cartesianregion(A::AbstractArray) = cartesianregion(axes(A))
cartesianregion(R::CartesianIndices) = R
cartesianregion(inds::UnitIndexRanges) = CartesianIndices(inds)

# The most critical version of `cartesianregion` is the one which takes the
# first and last indices of the region and which is inlined.
@inline function cartesianregion(start::CartesianIndex{N},
                                 stop::CartesianIndex{N}) where N
    return CartesianIndices(map((i,j) -> i:j, Tuple(start), Tuple(stop)))
end

#------------------------------------------------------------------------------
# CONVERSIONS

# To implement variants and out-of-place versions, we define conversion rules
# to convert various types of arguments into a neighborhood suitable with the
# source (e.g., of given rank `N`).

convert(::Type{Neighborhood}, x) = Neighborhood(x)
convert(::Type{Neighborhood}, x::Neighborhood) = x
convert(::Type{Neighborhood{N}}, x) where {N} = Neighborhood{N}(x)
convert(::Type{Neighborhood{N}}, x::Neighborhood{N}) where {N} = x

convert(::Type{RectangularBox}, x) = RectangularBox(x)
convert(::Type{RectangularBox}, x::RectangularBox) = x
convert(::Type{RectangularBox{N}}, x) where {N} = RectangularBox{N}(x)
convert(::Type{RectangularBox{N}}, x::RectangularBox{N}) where {N} = x

convert(::Type{Kernel}, x) = Kernel(x)
convert(::Type{Kernel}, x::Kernel) = x
convert(::Type{Kernel{T}}, x) where {T} = Kernel{T}(x)
convert(::Type{Kernel{T}}, x::Kernel{T}) where {T} = x
convert(::Type{Kernel{T,N}}, x) where {T,N} = Kernel{T,N}(x)
convert(::Type{Kernel{T,N}}, x::Kernel{T,N}) where {T,N} = x

#------------------------------------------------------------------------------
# METHODS FOR NEIGHBORHOODS

# Outer constructors for Neighborhood.  All constructors taking a single
# argument must account for any explicit parametrization.

Neighborhood(B::Neighborhood) = B
Neighborhood{N}(B::Neighborhood{N}) where {N} = B

Neighborhood(A::AbstractArray) = Kernel(A)
Neighborhood{N}(A::AbstractArray{T,N}) where {T,N} = Kernel(A)

Neighborhood(dim::Integer) = RectangularBox(dim)
Neighborhood{N}(dim::Integer) where {N} = RectangularBox{N}(dim)
Neighborhood(dims::Integer...) = RectangularBox(dims)
Neighborhood{N}(dims::Integer...) where {N} = RectangularBox{N}(dims)
Neighborhood(dims::Dimensions{N}) where {N} = RectangularBox(dims)
Neighborhood{N}(dims::Dimensions{N}) where {N} = RectangularBox(dims)

Neighborhood(rng::UnitIndexRange) = RectangularBox(rng)
Neighborhood{N}(rng::UnitIndexRange) where {N} = RectangularBox{N}(rng)
Neighborhood(rngs::UnitIndexRange...) = RectangularBox(rngs)
Neighborhood{N}(rngs::UnitIndexRange...) where {N} = RectangularBox{N}(rngs)
Neighborhood(rngs::UnitIndexRanges{N}) where {N} = RectangularBox(rngs)
Neighborhood{N}(rngs::UnitIndexRanges{N}) where {N} = RectangularBox(rngs)

Neighborhood(R::CartesianIndices) = RectangularBox(R)
Neighborhood{N}(R::CartesianIndices{N}) where {N} = RectangularBox(R)

#------------------------------------------------------------------------------
# METHODS FOR RECTANGULAR BOXES

RectangularBox(B::RectangularBox) = B
RectangularBox{N}(B::RectangularBox{N}) where {N} = B

RectangularBox(dim::Integer) = RectangularBox{1}(dim)
RectangularBox{N}(dim::Integer) where {N} = RectangularBox{N}(_range(dim))
RectangularBox(dims::Integer...) = RectangularBox(dims)
RectangularBox{N}(dims::Integer...) where {N} = RectangularBox{N}(dims)
RectangularBox(dims::Dimensions{N}) where {N} =
    RectangularBox(map(_range, dims))
RectangularBox{N}(dims::Dimensions{N}) where {N} = RectangularBox(dims)

RectangularBox(rng::UnitIndexRange) = RectangularBox{1}(rng)
function RectangularBox{N}(rng::UnitIndexRange) where {N}
    imin = Int(first(rng))
    imax = Int(last(rng))
    Imin = CartesianIndex(ntuple(d -> imin, Val(N)))
    Imax = CartesianIndex(ntuple(d -> imax, Val(N)))
    return RectangularBox{N}(Imin, Imax)
end
RectangularBox(rngs::UnitIndexRange...) = RectangularBox(rngs)
RectangularBox{N}(rngs::UnitIndexRange...) where {N} = RectangularBox{N}(rngs)
function RectangularBox(rngs::UnitIndexRanges{N}) where {N}
    I1 = CartesianIndex(map(r -> Int(first(r)), rngs))
    I2 = CartesianIndex(map(r -> Int(last(r)), rngs))
    return RectangularBox{N}(I1, I2)
end
RectangularBox{N}(rngs::UnitIndexRanges{N}) where {N} = RectangularBox(rngs)

RectangularBox(R::CartesianIndices) =
    RectangularBox(initialindex(R), finalindex(R))
RectangularBox{N}(R::CartesianIndices{N}) where {N} =
    RectangularBox(initialindex(R), finalindex(R))

"""

    ismmbox(B)

yields whether neighborhood `B` has the same effect as a rectangular box for
mathematical morphology operations.  This may be used to use fast separable
versions of mathematical morphology operations like the van Herk-Gil-Werman
algorithm.

See also: [`LocalFilters.RectangularBox`](@ref).
"""
ismmbox(::RectangularBox) = true
ismmbox(B::Kernel{Bool}) = all(identity, coefs(B))
ismmbox(B::Kernel{T}) where {T<:AbstractFloat} =
    all(x -> x == zero(T), coefs(B))
ismmbox(::Neighborhood) = false

_range(dim::Integer) = _range(Int(dim))

function _range(dim::Int)
    dim ≥ 1 ||
        throw(ArgumentError("neighborhood dimension(s) must be at least one"))
    imin = -(dim >> 1)
    imax = dim + imin - 1
    return imin:imax
end

#------------------------------------------------------------------------------
# METHODS FOR KERNELS

eltype(B::Kernel{T,N}) where {T,N} = T
length(B::Kernel) = length(coefs(B))
size(B::Kernel) = size(coefs(B))
size(B::Kernel, d) = size(coefs(B), d)
getindex(B::Kernel, I::CartesianIndex) = getindex(coefs(B), I + offset(B))
setindex!(B::Kernel, val, I::CartesianIndex) =
    setindex!(coefs(B), val, I + offset(B))

"""

`LocalFilters.coefs(B)` yields the array of coefficients embedded in
kernel `B`.

See also: [`LocalFilters.offset`](@ref).

"""
coefs(B::Kernel) = B.coefs

"""

`LocalFilters.offset(B)` yields the index offset of the array of coefficients
embedded in kernel `B`.   That is, `B[k] ≡ coefs(B)[k + offset(B)]`.

See also: [`LocalFilters.coefs`](@ref).

"""
offset(B::Kernel) = B.offset

# Kernel constructors given a Kernel instance.
Kernel(K::Kernel) = K
Kernel{T}(K::Kernel{T}) where {T} = K
Kernel{T}(K::Kernel{S}) where {S,T} = Kernel(convertcoefs(T, coefs(K)),
                                             initialindex(K))
Kernel{T,N}(K::Kernel{T,N}) where {T,N} = K
Kernel{T,N}(K::Kernel{S,N}) where {S,T,N} = Kernel{T}(K)

# Kernel constructors given a RectangularBox instance.
Kernel(B::RectangularBox{N}) where {N} = Kernel{Bool,N}(B)
Kernel{T}(B::RectangularBox{N}) where {T,N} = Kernel{T,N}(B)
Kernel{T,N}(B::RectangularBox{N}) where {T,N} =
    Kernel{T,N,Array{T,N}}(ones(T, size(B)), initialindex(B))

# Kernel constructors given an array of coefficients.
Kernel(A::AbstractArray) = Kernel(A, defaultstart(A))
Kernel{T}(A::AbstractArray{T}) where {T} = Kernel(A)
Kernel{T}(A::AbstractArray{S}) where {S,T} = Kernel(convertcoefs(T, A))
Kernel{T,N}(A::AbstractArray{T,N}) where {T,N} = Kernel(A)
Kernel{T,N}(A::AbstractArray{S,N}) where {S,T,N} = Kernel{T}(A)

# Kernel constructors given an array of coefficients and starting indices.
Kernel(A::AbstractArray{T,N}, start::CartesianIndex{N}) where {T,N} =
    Kernel{T,N,typeof(A)}(A, start)
Kernel{T}(A::AbstractArray{T,N}, start::CartesianIndex{N}) where {T,N} =
    Kernel{T,N,typeof(A)}(A, start)
Kernel{T,N}(A::AbstractArray{T,N}, start::CartesianIndex{N}) where {T,N} =
    Kernel{T,N,typeof(A)}(A, start)

# Kernel constructors given an array of coefficients and any argument
# suitable to define a Cartesian region.

function Kernel(A::AbstractArray{T,N}, bnds::CartesianRegion{N}) where {T,N}
    # Bounds for indexing the kernel.
    kmin, kmax = limits(bnds)

    # Bounds for indexing the array of coefficients.
    jmin, jmax = limits(A)

    # Check size is identical for all dimensions.
    all(d -> (jmax[d] - jmin[d]) == (kmax[d] - kmin[d]), 1:N) ||
        throw(DimensionMismatch("dimensions must be the same"))

    # Make a kernel with the correct initial index.
    return Kernel(A, kmin)
end

Kernel{T}(A::AbstractArray{T,N}, bnds::CartesianRegion{N}) where {T,N} =
    Kernel(A, bnds)
Kernel{T}(A::AbstractArray{S,N}, bnds::CartesianRegion{N}) where {S,T,N} =
    Kernel(convertcoefs(T, A), bnds)
Kernel{T,N}(A::AbstractArray{S,N}, bnds::CartesianRegion{N}) where {S,T,N} =
    Kernel{T}(A, bnds)

Kernel(A::AbstractArray, rngs::UnitIndexRange...) = Kernel(A, rngs)
Kernel{T}(A::AbstractArray{S}, rngs::UnitIndexRange...) where {S,T} =
    Kernel{T}(A, bnds)
Kernel{T,N}(A::AbstractArray{S,N}, rngs::UnitIndexRange...) where {S,T,N} =
    Kernel{T}(A, rngs)

# Another ways to specify the element type of the kernel coefficients is to
# have their type the first parameter.
Kernel(::Type{T}, B, args...) where {T} = Kernel{T}(B, args...)

# Methods to convert other neighborhoods.  Beware that booleans mean something
# specific, i.e. the result is a so-called *flat* structuring element
# when a kernel whose coefficients are boolean is converted to some
# floating-point type.  See [`convertcoefs`](@ref).

# Make a flat structuring element from a boolean mask.
Kernel(tup::Tuple{<:Any,<:Any}, msk::AbstractArray{Bool}, args...) =
    Kernel(convertcoefs(tup, msk), args...)

Kernel(tup::Tuple{<:Any,<:Any}, B::Kernel{Bool}) =
    Kernel(convertcoefs(tup, coefs(B)), initialindex(B))

# Make a kernel from a function and anything suitable to define a Cartesian
# region.  The element type of the kernel coefficients can be imposed.
Kernel(f::Function, bnds::CartesianRegion) =
    Kernel(map(f, cartesianregion(bnds)), initialindex(bnds))
function Kernel{T}(f::Function, bnds::CartesianRegion{N}) where {T,N}
    kmin, kmax = limits(bnds)
    dims = map(_length, Tuple(kmin), Tuple(kmax))
    W = Array{T,N}(undef, dims)
    offs = initialindex(W) - kmin
    @inbounds for i in cartesianregion(bnds)
        W[i + offs] = f(i)
    end
    return Kernel{T,N,Array{T,N}}(W, kmin)
end
Kernel{T,N}(f::Function, bnds::CartesianRegion{N}) where {T,N} =
    Kernel{T}(f, bnds)

# Conversion of the data type of the kernel coefficients.
for F in (:Float64, :Float32, :Float16)
    @eval begin
        Base.$F(K::Kernel{$F}) = K
        Base.$F(K::Kernel{T}) where {T} = Kernel{$F}(K)
    end
end

"""

    strel(T, A)

yields a *structuring element* suitable for mathematical morphology operations.
The result is a `Kernel` whose elements have type `T` (which can be `Bool` or a
floating-point type).  Argument `A` can be a rectangular box or a `Kernel`
with boolean elements.

If `T` is a floating-point type, then the result is a so-called *flat*
structuring element whose coefficients are `zero(T)` inside the shape defined
by `A` and `-T(Inf)` elsewhere.

See also: [`convertcoefs`](@ref).

"""
strel(::Type{Bool}, K::Kernel{Bool,N}) where {N} = K
strel(::Type{T}, K::Kernel{Bool}) where {T<:AbstractFloat} =
    Kernel(convertcoefs((zero(T), -T(Inf)), coefs(K)), initialindex(K))
strel(::Type{Bool}, B::RectangularBox) = Kernel{Bool}(B)
strel(::Type{T}, B::RectangularBox) where {T<:AbstractFloat} =
    Kernel(zeros(T, size(B)), initialindex(B))

"""

    convertcoefs(T, A)

yields an array of kernel coefficients equivalent to array `A` but whose
elements have type `T`.

If `T` is a floating-point type and `A` is a boolean array, then the values of
the result are `one(T)` where `A` is `true` and `zero(T)` elsewhere.  To use
different values (for instance, to define *flat* *structuring* *elements*), you
may call:

    convertcoefs((vtrue, vfalse), A)

with `A` a boolean array to get an array whose elements are equal to `vtrue`
where `A` is `true` and to `vfalse` otherwise.

See also: [`strel`](@ref).

"""
convertcoefs(::Type{T}, A::AbstractArray{T,N}) where {T,N} = A

function convertcoefs(::Type{T}, A::AbstractArray{S,N}) where {S,T,N}
    B = similar(Array{T,N}, axes(A))
    @inbounds @simd for i in eachindex(A, B)
        B[i] = A[i]
    end
    return B
end

convertcoefs(tup::Tuple{<:Any,<:Any}, A::AbstractArray{Bool}) =
    convertcoefs(promote(tup...), A)

function convertcoefs(tup::Tuple{T,T}, A::AbstractArray{Bool,N}) where {T,N}
    B = similar(Array{T,N}, axes(A))
    vtrue, vfalse = tup
    @inbounds for i in eachindex(A, B)
        B[i] = A[i] ? vtrue : vfalse
    end
    return B
end

"""

    reverse(B::LocalFilters.Neighborhood)

yields neighborhood `B` reversed along all its dimensions.  This can be used
to correlate by `B` rather than convolving by `B`.

"""
reverse(box::RectangularBox) =
    RectangularBox(-finalindex(box), -initialindex(box))

function reverse(ker::Kernel{T,N}) where {T,N}
    kmin = -finalindex(ker)
    kmax = -initialindex(ker)
    rev = Kernel(Array{T}(undef, size(ker)), kmin)
    @inbounds for k in cartesianregion(kmin, kmax)
        rev[k] = ker[-k]
    end
    return rev
end

function strictfloor(::Type{T}, x)::T where {T}
    n = floor(T, x)
    return (n < x ? n : n - one(T))
end

"""

    LocalFilters.ball(N, r)

yields a boolean mask which is a `N`-dimensional array with all dimensions odd
and equal and set to true where position is inside a `N`-dimesional ball of
radius `r`.

"""
function ball(rank::Integer, radius::Real)
    b = radius + 1/2
    r = strictfloor(Int, b)
    dim = 2*r + 1
    dims = ntuple(d->dim, rank)
    arr = Array{Bool}(undef, dims)
    qmax = strictfloor(Int, b^2)
    _ball!(arr, 0, qmax, r, 1:dim, tail(dims))
    return arr
end

@inline function _ball!(arr::AbstractArray{Bool,N},
                        q::Int, qmax::Int, r::Int,
                        range::AbstractUnitRange{Int},
                        dims::Tuple{Int}, I::Int...) where {N}
    nextdims = tail(dims)
    x = -r
    for i in range
        _ball!(arr, q + x*x, qmax, r, range, nextdims, I..., i)
        x += 1
    end
end

@inline function _ball!(arr::AbstractArray{Bool,N},
                        q::Int, qmax::Int, r::Int,
                        range::AbstractUnitRange{Int},
                        ::Tuple{}, I::Int...) where {N}
    x = -r
    for i in range
        arr[I...,i] = (q + x*x ≤ qmax)
        x += 1
    end
end