text:
# Div is truncating by default

"""
    div(x, y, r::RoundingMode=RoundToZero)

Compute the remainder of `x` after integer division by `y`, with the quotient rounded
according to the rounding mode `r`. In other words, the quantity

    y*round(x/y,r)

without any intermediate rounding.

See also: [`fld`](@ref), [`cld`](@ref) which are special cases of this function

# Examples:
```jldoctest
julia> div(4, 3, RoundDown) # Matches fld(4, 3)
1
julia> div(4, 3, RoundUp) # Matches cld(4, 3)
2
julia> div(5, 2, RoundNearest)
2
julia> div(5, 2, RoundNearestTiesAway)
3
julia> div(-5, 2, RoundNearest)
-2
julia> div(-5, 2, RoundNearestTiesAway)
-3
julia> div(-5, 2, RoundNearestTiesUp)
-2
```
"""
div(x, y, r::RoundingMode)

div(a, b) = div(a, b, RoundToZero)

"""
    rem(x, y, r::RoundingMode=RoundToZero)

Compute the remainder of `x` after integer division by `y`, with the quotient rounded
according to the rounding mode `r`. In other words, the quantity

    x - y*round(x/y,r)

without any intermediate rounding.

- if `r == RoundNearest`, then the result is exact, and in the interval
  ``[-|y|/2, |y|/2]``. See also [`RoundNearest`](@ref).

- if `r == RoundToZero` (default), then the result is exact, and in the interval
  ``[0, |y|)`` if `x` is positive, or ``(-|y|, 0]`` otherwise. See also [`RoundToZero`](@ref).

- if `r == RoundDown`, then the result is in the interval ``[0, y)`` if `y` is positive, or
  ``(y, 0]`` otherwise. The result may not be exact if `x` and `y` have different signs, and
  `abs(x) < abs(y)`. See also[`RoundDown`](@ref).

- if `r == RoundUp`, then the result is in the interval `(-y,0]` if `y` is positive, or
  `[0,-y)` otherwise. The result may not be exact if `x` and `y` have the same sign, and
  `abs(x) < abs(y)`. See also [`RoundUp`](@ref).

"""
rem(x, y, r::RoundingMode)

# TODO: Make these primitive and have the two-argument version call these
rem(x, y, ::RoundingMode{:ToZero}) = rem(x,y)
rem(x, y, ::RoundingMode{:Down}) = mod(x,y)
rem(x, y, ::RoundingMode{:Up}) = mod(x,-y)

"""
    fld(x, y)

Largest integer less than or equal to `x/y`. Equivalent to `div(x, y, RoundDown)`.

See also: [`div`](@ref)

# Examples
```jldoctest
julia> fld(7.3,5.5)
1.0
```
"""
fld(a, b) = div(a, b, RoundDown)

"""
    cld(x, y)

Smallest integer larger than or equal to `x/y`. Equivalent to `div(x, y, RoundUp)`.

See also: [`div`](@ref)

# Examples
```jldoctest
julia> cld(5.5,2.2)
3.0
```
"""
cld(a, b) = div(a, b, RoundUp)

# divrem
"""
    divrem(x, y, r::RoundingMode=RoundToZero)

The quotient and remainder from Euclidean division.
Equivalent to `(div(x,y,r), rem(x,y,r))`. Equivalently, with the the default
value of `r`, this call is equivalent to `(x÷y, x%y)`.

# Examples
```jldoctest
julia> divrem(3,7)
(0, 3)

julia> divrem(7,3)
(2, 1)
```
"""
divrem(x, y) = divrem(x, y, RoundToZero)
divrem(a, b, r::RoundingMode) = (div(a, b, r), rem(a, b, r))
function divrem(x::Integer, y::Integer, rnd::typeof(RoundNearest))
    (q, r) = divrem(x, y)
    if x >= 0
        if y >= 0
            r >=        (y÷2) + (isodd(y) | iseven(q)) ? (q+true, r-y) : (q, r)
        else
            r >=       -(y÷2) + (isodd(y) | iseven(q)) ? (q-true, r+y) : (q, r)
        end
    else
        if y >= 0
            r <= -signed(y÷2) - (isodd(y) | iseven(q)) ? (q-true, r+y) : (q, r)
        else
            r <=        (y÷2) - (isodd(y) | iseven(q)) ? (q+true, r-y) : (q, r)
        end
    end
end
function divrem(x::Integer, y::Integer, rnd:: typeof(RoundNearestTiesAway))
    (q, r) = divrem(x, y)
    if x >= 0
        if y >= 0
            r >=        (y÷2) + isodd(y) ? (q+true, r-y) : (q, r)
        else
            r >=       -(y÷2) + isodd(y) ? (q-true, r+y) : (q, r)
        end
    else
        if y >= 0
            r <= -signed(y÷2) - isodd(y) ? (q-true, r+y) : (q, r)
        else
            r <=        (y÷2) - isodd(y) ? (q+true, r-y) : (q, r)
        end
    end
end
function divrem(x::Integer, y::Integer, rnd::typeof(RoundNearestTiesUp))
    (q, r) = divrem(x, y)
    if x >= 0
        if y >= 0
            r >=        (y÷2) + isodd(y) ? (q+true, r-y) : (q, r)
        else
            r >=       -(y÷2) + true     ? (q-true, r+y) : (q, r)
        end
    else
        if y >= 0
            r <= -signed(y÷2) - true     ? (q-true, r+y) : (q, r)
        else
            r <=        (y÷2) - isodd(y) ? (q+true, r-y) : (q, r)
        end
    end
end

"""
    fldmod(x, y)

The floored quotient and modulus after division. A convenience wrapper for
`divrem(x, y, RoundDown)`. Equivalent to `(fld(x,y), mod(x,y))`.
"""
fldmod(x,y) = divrem(x, y, RoundDown)

# We definite generic rounding methods for other rounding modes in terms of
# RoundToZero.
function div(x::Signed, y::Unsigned, ::typeof(RoundDown))
    (q, r) = divrem(x, y)
    q - (signbit(x) & (r != 0))
end
function div(x::Unsigned, y::Signed, ::typeof(RoundDown))
    (q, r) = divrem(x, y)
    q - (signbit(y) & (r != 0))
end

function div(x::Signed, y::Unsigned, ::typeof(RoundUp))
    (q, r) = divrem(x, y)
    q + (!signbit(x) & (r != 0))
end
function div(x::Unsigned, y::Signed, ::typeof(RoundUp))
    (q, r) = divrem(x, y)
    q + (!signbit(y) & (r != 0))
end

function div(x::Integer, y::Integer, rnd::Union{typeof(RoundNearest),
                                              typeof(RoundNearestTiesAway),
                                              typeof(RoundNearestTiesUp)})
    divrem(x,y,rnd)[1]
end

# For bootstrapping purposes, we define div for integers directly. Provide the
# generic signature also
div(a::T, b::T, ::typeof(RoundToZero)) where {T<:Union{BitSigned, BitUnsigned64}} = div(a, b)
div(a::Bool, b::Bool, r::RoundingMode) = div(a, b)
# Prevent ambiguities
for rm in (RoundUp, RoundDown, RoundToZero)
    @eval div(a::Bool, b::Bool, r::$(typeof(rm))) = div(a, b)
end
function div(x::Bool, y::Bool, rnd::Union{typeof(RoundNearest),
                                        typeof(RoundNearestTiesAway),
                                        typeof(RoundNearestTiesUp)})
    div(x, y)
end
fld(a::T, b::T) where {T<:Union{Integer,AbstractFloat}} = div(a, b, RoundDown)
cld(a::T, b::T) where {T<:Union{Integer,AbstractFloat}} = div(a, b, RoundUp)
div(a::Int128, b::Int128, ::typeof(RoundToZero)) = div(a, b)
div(a::UInt128, b::UInt128, ::typeof(RoundToZero)) = div(a, b)
rem(a::Int128, b::Int128, ::typeof(RoundToZero)) = rem(a, b)
rem(a::UInt128, b::UInt128, ::typeof(RoundToZero)) = rem(a, b)

# These are kept for compatibility with external packages overriding fld/cld.
# In 2.0, packages should extend div(a,b,r) instead, in which case, these can
# be removed.
fld(x::Real, y::Real) = div(promote(x,y)..., RoundDown)
cld(x::Real, y::Real) = div(promote(x,y)..., RoundUp)
fld(x::Signed, y::Unsigned) = div(x, y, RoundDown)
fld(x::Unsigned, y::Signed) = div(x, y, RoundDown)
cld(x::Signed, y::Unsigned) = div(x, y, RoundUp)
cld(x::Unsigned, y::Signed) = div(x, y, RoundUp)
fld(x::T, y::T) where {T<:Real} = throw(MethodError(div, (x, y, RoundDown)))
cld(x::T, y::T) where {T<:Real} = throw(MethodError(div, (x, y, RoundUp)))

# Promotion
function div(x::Real, y::Real, r::RoundingMode)
    typeof(x) === typeof(y) && throw(MethodError(div, (x, y, r)))
    if r == RoundToZero
        # For compat. Remove in 2.0.
        div(promote(x, y)...)
    else
        div(promote(x, y)..., r)
    end
end

# Integers
# fld(x,y) == div(x,y) - ((x>=0) != (y>=0) && rem(x,y) != 0 ? 1 : 0)
div(x::T, y::T, ::typeof(RoundDown)) where {T<:Unsigned} = div(x,y)
function div(x::T, y::T, ::typeof(RoundDown)) where T<:Integer
    d = div(x, y, RoundToZero)
    return d - (signbit(x ⊻ y) & (d * y != x))
end

# cld(x,y) = div(x,y) + ((x>0) == (y>0) && rem(x,y) != 0 ? 1 : 0)
function div(x::T, y::T, ::typeof(RoundUp)) where T<:Unsigned
    d = div(x, y, RoundToZero)
    return d + (d * y != x)
end
function div(x::T, y::T, ::typeof(RoundUp)) where T<:Integer
    d = div(x, y, RoundToZero)
    return d + (((x > 0) == (y > 0)) & (d * y != x))
end

# Real
# NOTE: C89 fmod() and x87 FPREM implicitly provide truncating float division,
# so it is used here as the basis of float div().
div(x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} = convert(T,round((x-rem(x,y,r))/y))
rem(x::T, y::T, ::typeof(RoundUp)) where {T<:AbstractFloat} = convert(T,x-y*ceil(x/y))

max_stars_repo_name: derijkp/julia

max_issues_repo_name: derijkp/julia

max_forks_repo_name: derijkp/julia


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.SymmetricSurface
using GeometricExamples

const Δt     = 5.
const ntime  = 2000
const nx     = 100
const ny     = 100
const nsave  = 4
const nplot  = 4
const run_id = "poincare_2nd_symmetric_vprk_dt5"

include("../guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_symmetric_poincare_invariant_2nd.jl")

ext: jl


text:
import Base: view

using Base: ScalarIndex, ViewIndex, Slice, @boundscheck,
            to_indices, compute_offset1, unsafe_length, _maybe_reshape_parent, index_ndims

## traits and properties

struct Contiguous end
struct NonContiguous end

# Detect whether the view is contiguous or not
CuIndexStyle() = Contiguous()
CuIndexStyle(I...) = NonContiguous()
CuIndexStyle(i1::Colon, ::ScalarIndex...) = Contiguous()
CuIndexStyle(i1::AbstractUnitRange, ::ScalarIndex...) = Contiguous()
CuIndexStyle(i1::Colon, I...) = CuIndexStyle(I...)

cuviewlength() = ()
@inline cuviewlength(::Real, I...) = cuviewlength(I...) # skip scalar
@inline cuviewlength(i1::AbstractUnitRange, I...) = (unsafe_length(i1), cuviewlength(I...)...)
@inline cuviewlength(i1::AbstractUnitRange, ::ScalarIndex...) = (unsafe_length(i1),)

## construction

@inline function view(A::CuArray, I::Vararg{Any,N}) where {N}
    J = to_indices(A, I)
    @boundscheck begin
        # Base's boundscheck accesses the indices, so make sure they reside on the CPU.
        # this is expensive, but it's a bounds check after all.
        J_cpu = map(j->adapt(Array, j), J)
        checkbounds(A, J_cpu...)
    end
    J_gpu = map(j->adapt(CuArray, j), J)
    unsafe_view(A, J_gpu, CuIndexStyle(I...))
end

# for contiguous views just return a new CuArray
@inline function unsafe_view(A, I, ::Contiguous)
    unsafe_contiguous_view(_maybe_reshape_parent(A, index_ndims(I...)), I, cuviewlength(I...))
end
@inline function unsafe_contiguous_view(A::CuArray{T}, I::NTuple{N,ViewIndex}, dims::NTuple{M,Integer}) where {T,N,M}
    offset = compute_offset1(A, 1, I) * sizeof(T)
    CuArray{T,M}(pointer(A) + offset, dims, A)
end

# fallback to SubArray when the view is not contiguous
@inline function unsafe_view(A, I, ::NonContiguous)
    Base.unsafe_view(_maybe_reshape_parent(A, index_ndims(I...)), I...)
end

## operations

# upload the SubArray indices when adapting to the GPU
# (can't do this eagerly or the view constructor wouldn't be able to boundscheck)
Adapt.adapt_structure(to::Adaptor, A::SubArray) =
    SubArray(adapt(to, parent(A)), adapt(to, adapt(CuArray, parentindices(A))))


text:
export RNN, rnninit, rnnforw, rnnparam, rnnparams # TODO: we shouldn't export structs like RNN from ops
import Base: show
using AutoGrad: Param, value # TODO: ops should not use Param
using Knet: atype

"""
    rnn = RNN(inputSize, hiddenSize; opts...)
    rnn(x; batchSizes) => y
    rnn.h, rnn.c  # hidden and cell states

`RNN` returns a callable RNN object `rnn`. Given a minibatch of sequences `x`, `rnn(x)`
returns `y`, the hidden states of the final layer for each time step. `rnn.h` and `rnn.c`
fields can be used to set the initial hidden states and read the final hidden states of all
layers.  Note that the final time step of `y` always contains the final hidden state of the
last layer, equivalent to `rnn.h` for a single layer network.

**Dimensions:** The input `x` can be 1, 2, or 3 dimensional and `y` will have the same number
of dimensions as `x`. size(x)=(X,[B,T]) and size(y)=(H/2H,[B,T]) where X is inputSize, B is
batchSize, T is seqLength, H is hiddenSize, 2H is for bidirectional RNNs. By default a 1-D `x`
represents a single instance for a single time step, a 2-D `x` represents a single minibatch
for a single time step, and a 3-D `x` represents a sequence of identically sized minibatches
for multiple time steps. The output `y` gives the hidden state (of the final layer for
multi-layer RNNs) for each time step. The fields `rnn.h` and `rnn.c` represent the hidden
states of all layers in a single time step and have size (H,B,L/2L) where L is numLayers and
2L is for bidirectional RNNs.

**batchSizes:** If `batchSizes=nothing` (default), all sequences in a minibatch are assumed to
be the same length. If `batchSizes` is an array of (non-increasing) integers, it gives us the
batch size for each time step (allowing different sequences in the minibatch to have different
lengths). In this case `x` will typically be 2-D with the second dimension representing
variable size batches for time steps. If `batchSizes` is used, `sum(batchSizes)` should equal
`length(x) ÷ size(x,1)`. When the batch size is different in every time step, hidden states
will have size (H,B,L/2L) where B is always the size of the first (largest) minibatch.

**Hidden states:** The hidden and cell states are kept in `rnn.h` and `rnn.c` fields (the cell
state is only used by LSTM). They can be initialized during construction using the `h` and `c`
keyword arguments, or modified later by direct assignment. Valid values are `nothing`
(default), `0`, or an array of the right type and size possibly wrapped in a `Param`. If the
value is `nothing` the initial state is assumed to be zero and the final state is discarded
keeping the value `nothing`. If the value is `0` the initial state is assumed to be zero and
`0` is replaced by the final state on return. If the value is a valid state, it is used as the
initial state and is replaced by the final state on return.

In a differentiation context the returned final hidden states will be wrapped in `Result`
types. This is necessary if the same RNN object is to be called multiple times in a single
iteration. Between iterations (i.e. after diff/update) the hidden states need to be unboxed
with e.g. `rnn.h = value(rnn.h)` to prevent spurious dependencies. This happens automatically
during the backward pass for GPU RNNs but needs to be done manually for CPU RNNs. See the
[CharLM Tutorial](https://github.com/denizyuret/Knet.jl/blob/master/tutorial/80.charlm.ipynb)
for an example.

**Keyword arguments for RNN:**
- `h=nothing`: Initial hidden state.
- `c=nothing`: Initial cell state.
- `rnnType=:lstm` Type of RNN: One of :relu, :tanh, :lstm, :gru.
- `numLayers=1`: Number of RNN layers.
- `bidirectional=false`: Create a bidirectional RNN if `true`.
- `dropout=0`: Dropout probability. Applied to input and between layers.
- `skipInput=false`: Do not multiply the input with a matrix if `true`.
- `algo=0`: Algorithm to use, see CUDNN docs for details.
- `seed=0`: Random number seed for dropout. Uses `time()` if 0.
- `winit=xavier`: Weight initialization method for matrices.
- `binit=zeros`: Weight initialization method for bias vectors.
- `finit=ones`: Weight initialization method for the bias of forget gates.
- `atype=Knet.atype()`: array type for model weights.

**Formulas:** RNNs compute the output h[t] for a given iteration from the recurrent input
h[t-1] and the previous layer input x[t] given matrices W, R and biases bW, bR from the
following equations:

`:relu` and `:tanh`: Single gate RNN with activation function f:

    h[t] = f(W * x[t] .+ R * h[t-1] .+ bW .+ bR)

`:gru`: Gated recurrent unit:

    i[t] = sigm(Wi * x[t] .+ Ri * h[t-1] .+ bWi .+ bRi) # input gate
    r[t] = sigm(Wr * x[t] .+ Rr * h[t-1] .+ bWr .+ bRr) # reset gate
    n[t] = tanh(Wn * x[t] .+ r[t] .* (Rn * h[t-1] .+ bRn) .+ bWn) # new gate
    h[t] = (1 - i[t]) .* n[t] .+ i[t] .* h[t-1]

`:lstm`: Long short term memory unit with no peephole connections:

    i[t] = sigm(Wi * x[t] .+ Ri * h[t-1] .+ bWi .+ bRi) # input gate
    f[t] = sigm(Wf * x[t] .+ Rf * h[t-1] .+ bWf .+ bRf) # forget gate
    o[t] = sigm(Wo * x[t] .+ Ro * h[t-1] .+ bWo .+ bRo) # output gate
    n[t] = tanh(Wn * x[t] .+ Rn * h[t-1] .+ bWn .+ bRn) # new gate
    c[t] = f[t] .* c[t-1] .+ i[t] .* n[t]               # cell output
    h[t] = o[t] .* tanh(c[t])
"""
mutable struct RNN
    w
    h
    c
    inputSize::Cint
    hiddenSize::Cint
    numLayers::Cint
    dropout::Float64
    seed::Culonglong
    inputMode::Cint    # CUDNN_LINEAR_INPUT = 0, CUDNN_SKIP_INPUT = 1
    direction::Cint    # CUDNN_UNIDIRECTIONAL = 0, CUDNN_BIDIRECTIONAL = 1
    mode::Cint         # CUDNN_RNN_RELU = 0, CUDNN_RNN_TANH = 1, CUDNN_LSTM = 2, CUDNN_GRU = 3
    algo::Cint         # CUDNN_RNN_ALGO_STANDARD = 0, CUDNN_RNN_ALGO_PERSIST_STATIC = 1, CUDNN_RNN_ALGO_PERSIST_DYNAMIC = 2
    dataType::DataType # CUDNN_DATA_FLOAT  = 0, CUDNN_DATA_DOUBLE = 1, CUDNN_DATA_HALF   = 2
    rnnDesc
    dropoutDesc
    dx
    dhx
    dcx
end

function RNN(inputSize, hiddenSize; h=nothing, c=nothing,
             rnnType=:lstm,       # CUDNN_RNN_RELU = 0, CUDNN_RNN_TANH = 1, CUDNN_LSTM = 2, CUDNN_GRU = 3
             bidirectional=false, # CUDNN_UNIDIRECTIONAL = 0, CUDNN_BIDIRECTIONAL = 1
             skipInput=false,     # CUDNN_LINEAR_INPUT = 0, CUDNN_SKIP_INPUT = 1
             numLayers=1,
             dropout=0.0,
             winit=xavier,
             binit=zeros,
             finit=ones,          # forget bias for lstm
             algo=0,              # CUDNN_RNN_ALGO_STANDARD = 0, CUDNN_RNN_ALGO_PERSIST_STATIC = 1, CUDNN_RNN_ALGO_PERSIST_DYNAMIC = 2
             seed=0,              # seed=0 for random init, positive integer for replicability
             atype=atype(),
             # deprecated
             dataType=nothing,    # CUDNN_DATA_FLOAT  = 0, CUDNN_DATA_DOUBLE = 1, CUDNN_DATA_HALF   = 2
             usegpu=nothing,
             handle=nothing,
             )
    if handle !== nothing
        @warn "The handle option is deprecated, GPU implementation using CUDNN.handle()" maxlog=1
    end
    if dataType !== nothing || usegpu !== nothing
        @warn "dataType and usegpu options are deprecated, using atype=$atype" maxlog=1
    end
    dataType = eltype(atype)
    usegpu = !(atype <: Array)

    w = dx = dhx = dcx = nothing
    inputSize = Cint(inputSize)
    hiddenSize = Cint(hiddenSize)
    numLayers = Cint(numLayers)
    dropout = Float64(dropout)
    seed = Culonglong(seed)
    inputMode = Cint(skipInput ? 1 : 0)
    direction = Cint(bidirectional ? 1 : 0)
    mode = findfirst(isequal(rnnType), (:relu,:tanh,:lstm,:gru))
    @assert mode !== nothing "RNN: Valid modes are :relu,:tanh,:lstm,:gru"
    mode = Cint(mode - 1)
    algo = Cint(algo)
    # These should be set by the gpu implementations during first rnnforw:
    dropoutDesc = nothing # usegpu ? DD(handle=handle,dropout=dropout,seed=seed) : nothing # Need to keep dropoutDesc in RNN so it does not get gc'ed.
    rnnDesc = nothing # usegpu ? RD(hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,dataType) : nothing
    r = RNN(w,h,c,inputSize,hiddenSize,numLayers,dropout,seed,inputMode,direction,mode,algo,dataType,rnnDesc,dropoutDesc,dx,dhx,dcx)

    r.w = Array{dataType}(undef,1,1,getRNNParamsSize(r))
    for a in rnnparams(r; useview=true)
        if a == nothing
            continue
        elseif ndims(a) == 2
            copyto!(a, winit(dataType, size(a)))
        elseif ndims(a) == 1
            copyto!(a, binit(dataType, size(a)))
        else
            error("Invalid RNN param $(summary(a))")
        end
    end
    if rnnType == :lstm         # separate initialization for lstm forget biases
        for layer in 1:(numLayers*(bidirectional ? 2 : 1)), id in (2,6), param in (2,)
            a = rnnparam(r, layer, id, param, useview=true)
            if a != nothing
                copyto!(a, finit(dataType, size(a)))
            end
        end
    end
    # get rid of dimension keeping array/element type in case someone hands us atype=typeof(x)=KnetArray{Float32,2}
    if atype isa DataType; atype = (atype.name.wrapper){dataType}; end
    # many copyto! ops to gpu is expensive (~20s), so we init on cpu and copy it over once here.
    r.w = convert(atype, r.w)
    r.w = Param(r.w)
    return r
end

"""
    rnnparam(r::RNN, layer, id, param)

Return a single weight matrix or bias vector as a slice of RNN weights.

Valid `layer` values:
* For unidirectional RNNs 1:numLayers
* For bidirectional RNNs 1:2*numLayers, forw and back layers alternate.

Valid `id` values:
* For RELU and TANH RNNs, input = 1, hidden = 2.
* For GRU reset = 1,4; update = 2,5; newmem = 3,6; 1:3 for input, 4:6 for hidden
* For LSTM inputgate = 1,5; forget = 2,6; newmem = 3,7; output = 4,8; 1:4 for input, 5:8 for hidden

Valid `param` values:
* Return the weight matrix (transposed!) if `param==1`.
* Return the bias vector if `param==2`.

The effect of skipInput: Let I=1 for RELU/TANH, 1:3 for GRU, 1:4 for LSTM
* For skipInput=false (default), rnnparam(r,1,I,1) is a (inputSize,hiddenSize) matrix.
* For skipInput=true, rnnparam(r,1,I,1) is `nothing`.
* For bidirectional, the same applies to rnnparam(r,2,I,1): the first back layer.
* The input biases (par=2) are returned even if skipInput=true.
"""
function rnnparam(r::RNN, layer::Integer, id::Integer, par::Integer; useview=true)
    params_are_good =
    ((1 <= par <= 2) &&
     ((r.direction == 0 && 1 <= layer <= r.numLayers) ||
      (r.direction == 1 && 1 <= layer <= 2*r.numLayers)) &&
     ((r.mode == 0 && 1 <= id <= 2) ||
      (r.mode == 1 && 1 <= id <= 2) ||
      (r.mode == 2 && 1 <= id <= 8) ||
      (r.mode == 3 && 1 <= id <= 6)))
    params_are_good || throw(ArgumentError("Bad arguments for rnnparam, please see `@doc rnnparam`."))
    should_return_nothing =
        ((r.inputMode == 1) &&
         (par == 1) &&
         ((r.mode == 0 && id == 1) ||
          (r.mode == 1 && id == 1) ||
          (r.mode == 2 && 1 <= id <= 4) ||
          (r.mode == 3 && 1 <= id <= 3)) &&
         ((layer == 1) ||
          (r.direction == 1 && layer == 2)))

    i1 = i2 = len = 0
    w = value(r.w)
    # guess i1,i2,len from layer,id,par and rnn specs
    ids = rnnids(r)
    if par == 1 # matrix
        # input matrices are (inputSize,hiddenSize) all others (hiddenSize,hiddenSize)
        # if inputMode==1 then input matrices are empty
        # I=1 for RELU/TANH, 1:3 for GRU, 1:4 for LSTM are input matrices with L=1 for uni and L=1,2 for bi
        inputLayer(r,l)=(l==1 || (l==2 && r.direction==1))
        X, H = r.inputSize, r.hiddenSize
        XH, HH = X*H, H*H
        inputIds = div(ids,2)
        for l = 1:layer, i = 1:ids
            if inputLayer(r,l) && i <= inputIds
                len = (r.inputMode==0 ? XH : 0)
            elseif l>2 && r.direction==1 && i <= div(ids, 2)
                # bidirectional weight
                len = 2HH
            else
                len = HH
            end
            i2 += len
            if l==layer && i==id; break; end
        end
        i1 = i2 - len + 1
    else # bias
        # all biases are length=hidden and there are always numLayers * numIds of them
        len = r.hiddenSize
        layers = r.numLayers * (r.direction == 1 ? 2 : 1)
        i1 = 1 + length(w) - layers * ids * len + ((id-1) + ids * (layer-1)) * len
        i2 = i1 + len - 1
    end
    if i1 > i2
        @assert should_return_nothing
        nothing
    elseif par == 1 # matrix; weights are transposed
        h = Int(r.hiddenSize)
        reshape(view(r.w, i1:i2),:,h)
    else # bias
        view(r.w, i1:i2)
    end
end

function show(io::IO, r::RNN)
    print(io, ("RNNRELU","RNNTANH","LSTM","GRU")[r.mode+1], "(input=", r.inputSize, ",hidden=", r.hiddenSize)
    if r.direction == 1; print(io, ",bidirectional"); end
    if r.numLayers > 1; print(io, ",layers=", r.numLayers); end
    if r.dropout > 0; print(io, ",dropout=", r.dropout); end
    if r.inputMode == 1; print(io, ",skipinput"); end
    if r.dataType != Float32; print(io, ',', r.dataType); end
    print(io, ')')
end

"""
    rnnparams(r::RNN)

Return the RNN parameters as an Array{Any}.

The order of params returned (subject to change):
* All weight matrices come before all bias vectors.
* Matrices and biases are sorted lexically based on (layer,id).
* See @doc rnnparam for valid layer and id values.
* Input multiplying matrices are `nothing` if r.inputMode = 1.
"""
function rnnparams(r::RNN; useview=false)
    layers = r.numLayers * (r.direction == 1 ? 2 : 1)
    ids = rnnids(r)
    ws = []
    for m in (1,2)
        for l in 1:layers
            for i in 1:ids
                push!(ws, rnnparam(r, l, i, m; useview=useview))
            end
        end
    end
    return ws
end

# TODO: instead of rnnparam, rnnparams, have the views represented inside the struct.

rnnids(r) = (r.mode == 2 ? 8 : r.mode == 3 ? 6 : 2)

function getRNNParamsSize(r::RNN)
    whidden = r.hiddenSize * r.hiddenSize
    winput =  r.inputMode == 1 ? 0 : r.hiddenSize * r.inputSize
    bhidden = r.hiddenSize
    binput =  bhidden
    coef = (r.mode == 2 ? 4 : r.mode == 3 ? 3 : 1) * (1 + r.direction)
    nparams = 0
    for i = 1:r.numLayers
        nparams += coef * (whidden + winput + bhidden + binput)
        winput = (1 + r.direction) * whidden
        binput = bhidden
    end
    # TODO: do this check on unit testing
    # if r.rnnDesc != nothing
    #     @assert nparams == cudnnGetRNNParamsSize(r)
    # end
    return nparams
end

function (r::RNN)(x; batchSizes=nothing)
    # Check type/dims of inputs
    XTYPE = typeof(vec(value(x)))
    WTYPE = typeof(vec(value(r.w)))
    # We initialize r.w to Array and then fix it during first rnnforw
    @assert (WTYPE === XTYPE) || (r.rnnDesc === nothing && WTYPE <: Array)
    @assert length(x) > 0
    @assert ndims(x) <= 3
    @assert size(x,1) == r.inputSize
    HSIZE = (r.hiddenSize, batchSizes == nothing ? size(x,2) : batchSizes[1], r.numLayers * (r.direction + 1))
    @assert r.h == nothing || r.h == 0 || (vec(value(r.h)) isa XTYPE && ndims(r.h) <= 3 && (size(r.h,1),size(r.h,2),size(r.h,3)) == HSIZE)
    @assert r.c == nothing || r.c == 0 || (vec(value(r.c)) isa XTYPE && ndims(r.c) <= 3 && (size(r.c,1),size(r.c,2),size(r.c,3)) == HSIZE)
    # apply dropout to input: rnnforw only applies it between layers.
    # TODO: the cpu implementation does not respect the seed parameter.
    # TODO: reconsider dropout for the input in next release.
    x = dropout(x,r.dropout)
    # use hidden inputs unless nothing or 0
    hx = (r.h == 0 ? nothing : r.h)
    cx = (r.c == 0 ? nothing : r.c)
    # produce hidden outputs if hidden inputs != nothing
    hy = (r.h != nothing)
    cy = (r.mode == 2 && r.c != nothing)
    (y, hyout, cyout, rs) = rnnforw(r, r.w, x, hx, cx; hy=hy, cy=cy, batchSizes=batchSizes)
    # In a @diff context hyout and cyout will be Results
    # This is necessary to keep dependencies if r is run several times during one forw pass
    # However unless Results are cleared up before next iteration we will run out of memory
    # I do the clearing in the back function when I know the iteration is done for sure.
    if hy; r.h = hyout; end
    if cy; r.c = cyout; end
    return y
end

# CPU version: need to keep the name rnntest for unit testing
rnnforw(r::RNN, w...; o...) = rnntest(r, w...; o...)

# TODO: interface consistency, deprecate all r,w signatures, have r, signatures. Right now
# rnninit, rnnparam, rnnparams have r,w deprecated. rnnforw is a more difficult case because
# of the optional hx, cx parameters.

function rnntest(r::RNN, ws, x, hx=nothing, cx=nothing;
                 batchSizes=nothing,
                 hy = (hx != nothing),
                 cy = (cx != nothing && r.mode == 2))
    @assert batchSizes == nothing "batchSizes option has not been implemented for the CPU yet."
    @assert value(r.w) === value(ws)
    w = rnnparams(r)
    X,B,T = (size(x,i) for i=1:3) # ndims(x) may be 1,2 or 3
    @assert X == r.inputSize
    Y = Int(r.hiddenSize * (r.direction == 1 ? 2 : 1))
    ysize = ntuple(i->(i==1 ? Y : size(x,i)), ndims(x)) # to match ndims(y) to ndims(x)
    H = Int(r.hiddenSize)
    #@assert (r.inputMode == 0 || H == X)
    L = Int(r.numLayers) * (r.direction == 1 ? 2 : 1)
    hsize = (H, B, L)
    @assert hx == nothing || eqsize(size(hx), hsize)
    @assert cx == nothing || eqsize(size(cx), hsize)
    h = hx==nothing ? fill!(similar(value(x),hsize),0) : hx
    #  hs = Array{Any}[ h[:,:,l] for l=1:L ]
    hs = Array{Any}(undef,L)
    for l = 1:L
        hs[l] = h[:,:,l]
    end
    ys = []
    direction = r.direction
    pdrop = r.dropout
    #=
    All complexity of bidirectional execution
    is packed inside this inline function.
    This causes code repetition, but  works w/o
    touching the existing unidirectional test code
    =#
    @inline bidirect(update_h!) = begin
        xl = x
        for l = 1:(1+direction):L
            skip = l==1 && r.inputMode==1
            hts = []
            if l>1; xl = dropout(xl, pdrop); end
            for t = 1:T
                for (i,ti) in zip([l, l+1], [t, T-t+1])
                    # this function updates h[i]
                    update_h!(xl, i, ti, skip)
                    push!(hts, hs[i])
                end
            end
            # construct the next layer input
            yforw = Array{Any}(hts[1:2:end-1])
            yback = Array{Any}(reverse(hts[2:2:end]))
            ybs = []
            for (yf, yb) in zip(yforw, yback)
                push!(ybs, vcat(yf, yb))
            end
            # now ybs contans (2 * hiddenSize, batchSize) matrices
            # so cat them to add time dimension
            xl = reshape(hcat(ybs...), (2r.hiddenSize, size(x,2), length(ybs)))
        end
        ys = xl
    end

    if r.mode <= 1
        #@assert r.inputMode == 0 || all(w[1:1+r.direction] .== nothing)
        f = r.mode == 0 ? relu : tanh
        if direction == 0
            for t = 1:T
                for l = 1:L
                    xl = l>1 ? dropout(hs[l-1], pdrop) : nothing
                    wx,wh,bx,bh = w[2l-1],w[2l],w[2L+2l-1],w[2L+2l]
                    wxt = (l > 1 ? wx' * xl : r.inputMode==0 ? wx' * x[:,:,t] : x[:,:,t])
                    hs[l] = f.(wxt .+ wh' * hs[l] .+ bx .+ bh)
                end
                push!(ys, hs[L])
            end
        else
            bidirect() do xl, i, ti, skip
                wx,wh,bx,bh = w[2i-1],w[2i],w[2L+2i-1],w[2L+2i]
                wxt =  skip ? xl[:,:,ti] : wx' * xl[:,:,ti]
                hs[i] = f.(wxt .+ wh' * hs[i] .+ bx .+ bh)
            end
        end
    elseif r.mode == 2           # LSTM
        #@assert r.inputMode == 0 || all(w[1:4*(1+r.direction)] .== nothing)
        # it = σ(Wixt + Riht-1 + bWi + bRi)
        # ft = σ(Wfxt + Rfht-1 + bWf + bRf)
        # ot = σ(Woxt + Roht-1 + bWo + bRo)
        # c't = tanh(Wcxt + Rcht-1 + bWc + bRc)
        # ct = ft◦ct-1 + it◦c't
        # ht = ot◦tanh(ct)
        c = cx==nothing ? fill!(similar(value(x),hsize),0) : cx
        cs = Array{Any}(undef,L)
        for l = 1:L
            cs[l] = c[:,:,l]
        end
        if direction == 0
            for t = 1:T
                for l = 1:L
                    xl = l>1 ? dropout(hs[l-1], pdrop) : nothing
                    Wi,Wf,Wc,Wo,Ri,Rf,Rc,Ro = w[1+8*(l-1):8l]
                    bWi,bWf,bWc,bWo,bRi,bRf,bRc,bRo = w[8L+1+8*(l-1):8L+8l]
                    Wixt = (l > 1 ? Wi' * xl : r.inputMode==0 ? Wi' * x[:,:,t] : x[:,:,t])
                    Wfxt = (l > 1 ? Wf' * xl : r.inputMode==0 ? Wf' * x[:,:,t] : x[:,:,t])
                    Wcxt = (l > 1 ? Wc' * xl : r.inputMode==0 ? Wc' * x[:,:,t] : x[:,:,t])
                    Woxt = (l > 1 ? Wo' * xl : r.inputMode==0 ? Wo' * x[:,:,t] : x[:,:,t])
                    it = sigm.(Wixt .+ Ri' * hs[l] .+ bWi .+ bRi)
                    ft = sigm.(Wfxt .+ Rf' * hs[l] .+ bWf .+ bRf)
                    ot = sigm.(Woxt .+ Ro' * hs[l] .+ bWo .+ bRo)
                    cn = tanh.(Wcxt .+ Rc' * hs[l] .+ bWc .+ bRc)
                    cs[l] = ft .* cs[l] .+ it .* cn
                    hs[l] = ot .* tanh.(cs[l])
                end
                push!(ys, hs[L])
            end
        else
            bidirect() do xl, i, ti, skip
                Wi,Wf,Wc,Wo,Ri,Rf,Rc,Ro = w[1+8*(i-1):8i]
                bWi,bWf,bWc,bWo,bRi,bRf,bRc,bRo = w[8L+1+8*(i-1):8L+8i]
                Wixt = skip ? xl[:,:,ti] : Wi' * xl[:,:,ti]
                Wfxt = skip ? xl[:,:,ti] : Wf' * xl[:,:,ti]
                Wcxt = skip ? xl[:,:,ti] : Wc' * xl[:,:,ti]
                Woxt = skip ? xl[:,:,ti] : Wo' * xl[:,:,ti]
                it = sigm.(Wixt .+ Ri' * hs[i] .+ bWi .+ bRi)
                ft = sigm.(Wfxt .+ Rf' * hs[i] .+ bWf .+ bRf)
                ot = sigm.(Woxt .+ Ro' * hs[i] .+ bWo .+ bRo)
                cn = tanh.(Wcxt .+ Rc' * hs[i] .+ bWc .+ bRc)
                cs[i] = ft .* cs[i] .+ it .* cn
                hs[i] = ot .* tanh.(cs[i])
            end
        end
    elseif r.mode == 3           # GRU
        #@assert r.inputMode == 0 || all(w[1:3*(1+r.direction)] .== nothing)
        # rt = σ(Wrxt + Rrht-1 + bWr + bRr)
        # it = σ(Wixt + Riht-1 + bWi + bRu)
        # h't = tanh(Whxt + rt◦(Rhht-1 + bRh) + bWh)
        # ht = (1 - it)◦h't + it◦ht-1
        if direction == 0
            for t = 1:T
                for l = 1:L
                    xl = l>1 ? dropout(hs[l-1], pdrop) : nothing
                    Wr,Wi,Wh,Rr,Ri,Rh = w[1+6*(l-1):6l]
                    bWr,bWi,bWh,bRr,bRi,bRh = w[6L+1+6*(l-1):6L+6l]
                    Wrxt = (l > 1 ? Wr' * xl : r.inputMode==0 ? Wr' * x[:,:,t] : x[:,:,t])
                    Wixt = (l > 1 ? Wi' * xl : r.inputMode==0 ? Wi' * x[:,:,t] : x[:,:,t])
                    Whxt = (l > 1 ? Wh' * xl : r.inputMode==0 ? Wh' * x[:,:,t] : x[:,:,t])
                    rt = sigm.(Wrxt .+ Rr' * hs[l] .+ bWr .+ bRr)
                    it = sigm.(Wixt .+ Ri' * hs[l] .+ bWi .+ bRi)
                    ht = tanh.(Whxt .+ rt .* (Rh' * hs[l] .+ bRh) .+ bWh)
                    hs[l] = (1 .- it) .* ht .+ it .* hs[l]
                end
                push!(ys, hs[L])
            end
        else
            bidirect() do xl, i, ti, skip
                Wr,Wi,Wh,Rr,Ri,Rh = w[1+6*(i-1):6i]
                bWr,bWi,bWh,bRr,bRi,bRh = w[6L+1+6*(i-1):6L+6i]
                Wrxt = skip ? xl[:, :, ti] : Wr' * xl[:, :, ti]
                Wixt = skip ? xl[:, :, ti] : Wi' * xl[:, :, ti]
                Whxt = skip ? xl[:, :, ti] : Wh' * xl[:, :, ti]
                rt = sigm.(Wrxt .+ Rr' * hs[i] .+ bWr .+ bRr)
                it = sigm.(Wixt .+ Ri' * hs[i] .+ bWi .+ bRi)
                ht = tanh.(Whxt .+ rt .* (Rh' * hs[i] .+ bRh) .+ bWh)
                hs[i] = (1 .- it) .* ht .+ it .* hs[i]
            end
        end
    else
        error("RNN not supported")
    end
    y = r.direction == 0 ? reshape(hcat(ys...), ysize) : reshape(ys,ysize)
    hyout = hy ? reshape(hcat(hs...), hsize) : nothing
    cyout = cy && r.mode == 2 ? reshape(hcat(cs...), hsize) : nothing
    return (y,hyout,cyout,nothing,nothing)
end

# compare sizes ignoring trailing ones
function eqsize(a, b)
    na = length(a)
    nb = length(b)
    (na == nb ? a == b : na > nb ?
     a[1:nb] == b && all(a[nb+1:end] .== 1) :
     b[1:na] == a && all(b[na+1:end] .== 1))
end

## TO BE DEPRECATED:
function rnninit(x...; o...)
    # @warn "rnninit is deprecated, use RNN instead" maxlog=1
    r=RNN(x...; o...)
    return (r,value(r.w))
end

function rnnparams(r,w;o...)
    # @warn "rnnparams(r,w) is deprecated, use rnnparams(r) instead" maxlog=1
    @assert value(w)===value(r.w)
    rnnparams(r;o...)
end

function rnnparam(r,w,l,i,d;o...)
    # @warn "rnnparam(r,w,l,i,d) is deprecated, use rnnparam(r,l,i,d) instead" maxlog=1
    @assert value(w)===value(r.w)
    rnnparam(r,l,i,d;o...)
end

function xavier(a...; gain=1)
    w = rand(a...)
    @assert ndims(w) == 2
    fanout, fanin = size(w)
    s = convert(eltype(w), gain*sqrt(6 / (fanin + fanout)))
    return 2s .* w .- s
end


text:
#tablularReasoning
using Calculus, Symbolics, Latexify

mutable struct tabQ
	subQs::Vector{String}
	isHorizontal::Bool
	nFunct::Int
	qType::String
end

qTypes = ["approximation","expression manipulation"]
subApprx = ["mean value theorem","derivative approximation","Riemann sums"]
subExpMani = ["chain rule", "product rule", "quotient rule", "L'Hopital's rule"]
bothTypes = ["Average value","exact derivative","integrals"]

function makeExpr(c::Int,h,subExp::String)
	@variables x
	D = Differential(x)
	baseFunct = ["2^x","x+$c","sin(pi*x)","cos(pi*x)",]
	choice = rand(baseFunct)

end

#determine which data is needed
function table(q::tabQ)
	txt = """ """
	l = rand(UnitRange(-5,3))
	r = l+4
	vals = UnitRange(l,r)

	#determines which parts of question a,b,c,d are which types
	if q.qType == "approximation"
		q.nFunct = 1
		q.isHorizontal = true

	else
		#2-3 questions are expression manip, the others are from both types
		q.nFunct = 3

	end
	#generates display

end


text:
using Images, FactCheck, Base.Test, Colors

global checkboard

facts("Edge") do

    EPS = 1e-14

    ## Checkerboard array, used to test image gradients
    white{T}(::Type{T}) = one(T)
    black{T}(::Type{T}) = zero(T)
    white{T<:Unsigned}(::Type{T}) = typemax(T)
    black{T<:Unsigned}(::Type{T}) = typemin(T)

    function checkerboard{T}(::Type{T}, sq_width::Integer, count::Integer)
        wh = fill(white(T), (sq_width,sq_width))
        bk = fill(black(T), (sq_width,sq_width))
        bw = [wh bk; bk wh]
        vert = repmat(bw, (count>>1), 1)
        isodd(count) && (vert = vcat(vert, [wh bk]))
        cb = repmat(vert, 1, (count>>1))
        isodd(count) && (cb = hcat(cb, vert[:,1:sq_width]))
        cb
    end

    checkerboard(sq_width::Integer, count::Integer) = checkerboard(UInt8, sq_width, count)

    #Canny Edge Detection
    context("Canny Edge Detection") do
        #General Checks
        img = zeros(10, 10)
        edges = canny(img)
        @fact eltype(edges.data) --> Gray{U8}
        @fact all(edges .== 0.0) --> true

        #Box Edges

        img[2:end-1, 2:end-1] = 1
        edges = canny(img)
        @fact all(edges[2:end-1, 2] .== 1.0) --> true
        @fact all(edges[2:end-1, end-1] .== 1.0) --> true
        @fact all(edges[2, 2:end-1] .== 1.0) --> true
        @fact all(edges[end-1, 2:end-1] .== 1.0) --> true
        @fact all(edges[3:end-2, 3:end-2] .== 0.0) --> true

        edges = canny(img, 1.4, 0.9, 0.2, percentile = false)
        @fact all(edges[2:end-1, 2] .== 1.0) --> true
        @fact all(edges[2:end-1, end-1] .== 1.0) --> true
        @fact all(edges[2, 2:end-1] .== 1.0) --> true
        @fact all(edges[end-1, 2:end-1] .== 1.0) --> true
        @fact all(edges[3:end-2, 3:end-2] .== 0.0) --> true

        #Checkerboard - Corners are not detected as Edges!
        img = checkerboard(Gray, 5, 3)
        edges = canny(img, 1.4, 0.8, 0.2)
        id = [1,2,3,4,6,7,8,9,10,12,13,14,15]
        @fact all(edges[id, id] .== zero(Gray)) --> true
        id = [5, 11]
        id2 = [1,2,3,4,7,8,9,12,13,14,15]
        id3 = [5,6,10,11]
        @fact all(edges[id, id2] .== one(Gray)) --> true
        @fact all(edges[id2, id] .== one(Gray)) --> true
        @fact all(edges[id, id3] .== zero(Gray)) --> true
        @fact all(edges[id3, id] .== zero(Gray)) --> true

        #Diagonal Edge
        img = zeros(10,10)
        img[diagind(img)] = 1
        img[diagind(img, 1)] = 1
        img[diagind(img, -1)] = 1
        edges = canny(img)
        @fact all(edges[diagind(edges, 2)] .== 1.0) --> true
        @fact all(edges[diagind(edges, -2)] .== 1.0) --> true
        nondiags = setdiff(1:1:100, union(diagind(edges, 2), diagind(edges, -2)))
        @fact all(edges[nondiags] .== 0.0) --> true

        #Checks Hysteresis Thresholding
        img = ones(Gray{U8}, (10, 10))
        img[3:7, 3:7] = 0.0
        img[4:6, 4:6] = 0.7
        thresholded = Images.hysteresis_thresholding(img, 0.9, 0.8)
        @fact all(thresholded[3:7, 3:7] .== 0.0) --> true
        @fact all(thresholded[1:2, :] .== 0.9) --> true
        @fact all(thresholded[:, 1:2] .== 0.9) --> true
        @fact all(thresholded[8:10, :] .== 0.9) --> true
        @fact all(thresholded[:, 8:10] .== 0.9) --> true

        thresholded = Images.hysteresis_thresholding(img, 0.9, 0.6)
        @fact all(thresholded[4:6, 4:6] .== 0.5) --> true
        @fact all(thresholded[3:7, 3:7] .< 0.9) --> true
        @fact all(thresholded[1:2, :] .== 0.9) --> true
        @fact all(thresholded[:, 1:2] .== 0.9) --> true
        @fact all(thresholded[8:10, :] .== 0.9) --> true
        @fact all(thresholded[:, 8:10] .== 0.9) --> true

        img[3, 5] = 0.7
        thresholded = Images.hysteresis_thresholding(img, 0.9, 0.6)
        @fact all(thresholded[4:6, 4:6] .== 0.9) --> true
        @fact all(thresholded[3:7, 3] .== 0.0) --> true
        @fact all(thresholded[3:7, 7] .== 0.0) --> true
        @fact all(thresholded[7, 3:7] .== 0.0) --> true
        @fact all(vec(thresholded[3, 3:7]) .== [0.0, 0.0, 0.9, 0.0, 0.0]) --> true
        @fact all(thresholded[1:2, :] .== 0.9) --> true
        @fact all(thresholded[:, 1:2] .== 0.9) --> true
        @fact all(thresholded[8:10, :] .== 0.9) --> true
        @fact all(thresholded[:, 8:10] .== 0.9) --> true
    end

    SZ=5
    cb_array    = checkerboard(SZ,3)
    cb_image_xy = grayim(cb_array)
    cb_image_yx = grayim(cb_array)
    cb_image_yx["spatialorder"] = ["y","x"]

    cb_image_gray = convert(Image{Gray}, cb_image_xy)
    cb_image_rgb = convert(Image{RGB}, cb_image_xy)
    cb_image_rgb2 = convert(Image{RGB{Float64}}, cb_image_xy)

    @compat context("Checkerboard") do
        for method in ["sobel", "prewitt", "ando3", "ando4", "ando5", "ando4_sep", "ando5_sep"]
            ## Checkerboard array

            (agx, agy) = imgradients(cb_array, method)
            amag = magnitude(agx, agy)
            agphase = phase(agx, agy)
            @fact (amag, agphase) --> magnitude_phase(agx, agy)

            @fact agx[1,SZ]   --> less_than(0.0)      # white to black transition
            @fact agx[1,2*SZ] --> greater_than(0.0)   # black to white transition
            @fact agy[SZ,1]   --> less_than(0.0)      # white to black transition
            @fact agy[2*SZ,1] --> greater_than(0.0)   # black to white transition

            # Test direction of increasing gradient
            @fact cos(agphase[1,SZ])   - (-1.0) --> less_than(EPS)   # increasing left  (=  pi   radians)
            @fact cos(agphase[1,2*SZ]) -   1.0  --> less_than(EPS)   # increasing right (=   0   radians)
            @fact sin(agphase[SZ,1])   -   1.0  --> less_than(EPS)   # increasing up    (=  pi/2 radians)
            @fact sin(agphase[2*SZ,1]) - (-1.0) --> less_than(EPS)   # increasing down  (= -pi/2 radians)

            # Test that orientation is perpendicular to gradient
            aorient = orientation(agx, agy)
            @fact all((cos.(agphase) .* cos.(aorient) .+
                       sin.(agphase) .* sin.(aorient) .< EPS) |
                      ((agphase .== 0.0) & (aorient .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient

            ## Checkerboard Image with row major order

            (gx, gy) = imgradients(cb_image_xy, method)
            mag = magnitude(gx, gy)
            gphase = phase(gx, gy)
            @fact (mag, gphase) --> magnitude_phase(gx, gy)

            @fact gx[SZ,1]   --> less_than(0.0)      # white to black transition
            @fact gx[2*SZ,1] --> greater_than(0.0)   # black to white transition
            @fact gy[1,SZ]   --> less_than(0.0)      # white to black transition
            @fact gy[1,2*SZ] --> greater_than(0.0)   # black to white transition

            @fact cos(gphase[SZ,1])   - (-1.0) --> less_than(EPS)   # increasing left  (=  pi   radians)
            @fact cos(gphase[2*SZ,1]) -   1.0  --> less_than(EPS)   # increasing right (=   0   radians)
            @fact sin(gphase[1,SZ])   -   1.0  --> less_than(EPS)   # increasing up    (=  pi/2 radians)
            @fact sin(gphase[1,2*SZ]) - (-1.0) --> less_than(EPS)   # increasing down  (= -pi/2 radians)

            # Test that orientation is perpendicular to gradient
            orient = orientation(gx, gy)
            @fact all((cos.(gphase) .* cos.(orient) .+
                       sin.(gphase) .* sin.(orient) .< EPS) |
                      ((gphase .== 0.0) & (orient .== 0.0))) --> true
                       # this part is where both are zero because there is no gradient

            ## Checkerboard Image with column-major order

            (gx, gy) = imgradients(cb_image_yx, method)
            mag = magnitude(gx, gy)
            gphase = phase(gx, gy)
            @fact (mag, gphase) --> magnitude_phase(gx, gy)

            @fact gx[1,SZ]   --> less_than(0.0)      # white to black transition
            @fact gx[1,2*SZ] --> greater_than(0.0)   # black to white transition
            @fact gy[SZ,1]   --> less_than(0.0)      # white to black transition
            @fact gy[2*SZ,1] --> greater_than(0.0)   # black to white transition

            # Test direction of increasing gradient
            @fact cos(gphase[1,SZ])   - (-1.0) --> less_than(EPS)   # increasing left  (=  pi   radians)
            @fact cos(gphase[1,2*SZ]) -   1.0  --> less_than(EPS)   # increasing right (=   0   radians)
            @fact sin(gphase[SZ,1])   -   1.0  --> less_than(EPS)   # increasing up    (=  pi/2 radians)
            @fact sin(gphase[2*SZ,1]) - (-1.0) --> less_than(EPS)   # increasing down  (= -pi/2 radians)

            # Test that orientation is perpendicular to gradient
            orient = orientation(gx, gy)
            @fact all((cos.(gphase) .* cos.(orient) .+
                       sin.(gphase) .* sin.(orient) .< EPS) |
                      ((gphase .== 0.0) & (orient .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient

            ## Checkerboard Image with Gray pixels

            (gxg, gyg) = imgradients(cb_image_gray, method)
            magg = magnitude(gxg, gyg)
            gphaseg = phase(gxg, gyg)
            @fact (magg, gphaseg) --> magnitude_phase(gxg, gyg)

            @fact gxg[SZ,1]   --> less_than(0.0)      # white to black transition
            @fact gxg[2*SZ,1] --> greater_than(0.0)   # black to white transition
            @fact gyg[1,SZ]   --> less_than(0.0)      # white to black transition
            @fact gyg[1,2*SZ] --> greater_than(0.0)   # black to white transition

            @fact cos(gphaseg[SZ,1])   - (-1.0) --> less_than(EPS)   # increasing left  (=  pi   radians)
            @fact cos(gphaseg[2*SZ,1]) -   1.0  --> less_than(EPS)   # increasing right (=   0   radians)
            @fact sin(gphaseg[1,SZ])   -   1.0  --> less_than(EPS)   # increasing up    (=  pi/2 radians)
            @fact sin(gphaseg[1,2*SZ]) - (-1.0) --> less_than(EPS)   # increasing down  (= -pi/2 radians)

            # Test that orientation is perpendicular to gradient
            orientg = orientation(gxg, gyg)
            @fact all((cos.(gphaseg) .* cos.(orientg) .+
                       sin.(gphaseg) .* sin.(orientg) .< EPS) |
                      ((gphaseg .== 0.0) & (orientg .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient

            ## Checkerboard Image with RBG pixels

            (gx_rgb, gy_rgb) = imgradients(cb_image_rgb, method)
            mag_rgb = magnitude(gx_rgb, gy_rgb)
            gphase_rgb = phase(gx_rgb, gy_rgb)
            @fact (mag_rgb, gphase_rgb) --> magnitude_phase(gx_rgb, gy_rgb)

            @fact gx_rgb[1,SZ,1]   --> less_than(0.0)      # white to black transition
            @fact gx_rgb[1,2*SZ,1] --> greater_than(0.0)   # black to white transition
            @fact gy_rgb[1,1,SZ]   --> less_than(0.0)      # white to black transition
            @fact gy_rgb[1,1,2*SZ] --> greater_than(0.0)   # black to white transition

            @fact cos(gphase_rgb[1,SZ,1])   - (-1.0) --> less_than(EPS)   # increasing left  (=  pi   radians)
            @fact cos(gphase_rgb[1,2*SZ,1]) -   1.0  --> less_than(EPS)   # increasing right (=   0   radians)
            @fact sin(gphase_rgb[1,1,SZ])   -   1.0  --> less_than(EPS)   # increasing up    (=  pi/2 radians)
            @fact sin(gphase_rgb[1,1,2*SZ]) - (-1.0) --> less_than(EPS)   # increasing down  (= -pi/2 radians)

            # Test that orientation is perpendicular to gradient
            orient_rgb = orientation(gx_rgb, gy_rgb)
            @fact all((cos.(gphase_rgb) .* cos.(orient_rgb) .+
                       sin.(gphase_rgb) .* sin.(orient_rgb) .< EPS) |
                      ((gphase_rgb .== 0.0) & (orient_rgb .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient

            ## Checkerboard Image with RBG{Float64} pixels

            (gx_rgb, gy_rgb) = imgradients(cb_image_rgb2, method)
            mag_rgb = magnitude(gx_rgb, gy_rgb)
            gphase_rgb = phase(gx_rgb, gy_rgb)
            @fact (mag_rgb, gphase_rgb) --> magnitude_phase(gx_rgb, gy_rgb)

            @fact gx_rgb[1,SZ,1]   --> less_than(0.0)      # white to black transition
            @fact gx_rgb[1,2*SZ,1] --> greater_than(0.0)   # black to white transition
            @fact gy_rgb[1,1,SZ]   --> less_than(0.0)      # white to black transition
            @fact gy_rgb[1,1,2*SZ] --> greater_than(0.0)   # black to white transition

            @fact cos(gphase_rgb[1,SZ,1])   - (-1.0) --> less_than(EPS)   # increasing left  (=  pi   radians)
            @fact cos(gphase_rgb[1,2*SZ,1]) -   1.0  --> less_than(EPS)   # increasing right (=   0   radians)
            @fact sin(gphase_rgb[1,1,SZ])   -   1.0  --> less_than(EPS)   # increasing up    (=  pi/2 radians)
            @fact sin(gphase_rgb[1,1,2*SZ]) - (-1.0) --> less_than(EPS)   # increasing down  (= -pi/2 radians)

            # Test that orientation is perpendicular to gradient
            orient_rgb = orientation(gx_rgb, gy_rgb)
            @fact all((cos.(gphase_rgb) .* cos.(orient_rgb) .+
                       sin.(gphase_rgb) .* sin.(orient_rgb) .< EPS) |
                      ((gphase_rgb .== 0.0) & (orient_rgb .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient
        end
    end

    @compat context("Diagonals") do
        # Create an image with white along diagonals -2:2 and black elsewhere
        m = zeros(UInt8, 20,20)
        for i = -2:2; m[diagind(m,i)] = 0xff; end

        m_xy = grayim(m')
        m_yx = grayim(m)
        m_yx["spatialorder"] = ["y","x"]

        for method in ["sobel", "prewitt", "ando3", "ando4", "ando5", "ando4_sep", "ando5_sep"]
            ## Diagonal array

            (agx, agy) = imgradients(m, method)
            amag = magnitude(agx, agy)
            agphase = phase(agx, agy)
            @fact (amag, agphase) --> magnitude_phase(agx, agy)

            @fact agx[7,9]  --> less_than(0.0)      # white to black transition
            @fact agx[10,8] --> greater_than(0.0)   # black to white transition
            @fact agy[10,8] --> less_than(0.0)      # white to black transition
            @fact agy[7,9]  --> greater_than(0.0)   # black to white transition

            # Test direction of increasing gradient
            @fact abs(agphase[10,8] -    pi/4 ) --> less_than(EPS)   # lower edge (increasing up-right  =   pi/4 radians)
            @fact abs(agphase[7,9]  - (-3pi/4)) --> less_than(EPS)   # upper edge (increasing down-left = -3pi/4 radians)

            # Test that orientation is perpendicular to gradient
            aorient = orientation(agx, agy)
            @fact all((cos.(agphase) .* cos.(aorient) .+
                       sin.(agphase) .* sin.(aorient) .< EPS) |
                      ((agphase .== 0.0) & (aorient .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient

            ## Diagonal Image, row-major order

            (gx, gy) = imgradients(m_xy, method)
            mag = magnitude(gx, gy)
            gphase = phase(gx, gy)
            @fact (mag, gphase) --> magnitude_phase(gx, gy)

            @fact gx[9,7]  --> less_than(0.0)      # white to black transition
            @fact gx[8,10] --> greater_than(0.0)   # black to white transition
            @fact gy[8,10] --> less_than(0.0)      # white to black transition
            @fact gy[9,7]  --> greater_than(0.0)   # black to white transition

            # Test direction of increasing gradient
            @fact abs(gphase[8,10] -    pi/4 ) --> less_than(EPS)   # lower edge (increasing up-right  =   pi/4 radians)
            @fact abs(gphase[9,7]  - (-3pi/4)) --> less_than(EPS)   # upper edge (increasing down-left = -3pi/4 radians)

            # Test that orientation is perpendicular to gradient
            orient = orientation(gx, gy)
            @fact all((cos.(gphase) .* cos.(orient) .+
                       sin.(gphase) .* sin.(orient) .< EPS) |
                      ((gphase .== 0.0) & (orient .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient

            ## Diagonal Image, column-major order

            (gx, gy) = imgradients(m_yx, method)
            mag = magnitude(gx, gy)
            gphase = phase(gx, gy)
            @fact (mag, gphase) --> magnitude_phase(gx, gy)

            @fact gx[7,9]  --> less_than(0.0)      # white to black transition
            @fact gx[10,8] --> greater_than(0.0)   # black to white transition
            @fact gy[10,8] --> less_than(0.0)      # white to black transition
            @fact gy[7,9]  --> greater_than(0.0)   # black to white transition

            # Test direction of increasing gradient
            @fact abs(gphase[10,8] -    pi/4 ) --> less_than(EPS)   # lower edge (increasing up-right  =   pi/4 radians)
            @fact abs(gphase[7,9]  - (-3pi/4)) --> less_than(EPS)   # upper edge (increasing down-left = -3pi/4 radians)

            # Test that orientation is perpendicular to gradient
            orient = orientation(gx, gy)
            @fact all((cos.(gphase) .* cos.(orient) .+
                       sin.(gphase) .* sin.(orient) .< EPS) |
                      ((gphase .== 0.0) & (orient .== 0.0))) --> true
                      # this part is where both are zero because there is no gradient
        end
    end

    # Nonmaximal suppression

    function thin_edges(img)
        # Get orientation
        gx,gy = imgradients(img)
        orient = phase(gx,gy)

        # Do NMS thinning
        thin_edges_nonmaxsup_subpix(img, orient, radius=1.35)
    end

    function nms_test_horiz_vert(img, which)
        ## which = :horizontal or :vertical

        # Do NMS thinning
        t,s = thin_edges(img)

        # Calc peak location by hand

        # Interpolate values 1.35 pixels left and right
        # Orientation is zero radians -> to the right
        v1 = 6 - 0.35   # slope on right is 1
        v2 = 5 - 0.35*2 # slope on left is 2
        c = 7.0         # peak value

        # solve v = a*r^2 + b*r + c
        a = (v1 + v2)/2 - c
        b = a + c - v2
        r = -b/2a

        @fact abs(r - 1/6) --> less_than(EPS)

        # Location and value at peak
        peakloc = r*1.35 + 3
        peakval = a*r^2 + b*r + c

        transposed = spatialorder(img)[1] == "x"
        horizontal = which == :horizontal

        test_axis1 = transposed $ !horizontal

        if test_axis1
            @fact all(t[:,[1,2,4,5]] .== 0) --> true
            @fact all(t[:,3]   .== peakval) --> true
            @fact all(s[:,[1,2,4,5]] .== zero(Graphics.Point)) --> true
        else
            @fact all(t[[1,2,4,5],:] .== 0) --> true
            @fact all(t[3,:]   .== peakval) --> true
            @fact all(s[[1,2,4,5],:] .== zero(Graphics.Point)) --> true
        end

        if transposed
            if which == :horizontal
                @fact     [pt.x for pt in s[:,3]]  --> [1:5;]
                @fact all([pt.y for pt in s[:,3]] .== peakloc) --> true
            else
                @fact all([pt.x for pt in s[3,:]] .== peakloc) --> true
                @fact     [pt.y for pt in s[3,:]]  --> [1:5;]
            end
        else
            if which == :horizontal
                @fact     [pt.x for pt in s[3,:]]  --> [1:5;]
                @fact all([pt.y for pt in s[3,:]] .== peakloc) --> true
            else
                @fact all([pt.x for pt in s[:,3]] .== peakloc) --> true
                @fact     [pt.y for pt in s[:,3]]  --> [1:5;]
            end
        end
    end

    context("Nonmax suppression vertical edge") do
        # Test image: vertical edge
        m = [3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0]

        m_xy = grayim(m')
        m_yx = grayim(m)
        m_yx["spatialorder"] = ["y","x"]

        nms_test_horiz_vert(m, :vertical)
        nms_test_horiz_vert(m_xy, :vertical)
        nms_test_horiz_vert(m_yx, :vertical)
    end

    context("Nonmax suppression horizontal edge") do
        # Test image: horizontal edge
        m = [3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0
             3.0  5.0  7.0  6.0  5.0]
        m = m'
        m_xy = grayim(m')
        m_yx = grayim(m)
        m_yx["spatialorder"] = ["y","x"]

        nms_test_horiz_vert(m, :horizontal)
        nms_test_horiz_vert(m_xy, :horizontal)
        nms_test_horiz_vert(m_yx, :horizontal)
    end

function nms_test_diagonal(img)
    # Do NMS thinning
    t,s = thin_edges(img)

    # Calc peak location by hand

    # Interpolate values 1.35 pixels up and left, down and right
    # using bilinear interpolation
    # Orientation is π/4 radians -> 45 degrees up
    fr = 1.35*cos(π/4)
    lower = (7 + fr*(6-7))
    upper = (6 + fr*(5-6))
    v1 = lower + fr*(upper-lower)

    lower = (7 + fr*(5-7))
    upper = (5 + fr*(3-5))
    v2 = lower + fr*(upper-lower)

    c = 7.0         # peak value

    # solve v = a*r^2 + b*r + c
    a = (v1 + v2)/2 - c
    b = a + c - v2
    r = -b/2a

    @fact (r - 1/6) --> less_than(EPS)

    transposed = spatialorder(img)[1] == "x"

    # Location and value at peak

    x_peak_offset, y_peak_offset = r*fr, -r*fr
    peakval = a*r^2 + b*r + c

    @fact all(diag(data(t))[2:4] .== peakval) --> true  # Edge pixels aren't interpolated here
    @fact all(t - diagm(diag(data(t))) .== 0) --> true

    diag_s = copyproperties(s, diagm(diag(data(s))))
    @fact s --> diag_s

    @fact all([pt.x for pt in diag(data(s))[2:4]] - ((2:4) + x_peak_offset) .< EPS) --> true
    @fact all([pt.y for pt in diag(data(s))[2:4]] - ((2:4) + y_peak_offset) .< EPS) --> true

end

    context("Nonmax suppression diagonal") do
        # Test image: diagonal edge
        m = [7.0  6.0  5.0  0.0  0.0
             5.0  7.0  6.0  5.0  0.0
             3.0  5.0  7.0  6.0  5.0
             0.0  3.0  5.0  7.0  6.0
             0.0  0.0  3.0  5.0  7.0]

        m_xy = grayim(m')
        m_yx = grayim(m)
        m_yx["spatialorder"] = ["y","x"]

        nms_test_diagonal(m)
        nms_test_diagonal(m_xy)
        nms_test_diagonal(m_yx)
    end

end


text:
using Distributions

include("simutil.jl")

struct SAParameter
    tmin
    tmax
    steps
end

mutable struct SAModel
    SAParameter::SAParameter
    move::Function
    energy::Function # energy function
end

mutable struct Trial
    state
    energy::Float32
end

function fit!(SAModel, X, y)
    println("Simulated Annealing needs not to train")
    SAModel
end

function predict(
    SAModel,
    initstate,
    )

    tfactor = - log(SAModel.SAParameter.tmax / SAModel.SAParameter.tmin)
    nowtrial, besttrial = init(initstate, SAModel.energy)
    for step in 1:SAModel.SAParameter.steps
        t = SAModel.SAParameter.tmax * exp(tfactor * step / SAModel.SAParameter.steps)
        nowtrial, prevtrial = anneal_one(nowtrial, SAModel.move, SAModel.energy)
        if check_state_update(nowtrial, prevtrial, t)
            nowtrial, besttrial = update(nowtrial, prevtrial, besttrial)
        else
            nowtrial = prevtrial
        end

        if step % 1000 == 0
            println("step is over ", step)
            println("best energy is ", besttrial.energy)
        end
    end

    println("step is end")
    println("best energy is ", besttrial.energy)
    println("best state is ", besttrial.state)
    return besttrial
end

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
"""
    Code to estimated RV from centroid of core of CCF
Author: Eric Ford
Created: September 2020
"""

"""  Functor to estimate RV based on the centroid of the CCF.
Inputs:
* `vels`: Array of velocites where CCF was evaluated.
* `ccf`:  Array of values of CCF
* `ccf_var`:  Array of variances of CCF
Warning: Not well tested yet
"""
struct MeasureRvFromCCFCentroid <: AbstractMeasureRvFromCCF
    frac_of_width_to_fit::Float64
    measure_width_at_frac_depth::Float64
end

"""
Construct functor to estimate RV based on the CCF.
Optional Arguments:
* `frac_of_width_to_fit`: (0.5)
* `measure_width_at_frac_depth`: (0.5)
"""
function MeasureRvFromCCFCentroid(; frac_of_width_to_fit::Real = default_frac_of_width_to_fit,
                                    measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth )
    @assert 0.25 <= measure_width_at_frac_depth <= 0.75
    @assert 0.1 <= frac_of_width_to_fit <= 2.0
    MeasureRvFromCCFCentroid(frac_of_width_to_fit,measure_width_at_frac_depth)
end

#=
"""
Estimate RV based on centroid velocity of the CCF.
Inputs:
* `vels`: Array of velocites where CCF was evaluated.
* `ccf`:  Array of values of CCF
Only makes sense if the CCF has been variance normalized.  So don't use this version without ccf variances.
"""
function (mrv::MeasureRvFromCCFCentroid)(vels::A1, ccf::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} }
    # find the min and use only that part of the CCF for computing centroid
    amin, inds = find_idx_at_and_around_minimum(vels, ccf,
                        frac_of_width_to_fit=mrv.frac_of_width_to_fit,
                        measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)

    weights = exp.(-0.5.*(ccf[inds].-ccf[amin]))
    v_centroid = sum(weights.*vels[inds]) / sum(weights)
    return (rv=v_centroid, σ_rv=NaN)
end
=#

# See struct for functor due to bug in Julia/Documenter.jl
function (mrv::MeasureRvFromCCFCentroid)(vels::A1, ccf::A2, ccf_var::A3 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1} }
    # find the min and use only that part of the CCF for computing centroid
    amin, inds = find_idx_at_and_around_minimum(vels, ccf,
                        frac_of_width_to_fit=mrv.frac_of_width_to_fit,
                        measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)

    weights = exp.(-0.5.*(ccf[inds].-ccf[amin])) ./ ccf_var[inds]
    v_centroid = sum(weights.*vels[inds]) / sum(weights)
    return (rv=v_centroid, σ_rv=NaN)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockDiscretization(dims)

Discretize spatial region into `dims` blocks.
"""
struct BlockDiscretization{N} <: DiscretizationMethod
  dims::Dims{N}
end

discretize(r::Rectangle, method::BlockDiscretization) =
  RegularGrid(extrema(r)..., dims=method.dims)


text:
#==========================================================================================
                                Constructors
==========================================================================================#
function SymSemiseparableCholesky(U::AbstractArray, W::AbstractArray)
	if size(U,1) == size(W,1) && size(U,2) == size(W,2)
		return SymSemiseparableCholesky(size(U,2),size(U,1),U,W)
	else
		error("Dimension mismatch between the generators U and W")
	end
end
function SymSemiseparableCholesky(K::SymSemiseparableMatrix)
    return SymSemiseparableCholesky(K.n, K.p, K.Ut, ss_create_w(K.Ut, K.Vt))
end
#==========================================================================================
                    Defining Matrix Properties / Overloading Base
==========================================================================================#
Base.Matrix(K::SymSemiseparableCholesky) = getproperty(K,:L)
Base.size(K::SymSemiseparableCholesky) = (K.n, K.n)
function Base.size(K::SymSemiseparableCholesky, d::Int)
    return (1 <= d && d <=2) ? size(K)[d] : throw(ArgumentError("Invalid dimension $d"))
end
function Base.getindex(K::SymSemiseparableCholesky, i::Int, j::Int)
	i >= j && return dot(K.Ut[:,i], K.Wt[:,j])
	return 0
end
function Base.getproperty(K::SymSemiseparableCholesky, d::Symbol)
    if d === :U
        return UpperTriangular(K.Wt'*K.Ut)
    elseif d === :L
        return LowerTriangular(K.Ut'*K.Wt)
    else
        return getfield(K, d)
    end
end
function Base.propertynames(F::SymSemiseparableCholesky, private::Bool=false)
    return (:U, :L, (private ? fieldnames(typeof(F)) : ())...)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")},
		 K::SymSemiseparableCholesky{<:Any,<:AbstractMatrix,<:AbstractMatrix})
    summary(io, K); println(io)
    show(io, mime, K.L)
end

#==========================================================================================
                            Overloading LinearAlgebra routines
==========================================================================================#
function LinearAlgebra.det(L::SymSemiseparableCholesky)
    dd = one(eltype(L))
    @inbounds for i in 1:L.n
        dd *= dot(L.Ut[:,i],L.Wt[:,i])
    end
    return dd
end

function LinearAlgebra.logdet(L::SymSemiseparableCholesky)
    dd = zero(eltype(L))
    @inbounds for i in 1:L.n
        dd += log(dot(L.Ut[:,i],L.Wt[:,i]))
    end
    return dd
end

#==========================================================================================
                        Defining multiplication and inverse
==========================================================================================#
function LinearAlgebra.mul!(y::AbstractArray, L::SymSemiseparableCholesky, b::AbstractArray)
	ss_tri_mul!(y, L.Ut,   L.Wt,   b)
    return y
end
function LinearAlgebra.mul!(y::AbstractArray, L::Adjoint{<:Any,<:SymSemiseparableCholesky}, b::AbstractArray)
	ssa_tri_mul!(y, L.parent.Ut, L.parent.Wt, b)
    return y
end
function LinearAlgebra.ldiv!(y::AbstractArray,F::SymSemiseparableCholesky, b::AbstractArray)
	ss_forward!(y,F.Ut,F.Wt,b)
	return y
end
function LinearAlgebra.ldiv!(y::AbstractArray,F::Adjoint{<:Any,<:SymSemiseparableCholesky}, b::AbstractArray)
	ssa_backward!(y,F.parent.Ut,F.parent.Wt,b)
	return Y
end
function (Base.:\)(F::SymSemiseparableCholesky, b::AbstractVecOrMat)
    y = similar(b)
	ss_forward!(y,F.Ut,F.Wt,b)
	return y
end
function (Base.:\)(F::Adjoint{<:Any,<:SymSemiseparableCholesky}, b::AbstractVecOrMat)
    y = similar(b)
	ssa_backward!(y,F.parent.Ut,F.parent.Wt,b)
	return y
end

#### Inverse of a SymSemiseparableCholesky using ####
function inv(L::SymSemiseparableCholesky, b::AbstractArray)
	return L'\(L\b)
end

#==========================================================================================
                        Defining the computations
==========================================================================================#
"""
    ss_create_w(U,V)

Creates matrix `W` such that the Cholesky factorization of the SymSemiseparable(U,V) matrix
is `L = tril(U*W')` in linear complexity.
"""
function ss_create_w(U,V)
    p,n = size(U)
    wTw = zeros(p,p)
    W = zeros(p,n)
    @inbounds for (u,v,w) in zip(eachcol(U),eachcol(V),eachcol(W))
        w .= (v - wTw*u)/sqrt(abs(dot(u,(v - wTw*u))))
        add_outer_product!(wTw,w)
    end
    return W
end

"""
    ss_tri_mul!(Y,U,W,X)

Computes the multiplization `tril(U*W')*X = Y` in linear complexity.
"""
function ss_tri_mul!(Y,U,W,X)
    p     = size(U,1)
    n_rhs = size(X,2)
    Wbar  = zeros(p,n_rhs)
    @inbounds for (w,u,y,x) in zip(eachcol(W),eachcol(U),eachrow(Y),eachrow(X))
        add_inner_plus!(Wbar,w,x)
        add_Y_tri!(y,Wbar,u)
    end
end
"""
    ssa_tri_mul!(Y,U,W,X)

Computes the multiplization `triu(W*U')*X = Y` in linear complexity.
"""
function ssa_tri_mul!(Y,U,W,X)
    Ubar = U*X
    @inbounds for (w,u,y,x) in zip(eachcol(W),eachcol(U),eachrow(Y),eachrow(X))
        add_Y_tri!(y,Ubar,w)
        add_inner_minus!(Ubar,u,x)
    end
end
"""
    ss_forward!

Solves the linear system `tril(U*W')*X = B`.
"""
function ss_forward!(X,U,W, B)
    p     = size(U,1)
    n_rhs = size(B,2)
    Wbar  = zeros(p, n_rhs)
    @inbounds for (u,w,x,b) in zip(eachcol(U),eachcol(W),eachrow(X),eachrow(B))
        x .= (b - Wbar'*u)/dot(u,w)
        add_inner_plus!(Wbar,w,x)
    end
end
"""
    ssa_backward!(X, U, W, B)

Solves the linear system `triu(W*U')*X = Y`.
"""
function ssa_backward!(X, U, W, B)
    p     = size(U,1)
    n_rhs = size(B,2)
    Ubar  = zeros(p,n_rhs)
    @inbounds for (u,w,x,b) in zip(Iterators.reverse(eachcol(U)),
                                   Iterators.reverse(eachcol(W)),
                                   Iterators.reverse(eachrow(X)),
                                   Iterators.reverse(eachrow(B)))
        x .= (b - Ubar'w)/dot(u,w)
        add_inner_plus!(Ubar,u,x)
    end
end


text:
###Function to create a string that shows the size of the object in memory
function objectsizetostr(a::Any)

    bytes = Base.summarysize(a)
    result = ""
    if bytes < 10*1024
        result = @sprintf("%6d bytes  ",bytes)
    elseif bytes < 1024*1024
        result = @sprintf("%5.3f kB",bytes/(1024))
    else
        result = @sprintf("%5.3f MB",bytes/1024/1024)
    end
    result

end


text:
#!/usr/bin/env julia
using ArgParse

s = ArgParseSettings()

@add_arg_table s begin
    "--geometry"
    "--gamma"
        default = "1.0"
    "--edipole"
    "--initialstate_phi"
    "--initialstate_theta"
    "--opath"
    "--T"
end

parameters = parse_args(s)

const opath = parameters["opath"]
const T = float(eval(Meta.parse(parameters["T"])))
@assert T[1]<T[end]
const γ = float(parameters["gamma"])
@assert 0<γ
const edipole = float(eval(Meta.parse(parameters["edipole"])))

# Initial state (Bloch state)
const phi = float(parameters["initialstate_phi"])
const theta = float(parameters["initialstate_theta"])

const sx0 = cos(phi)*sin(theta)
const sy0 = sin(phi)*sin(theta)
const sz0 = cos(theta)

using QuantumOptics, CollectiveSpins
const geomstring = parameters["geometry"]
const systemgeometry = eval(Meta.parse("CollectiveSpins.geometry.$geomstring"))
const system = SpinCollection(systemgeometry, edipole, γ)
const N = length(system.spins)

const td_ind = Float64[]
const td_mf = Float64[]
const td_mpc = Float64[]

state_master_t = Vector{Float64}[]

# Independent
state0 = CollectiveSpins.independent.blochstate(phi,theta,N)
tout, state_ind_t = CollectiveSpins.independent.timeevolution(T, system, state0)

# Meanfield
state0 = CollectiveSpins.meanfield.blochstate(phi,theta,N)
tout, state_mf_t = CollectiveSpins.meanfield.timeevolution(T, system, state0)

# Meanfield + Correlations
state0 = CollectiveSpins.mpc.blochstate(phi,theta,N)
tout, state_mpc_t = CollectiveSpins.mpc.timeevolution(T, system, state0)

# Quantum: master equation

sx_operators = Operator[quantumoptics.embed(CollectiveSpins.quantum.basis(system), [i], [sigmax]) for i=1:N]
sy_operators = Operator[quantumoptics.embed(CollectiveSpins.quantum.basis(system), [i], [sigmay]) for i=1:N]
sz_operators = Operator[quantumoptics.embed(CollectiveSpins.quantum.basis(system), [i], [sigmaz]) for i=1:N]

function fout(t, rho::Operator)
    i = findfirst(T, t)
    rho_ind = CollectiveSpins.independent.densityoperator(state_ind_t[i])
    rho_mf = CollectiveSpins.meanfield.densityoperator(state_mf_t[i])
    rho_mpc = CollectiveSpins.mpc.densityoperator(state_mpc_t[i])
    push!(td_ind, tracedistance(rho, rho_ind))
    push!(td_mf, tracedistance(rho, rho_mf))
    push!(td_mpc, tracedistance(rho, rho_mpc))
    state = zeros(Float64, 3*N)
    state[0*N+1:1*N] = [abs(expect(sx_operators[i], rho)) for i=1:N]
    state[1*N+1:2*N] = [abs(expect(sy_operators[i], rho)) for i=1:N]
    state[2*N+1:3*N] = [abs(expect(sz_operators[i], rho)) for i=1:N]
    push!(state_master_t, state)
end

Ψ₀ = CollectiveSpins.quantum.blochstate(phi,theta,N)
ρ₀ = Ψ₀⊗dagger(Ψ₀)
CollectiveSpins.quantum.timeevolution(T, system, ρ₀, fout=fout)

f = open(opath, "w")
CollectiveSpins.io.write_head(f, system, parameters)
CollectiveSpins.io.write_state(f, "timevector", T; time=false)
CollectiveSpins.io.write_state(f, "td_ind", td_ind; time=false)
CollectiveSpins.io.write_state(f, "td_mf", td_mf; time=false)
CollectiveSpins.io.write_state(f, "td_mpc", td_mpc; time=false)

write(f, "<master>\n")
for i=1:length(T)
    CollectiveSpins.io.write_state(f, "meanfield", state_master_t[i]; time=false)
end
write(f, "</master>\n")
write(f, "<meanfield>\n")
for i=1:length(T)
    CollectiveSpins.io.write_state(f, "meanfield", state_mf_t[i]; time=false)
end
write(f, "</meanfield>\n")
write(f, "<mpc>\n")
for i=1:length(T)
    sx = CollectiveSpins.mpc.sx(state_mpc_t[i])
    sy = CollectiveSpins.mpc.sy(state_mpc_t[i])
    sz = CollectiveSpins.mpc.sz(state_mpc_t[i])
    state = zeros(Float64, 3*N)
    for i=1:N
        state[i] = sx[i]
        state[N+i] = sy[i]
        state[2*N+i] = sz[i]
    end
    CollectiveSpins.io.write_state(f, "meanfield", state; time=false)
end
write(f, "</mpc>\n")

close(f)

# Visualization
# using PyCall
# @pyimport matplotlib.pyplot as plt

# plt.figure(figsize=(6,4))
# plt.plot(T, td_ind, lw=3.0, color="darkorange", label="Independent")
# plt.plot(T, td_mf, lw=3.0, color="gray", label="meanfield")
# plt.plot(T, td_mpc, "-", lw=1.5, color="navy", label="MPC")
# plt.ylabel("\$-\\langle\\sigma_z\\rangle\$")
# plt.ylim(0,1)
# plt.legend()

# plt.show()


text:
export DeepQuadratic_fast,DeepQuadratic

"""
    DeepQuadratic_fast(param;label,dealias)

Same as `DeepQuadratic`, but faster.
"""
mutable struct DeepQuadratic_fast <: AbstractModel

    label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	param	:: NamedTuple
	kwargs	:: NamedTuple

    function DeepQuadratic_fast( param::NamedTuple; label="Deep quadratic", dealias=true)

		@warn("You should provide the initial value for v by ∂t η = - ∂x v.
		It equals the tangential velocity only when they are null.")

        ϵ = param.ϵ
        mesh  = Mesh(param)
		param = ( ϵ = ϵ, μ=Inf, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N )
		kwargs = (label=label, )

        x = mesh.x
        Γ = abs.(mesh.k)
        Dx    =  1im * mesh.k            # Differentiation
        H     = -1im * sign.(mesh.k)     # Hilbert transform
		if dealias == true || dealias == 1
			Π⅔    = Γ .< mesh.kmax * 2/3     # Dealiasing low-pass filter
		else
			Π⅔    = zero(Γ) .+ 1     # Dealiasing low-pass filter
		end

        hnew = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

        Px  = plan_fft(hnew; flags = FFTW.MEASURE)

		function f!(U)

		    ldiv!(hnew, Px , view(U,:,1))

		    I₁  .= view(U,:,2) .* Γ
		    ldiv!(unew, Px , I₁)
		    unew .^= 2
		    mul!(I₁, Px , unew)
		    I₁ .*= H

		    I₂  .= view(U,:,1) .* Dx
		    ldiv!(unew, Px , I₂)
		    unew .*= hnew
		    mul!(I₂, Px , unew)

		    I₃  .= view(U,:,1) .* Γ
		    ldiv!(unew, Px, I₃)
		    unew .*= hnew
		    mul!(I₃ , Px , unew)
		    I₃ .*= H

		    hnew  .= .- view(U,:,2) .* Dx

		    I₁ .-= I₂
		    I₁ .-= I₃
		    I₁ .*= Π⅔
		    I₁ .*= ϵ

		    U[:,2] .= view(U,:,1) .* H .+ I₁
		    U[:,1] .= hnew

		end

		# This function is correct only when the initial data for v is zero.
		function mapto(data::InitialData)
		    [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is detemined by ∂t η + ∂x v = 0.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2]))
		end

		new(label, f!, mapto, mapfro, param, kwargs )

    end
end

"""
    DeepQuadratic(param;label,dealias)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the quadratic deep-water model proposed by [Akers and Milewski](https://doi.org/10.1137/090758386)
and [Cheng, Granero-Belinchón, Shkoller and Milewski](https://dx.doi.org/10.1007/s42286-019-00005-w)

# Arguments
`param` is of type `NamedTuple` and must contain
- the dimensionless parameters `ϵ` (nonlinearity);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`

`label` is optional and provides a label for future plottings ("Deep quadratic" is used if not provided).

`dealias`: dealiasing with `1/3` Orlicz rule if `true` (by default) or no dealiasing if `false`.

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!` and in particular
1. a function `DeepQuadratic.f!` to be called in the time-integration solver;
2. a function `DeepQuadratic.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `DeepQuadratic.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where
    - `η` is the surface deformation;
    - `v` is given by `∂t η = - ∂x v`.

"""
mutable struct DeepQuadratic <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	param	:: NamedTuple
	kwargs	:: NamedTuple

    function DeepQuadratic( param::NamedTuple; label="Deep quadratic", dealias=true)

		@warn("You should provide the initial value for v by ∂t η = - ∂x v.
		It equals the tangential velocity only when they are null.")

		ϵ = param.ϵ
		mesh  = Mesh(param)
		param = ( ϵ = ϵ, μ=Inf, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N )
		kwargs = (label=label, )

		x = mesh.x
        Γ = abs.(mesh.k)
        Dx    =  1im * mesh.k           # Differentiation
        H     = -1im * sign.(mesh.k)    # Hilbert transform
		if dealias == true || dealias == 1
			Π⅔    = Γ .< mesh.kmax * 2/3     # Dealiasing low-pass filter
		else
			Π⅔    = zero(Γ) .+ 1     # Dealiasing low-pass filter
		end

        hnew = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

		function f!(U)

			hnew .=ifft(U[:,1]);
			I₁ .=H.*fft(ifft(Γ.*U[:,2]).^2);
			I₂ .=fft(hnew.*ifft(Dx.*U[:,1])) ;
			I₃ .=H.*fft(hnew.*ifft(Γ.*U[:,1]));
			hnew .= U[:,1] ;
			U[:,1] .= -(Dx.*U[:,2]) ;
			U[:,2] .= H.*hnew+ϵ*Π⅔.*(I₁-I₂-I₃) ;

		end

		# This function is correct only when the initial data for v is zero.
		function mapto(data::InitialData)
			[Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is detemined by ∂t η + ∂x v = 0.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2]))
		end

		new(label, f!, mapto, mapfro, param, kwargs )

    end
end


text:
TABLE_CAST = """id,height,age,name,occupation
1.1,10,1,string1,2012
"""

@testset "Read a Table from file" begin

    @testset "Basic data reading" begin
        t = Table("../data/data_types.csv")
        trs = TableSchema.read(t, cast=false)
        # check the headers
        @test length(t.headers) == 5
        HEIGHT_COLUMN = 2
        @test t.headers[HEIGHT_COLUMN] == "height"
        # check the number of rows
        @test length(trs[:,1]) == 5
        # check the bottom left index
        @test trs[5,1] == 5
        # iterate over the rows
        @test sum([ row[HEIGHT_COLUMN] for row in t ]) == 51
        # no schema, hence exception
        @test_throws TableValidationException validate(t)
    end

    @testset "Passing an IO buffered file" begin
        t = Table()
        t.schema = Schema("../data/schema_valid_missing.json")
        @test TableSchema.is_valid(t.schema)
        f = readdlm("../data/data_types.csv", ',')
        tr = TableSchema.read(t, data=f)
        @test validate(t)
    end

    @testset "With schema validation" begin
        s = Schema("../data/schema_valid_missing.json")
        t = Table("../data/data_types.csv", s)
        @test validate(t)
        t2 = Table("../data/data_constraints.csv", s)
        @test !validate(t2)
        @test length(t2.errors) == 2
    end

    @testset "With row casting" begin
        s = Schema("../data/schema_valid_missing.json")
        t = Table(IOBuffer(TABLE_CAST), s)
        tr = TableSchema.read(t)
        @test TableSchema.is_valid(t.schema)
        @test isa(tr[1,1], String)
        @test isa(tr[1,2], Float64)
    end

    @testset "Remote data and schema reading" begin
        t = Table("https://raw.githubusercontent.com/frictionlessdata/tableschema-jl/master/data/data_types.csv")
        trs = TableSchema.read(t, cast=false)
        # check the headers
        @test length(t.headers) == 5
        @test t.headers[2] == "height"
        # validate from remote schema
        t.schema = Schema("https://raw.githubusercontent.com/frictionlessdata/tableschema-jl/master/data/schema_valid_missing.json")
        @test validate(t)
    end

end


text:
module MatrixMath

using StaticArrays
using Quaternions

export array
export avec

# converts to normal array
array(x::AbstractArray) = convert(array(typeof(x)),vec(x))
array{T<:AbstractArray}(::Type{T}) = Array{array(eltype(T)),1};
array{T}(::Type{T}) = T;

# converts to one-dim normal array (faster than using vcat(x...) only! you need to convert into normal array first!)
avec(x::AbstractArray) = vcat(array(x)...)

export Vec
export Mat

const Vec = SVector
const Mat = SMatrix

export Vec2f
export Vec3f
export Vec4f

const Vec2f = Vec{2, Float32}
const Vec3f = Vec{3, Float32}
const Vec4f = Vec{4, Float32}

export Vec2i
export Vec3i
export Vec4i

const Vec2i = Vec{2, Int32}
const Vec3i = Vec{3, Int32}
const Vec4i = Vec{4, Int32}

export Vec2u
export Vec3u
export Vec4u

const Vec2u = Vec{2, UInt32}
const Vec3u = Vec{3, UInt32}
const Vec4u = Vec{4, UInt32}

export Mat2x2f
export Mat3x3f
export Mat4x4f

const Mat2x2f = Mat{2,2,Float32,4} #StaticArrays.SArray{Tuple{2,2},Float32,2,4}
const Mat3x3f = Mat{3,3,Float32,8} ##StaticArrays.SArray{Tuple{3,3},Float32,2,8}
const Mat4x4f = Mat{4,4,Float32,16} ##StaticArrays.SArray{Tuple{4,4},Float32,2,16}

export IMat2x2f
export IMat3x3f
export IMat4x4f

const IMat2x2f = @SMatrix [
	1f0 0 ;
	0 1f0 ;
]

const IMat3x3f = @SMatrix [
	1f0 0 0 ;
	0 1f0 0 ;
	0 0 1f0 ;
]

const IMat4x4f = @SMatrix [
	1f0 0 0 0 ;
	0 1f0 0 0 ;
	0 0 1f0 0 ;
	0 0 0 1f0 ;
]

export AVec2f
export AVec3f
export AVec4f

const AVec2f = Array{Vec2f,1}
const AVec3f = Array{Vec3f,1}
const AVec4f = Array{Vec4f,1}

#t = MatrixMath.translationmatrix(Vec3f(0,0,-5))
#c = MatrixMath.lookat(Vec3f(sin(-(cursorpos[1])*6.3),0,1-cos((cursorpos[1])*6.3)),Vec3f(0,0,1), Vec3f(0,1,0))
#c = MatrixMath.rotate(-sin(-0.5f0+cursorpos[1]),Vec3f(0,1,0))

export rotr90flipdim2
rotr90flipdim2(x::AbstractArray) = flipdim(rotr90(x),2)

function convertMatrixToArray(values::AbstractArray)
	dims=ndims(values)
	elems=dims > 1 ? size(values)[dims] : 1
	values=avec(values) #wow! < 0.05 sec for > 1000 vertices!
	#print("vcat: "); @time v=vcat(values...) #end # convert to one dimensional array
	#print("convert: "); @time c=convert(Array{typeof(v[1]),1},v)
	(values, elems)
end

#c = FPSViewRH(normalize(Vec3f(0,0,1)),(cursorpos[1]-0.5f0)*3,(cursorpos[2]-0.5f0) * 3)
function FPSViewRH(eye::Any, yaw::Float32, pitch::Float32)
    # If the pitch and yaw angles are in degrees,
    # they need to be converted to radians. Here
    # I assume the values are already converted to radians.
    cosPitch = cos(pitch)
    sinPitch = sin(pitch)
    cosYaw = cos(yaw)
    sinYaw = sin(yaw)

    xaxis = Vec3f( cosYaw, 0, -sinYaw )
    yaxis = Vec3f( sinYaw * sinPitch, cosPitch, cosYaw * sinPitch )
    zaxis = Vec3f( sinYaw * cosPitch, -sinPitch, cosPitch * cosYaw )

    # Create a 4x4 view matrix from the right, up, forward and eye position vectors
    Mat{4,4,Float32}([
      xaxis[1] yaxis[1] zaxis[1] -dot(xaxis,eye)*0;
      xaxis[2] yaxis[2] zaxis[2] -dot(yaxis,eye)*0;
      xaxis[3] yaxis[3] zaxis[3] -dot(zaxis,eye)*0;
      -dot(xaxis,eye) -dot(yaxis,eye) -dot(zaxis,eye) 1;
    ])
end

function perspectiveprojection{T}(fovy::T, aspect::T, znear::T, zfar::T)
		(znear == zfar) && error("znear ($znear) must be different from tfar ($zfar)")
		zfn = 1/(zfar-znear)

		t = tan(fovy * 0.5)
    h = T(tan(fovy * pi / 360) * znear)
    w = T(h * aspect)

		left = -w
		right = w
		bottom = -h
		top = h

		frustum(-w, w, -h, h, znear, zfar)

    #Mat{4,4,T}([
    #    1/(aspect*t) 0 0 0;
    #    0 1/t 0 0;
    #    0 0 -(zfar+znear)*zfn -(2*zfar*znear)*zfn;
		#		0 0 -1 0;
    # ])
end

function frustum{T}(left::T, right::T, bottom::T, top::T, znear::T, zfar::T)
    (right == left || bottom == top || znear == zfar) && return eye(Mat{4,4,T})
    Mat{4,4,T}([
        2*znear/(right-left) 0 0 0;
        0 2*znear/(top-bottom) 0 0;
        (right+left)/(right-left) (top+bottom)/(top-bottom) (-(zfar+znear)/(zfar-znear)) -(2*znear*zfar) / (zfar-znear);
        0 0 -1 0
    ])
end

function translationmatrix{T}(t::Vec{3, T})
    Mat{4,4,T}([
        1 0 0 t[1];
        0 1 0 t[2];
        0 0 1 t[3];
				0 0 0 1;
    ])
end

function inverse_translationmatrix{T}(t::Vec{3, T})
    Mat{4,4,T}([
        1 0 0 0;
        0 1 0 0;
        0 0 1 0;
				t[1] t[2] t[3] 1;
    ])
end

function rotationmatrix{T}(t::Vec{3, T})
    Mat{4,4,T}([
        1 0 0 0;
        0 1 0 0;
        0 0 1 0;
				t[1] t[2] t[3] 1;
    ])
end

function scalingmatrix{T}(t::Vec{3, T})
    Mat{4,4,T}([
        t[1] 0 0 0;
        0 t[2] 0 0;
        0 0 t[3] 0;
				0 0 0 1;
    ])
end

function orthographicprojection{T}(left::T,right::T,bottom::T,top::T,znear::T,zfar::T)
    (right==left || bottom==top || znear==zfar) && return eye(Mat{4,4,T})
    Mat{4,4,T}([
        2/(right-left) 0 0 0;
        0 2/(top-bottom) 0 0;
        0 0 -2/(zfar-znear) 0;
				-(right+left)/(right-left) -(top+bottom)/(top-bottom) -(zfar+znear)/(zfar-znear) 1;
    ])
end

function lookat{T}(eye::Vec{3, T}, lookAt::Vec{3, T}, up::Vec{3, T})
    zaxis  = normalize(eye-lookAt)
    xaxis  = normalize(cross(up,    zaxis))
    yaxis  = normalize(cross(zaxis, xaxis))
    Mat{4,4,T}([
        xaxis[1] yaxis[1] zaxis[1] 0;
        xaxis[2] yaxis[2] zaxis[2] 0;
        xaxis[3] yaxis[3] zaxis[3] 0;
        -dot(xaxis,eye) -dot(yaxis,eye) -dot(zaxis,eye) 1
    ])
end

rotate{T}(angle::T, axis::Vec{3, T}) = rotationmatrix4(Quaternions.qrotation(convert(Array, axis), angle))
rotate{T}(v::Vec{2, T}, angle::T) = Vec{2,T}(v[1] * cos(angle) - v[2] * sin(angle), v[1] * sin(angle) + v[1] * cos(angle))

function rotationmatrix4{T}(q::Quaternions.Quaternion{T})
    sx, sy, sz = 2q.s*q.v1,  2q.s*q.v2,   2q.s*q.v3
    xx, xy, xz = 2q.v1^2,    2q.v1*q.v2,  2q.v1*q.v3
    yy, yz, zz = 2q.v2^2,    2q.v2*q.v3,  2q.v3^2
    Mat{4,4,T}([
        1-(yy+zz)	xy+sz				xz-sy				0;
        xy-sz				1-(xx+zz)	yz+sx				0;
        xz+sy				yz-sx				1-(xx+yy)	0;
        0 0 0 1
    ])
end

forwardVector4{T}(m::Mat{4,4,T}) = Vec{3, T}(m[3,1],m[3,2],m[3,3])
rightVector4{T}(m::Mat{4,4,T}) = Vec{3, T}(m[1,1],m[1,2],m[1,3])
upVector4{T}(m::Mat{4,4,T}) = Vec{3, T}(m[2,1],m[2,2],m[2,3])

function computeRotation{T}(rot::Vec{3,T})
	dirBackwards= Vec{3, T}(-1,0,0)
	dirRight = Vec{3, T}(0,0,1)
	dirUp = Vec{3, T}(0,1,0) #cross(dirRight, dirBackwards)

	q = Quaternions.qrotation(convert(Array, dirRight), rot[3]) *
	Quaternions.qrotation(convert(Array, dirUp), rot[1]) *
	Quaternions.qrotation(convert(Array, dirBackwards), rot[2])

	rotationmatrix4(q)
end

#=
using GeometryTypes

function toOneArray{S1,S2,T}(m::Mat{S1,S2,T})
	a = T[]
	l(x,y) = (push!(a, x); y)
	push!(a, reduce(l, convert(Array{T,}, m)))
	a
end

function scalematrix{T}(s::Vec{3, T})
    T0, T1 = zero(T), one(T)
    Mat{4,4,T}(
        (s[1],T0,  T0,  T0),
        (T0,  s[2],T0,  T0),
        (T0,  T0,  s[3],T0),
        (T0,  T0,  T0,  T1),
    )
end

translationmatrix_x{T}(x::T) = translationmatrix(Vec{3, T}(x, 0, 0))
translationmatrix_y{T}(y::T) = translationmatrix(Vec{3, T}(0, y, 0))
translationmatrix_z{T}(z::T) = translationmatrix(Vec{3, T}(0, 0, z))

function translationmatrix{T}(t::Vec{3, T})
    T0, T1 = zero(T), one(T)
    Mat(
        (T1,  T0,  T0,  T0),
        (T0,  T1,  T0,  T0),
        (T0,  T0,  T1,  T0),
        (t[1],t[2],t[3],T1),
    )
end

rotate{T}(angle::T, axis::Vec{3, T}) = rotationmatrix4(Quaternions.qrotation(convert(Array, axis), angle))

function rotationmatrix_x{T}(angle::T)
    T0, T1 = zero(T), one(T)
    Mat{4,4,T}(
        (T1, T0, T0, T0),
        (T0, cos(angle), sin(angle), T0),
        (T0, -sin(angle), cos(angle),  T0),
        (T0, T0, T0, T1)
    )
end
function rotationmatrix_y{T}(angle::T)
    T0, T1 = zero(T), one(T)
    Mat{4,4,T}(
        (cos(angle), T0, -sin(angle),  T0),
        (T0, T1, T0, T0),
        (sin(angle), T0, cos(angle), T0),
        (T0, T0, T0, T1)
    )
end
function rotationmatrix_z{T}(angle::T)
    T0, T1 = zero(T), one(T)
    Mat{4,4,T}(
        (cos(angle), sin(angle), T0, T0),
        (-sin(angle), cos(angle),  T0, T0),
        (T0, T0, T1, T0),
        (T0, T0, T0, T1)
    )
end
#=
    Create view frustum
    Parameters
    ----------
        left : float
         Left coordinate of the field of view.
        right : float
         Left coordinate of the field of view.
        bottom : float
         Bottom coordinate of the field of view.
        top : float
         Top coordinate of the field of view.
        znear : float
         Near coordinate of the field of view.
        zfar : float
         Far coordinate of the field of view.
    Returns
    -------
        M : array
         View frustum matrix (4x4).
=#
function frustum{T}(left::T, right::T, bottom::T, top::T, znear::T, zfar::T)
    (right == left || bottom == top || znear == zfar) && return eye(Mat{4,4,T})
    T0, T1, T2 = zero(T), one(T), T(2)
    Mat{4,4,T}(
        (T2 * znear / (right - left), T0, T0, T0),
        (T0, T2 * znear / (top - bottom), T0, T0),
        ((right + left) / (right - left), (top + bottom) / (top - bottom), -(zfar + znear) / (zfar - znear), -T1),
        (T0, T0, (-T2 * znear * zfar) / (zfar - znear), T0)
    )
end

perspectiveprojection{T}(wh::SimpleRectangle, fov::T, near::T, far::T) = perspectiveprojection(fov, T(wh.w/wh.h), near, far)
function perspectiveprojection{T}(fovy::T, aspect::T, znear::T, zfar::T)
    (znear == zfar) && error("znear ($znear) must be different from tfar ($zfar)")
    h = T(tan(fovy / 360.0 * pi) * znear)
    w = T(h * aspect)
    frustum(-w, w, -h, h, znear, zfar)
end

orthographicprojection{T}(wh::SimpleRectangle, near::T, far::T) =
    orthographicprojection(zero(T), T(wh.w), zero(T), T(wh.h), near, far)

function orthographicprojection{T}(
        left  ::T, right::T,
        bottom::T, top  ::T,
        znear ::T, zfar ::T
    )
    (right==left || bottom==top || znear==zfar) && return eye(Mat{4,4,T})
    T0, T1, T2 = zero(T), one(T), T(2)
    Mat{4,4,T}(
        (T2/(right-left), T0, T0,  T0),
        (T0, T2/(top-bottom), T0,  T0),
        (T0, T0, -T2/(zfar-znear), T0),
        (-(right+left)/(right-left), -(top+bottom)/(top-bottom), -(zfar+znear)/(zfar-znear), T1)
    )
end
=#
end #MatrixMath

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
# # Generate upper left coordinates for specic index
# function coords(ga::GeoArray, p::SVector{2, Int})
#     ga.f(p.-1)
# end
# coords(ga::GeoArray, p::Vector{Int}) = coords(ga, SVector{2}(p))

# # Generate coordinates for complete GeoArray
# function coords(ga::GeoArray)
#     (ui, uj) = size(ga)[1:2]
#     ci = [coords(ga, SVector{2}(i,j)) for i in 1:ui+1, j in 1:uj+1]
# end
# function coords2(ga::GeoArray, p::SVector{2, Int}, mid::Vector{Int} = [1, 1])
#     ga.f(p .- 1 + mid./2)
# end
# coords2(ga::GeoArray, p::Vector{Int}, mid::Vector{Int} = [1, 1]) = coords2(ga, SVector{2}(p), mid)

# function coords2(ga::GeoArray, mid::Vector{Int} = [1, 1])
#     (ui, uj) = size(ga)[1:2]
#     ci = [coords2(ga, SVector{2}(i,j), mid) for i in 1:ui, j in 1:uj]
# end

# # Generate coordinates for one dimension of a GeoArray
# function coords(ga::GeoArray, dim::Symbol)
#     if is_rotated(ga)
#         error("This method cannot be used for a rotated GeoArray")
#     end
#     if dim==:x
#         ui = size(ga,1)
#         ci = [coords(ga, SVector{2}(i,1))[1] for i in 1:ui+1]
#     elseif dim==:y
#         uj = size(ga,2)
#         ci = [coords(ga, SVector{2}(1,j))[2] for j in 1:uj+1]
#     else
#         error("Use :x or :y as second argument")
#     end
#     return ci
# end

function meshgrid(x::AbstractArray{T,1}, y::AbstractArray{T,1}) where T <: Real
    X = x .* ones(1, length(y))
    Y = ones(length(x)) .* y'
    X, Y
end

"""
    st_dim(ga::GeoArray; mid::Vector{Int} = [1, 1])
    st_dim(b::bbox, cellsize::T) where {T <: Real}
    st_dim(ga::GeoArray, dim::Symbol; mid::Vector{Int} = [1, 1])

Get dimensions of x and y
"""
function st_dim(ga::GeoArray; mid::Vector{Int} = [1, 1])
    if length(mid) == 1; mid = [mid, mid]; end

    cellsize_x = ga.f.linear[1]
    cellsize_y = abs(ga.f.linear[4])
    cellsize_y2 = ga.f.linear[4]

    delta = [cellsize_x, cellsize_y]/2 .* mid

    rbbox = st_bbox(ga)
    x = rbbox.xmin + delta[1]:cellsize_x:rbbox.xmax
    y = rbbox.ymin + delta[2]:cellsize_y:rbbox.ymax
    if cellsize_y2 < 0; y = reverse(y); end

    x, y
end

function st_dim(b::bbox, cellsize::T) where {T <: Real}
    lon = b.xmin + cellsize/2 : cellsize : b.xmax
    lat = reverse(b.ymin + cellsize/2 : cellsize : b.ymax)
    lon, lat # return
end

st_dim_x(ga::GeoArray; mid::Vector{Int} = [1, 1]) = st_dim(ga; mid)[1]
st_dim_y(ga::GeoArray; mid::Vector{Int} = [1, 1]) = st_dim(ga; mid)[2]

function st_dim(ga::GeoArray, dim::Symbol; mid::Vector{Int} = [1, 1])
    if dim==:x
        ci = st_dim_x(ga; mid = mid)
    elseif dim==:y
        ci = st_dim_y(ga; mid = mid)
    else
        error("Use :x or :y as second argument")
    end
    return ci
end

export meshgrid, st_dim