text:
#=--------------------------------------------------------------------

transforms3d - Functions for performing 3D transformations

Copyright (c) 2015-2017 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK August    2015
   April     2017  Cleanups and optimisations
   September 2018  Updated for v0.7/v1.0

---------------------------------------------------------------------=#

export trans, rotx, roty, rotz, invht
export dhtrans, homotrans
export inveuler, invrpy
export angleaxis, angleaxis2matrix, angleaxisrotate, normaliseangleaxis
export matrix2angleaxis, matrix2quaternion, quaternion2matrix
export quaternion, quaternionconjugate, quaternionproduct, quaternionrotate
export vector2quaternion

using LinearAlgebra

#---------------------------------------------------------------------
"""
trans - Homogeneous transformation for a translation by x, y, z
```
Usage: T = trans(x, y, z)
       T = trans(v)

Arguments:  x,y,z - translations in x,y and z, or alternatively
            v     - 3x1 vector or array defining x, y and z.
Returns:    T     - 4x4 homogeneous transformation matrix.
```
See also: rotx(), roty(), rotz(), invht()
"""
function trans(x::Real, y::Real, z::Real)

  T = [ 1.0  0.0  0.0   float(x)
        0.0  1.0  0.0   float(y)
        0.0  0.0  1.0   float(z)
        0.0  0.0  0.0  1.0 ]
end

function trans(t::Array{T1}) where T1<:Real

    if length(t) != 3
        error("Translation must be a 3x1 vector or array")
    end

    T = [ 1.0  0.0  0.0  float(t[1])
          0.0  1.0  0.0  float(t[2])
          0.0  0.0  1.0  float(t[3])
          0.0  0.0  0.0  1.0 ]
end

#----------------------------------------------------------------------
"""
rotx - Homogeneous transformation for a rotation about the x axis
```
Usage: T = rotx(theta)

Argument:  theta  - rotation about x axis
Returns:    T     - 4x4 homogeneous transformation matrix
```
See also: trans(), roty(), rotz(), invht()
"""
function rotx(theta::Real)
    T = [ 1.0     0.0         0.0      0.0
          0.0  cos(theta) -sin(theta)  0.0
          0.0  sin(theta)  cos(theta)  0.0
          0.0     0.0         0.0      1.0 ]
end

#---------------------------------------------------------------------
"""
roty - Homogeneous transformation for a rotation about the y axis
```
Usage: T = roty(theta)

Argument:  theta  - rotation about y axis
Returns:    T     - 4x4 homogeneous transformation matrix
```
See also: trans(), rotx(), rotz(), invht()
"""
function roty(theta::Real)
    T = [ cos(theta)  0.0  sin(theta)  0.0
              0.0     1.0      0.0     0.0
         -sin(theta)  0.0  cos(theta)  0.0
              0.0     0.0      0.0     1.0 ]
end

#---------------------------------------------------------------------
"""
rotz - Homogeneous transformation for a rotation about the z axis
```
Usage: T = rotz(theta)

Argument:  theta  - rotation about z axis
Returns:    T     - 4x4 homogeneous transformation matrix
```
See also: trans(), rotx(), roty(), invht()
"""
function rotz(theta::Real)
    T = [ cos(theta) -sin(theta)  0.0   0.0
          sin(theta)  cos(theta)  0.0   0.0
              0.0        0.0      1.0   0.0
              0.0        0.0      0.0   1.0 ]
end

#---------------------------------------------------------------------
"""
invht - Inverse of a homogeneous transformation matrix
```
Usage:  Tinv = invht(T)

Argument:   T    - 4x4 homogeneous transformation matrix
Returns:    Tinv - inverse
```
See also: trans(), rotx(), roty(), rotz()
"""
function invht(T::Array{T1,2}) where T1 <: Real

    if size(T) != (4,4)
        error("T must be 4 x 4")
    end

    A = T[1:3,1:3]

    Tinv = [      A'        -A'*T[1:3,4]
            0.0  0.0  0.0      1.0      ]
end

#---------------------------------------------------------------------
"""
angleaxis2matrix - Converts angle-axis descriptor to 4x4 homogeneous
transformation  matrix
```
Usage:     T = angleaxis2matrix(t)

Argument:  t - 3x1 vector specifying the rotation axis and having
               magnitude equal to the rotation angle in radians.
Returns:   T - 4x4 Homogeneous transformation matrix.
```
See also: matrix2angleaxis(), angleaxisrotate(), angleaxis(), normaliseangleaxis()
"""
function  angleaxis2matrix(t::Vector{T1}) where T1 <: Real

    if length(t) != 3
        error("Angle-axis must be a 3-vector")
    end

    theta = norm(t)
    if theta < eps()    # If the rotation is very small...
        T = [  1.0  -t[3]   t[2]   0.0
              t[3]   1.0   -t[1]   0.0
             -t[2]   t[1]    1.0   0.0
               0.0   0.0     0.0   1.0]

        return T
    end

    # Otherwise set up standard matrix, first setting up some convenience
    # variables
    tn = t./theta;  x = tn[1]; y = tn[2]; z = tn[3];

    c = cos(theta); s = sin(theta); C = 1.0-c
    xs = x*s;   ys = y*s;   zs = z*s
    xC = x*C;   yC = y*C;   zC = z*C
    xyC = x*yC; yzC = y*zC; zxC = z*xC

    T = [ x*xC+c   xyC-zs   zxC+ys  0.0
          xyC+zs   y*yC+c   yzC-xs  0.0
          zxC-ys   yzC+xs   z*zC+c  0.0
           0.0       0.0     0.0    1.0]
end

#---------------------------------------------------------------------
"""
angleaxisrotate - Uses angle axis descriptor to rotate vectors
```
Usage: v2 = angleaxisrotate(t, v)

Arguments:  t  - 3-vector defining rotation axis and having magnitude
                 equal to the rotation angle in radians.
            v  - 4xn matrix of homogeneous 4-vectors to be rotated or
                 3xn matrix of inhomogeneous 3-vectors to be rotated.
Returns:    v2 - The rotated vectors.
```
See also: matrix2angleaxis(), angleaxis(), angleaxis2matrix(), normaliseangleaxis()
"""
function angleaxisrotate(t::Vector, v::Array)

    ndim = size(v,1)
    T = angleaxis2matrix(t)

    if ndim == 3
      v2 = T[1:3,1:3]*v

    elseif ndim == 4
      v2 = T*v

    else
      error("v must be 4xN or 3xN");
    end

    return v2
end

#---------------------------------------------------------------------
"""
dhtrans - Computes Denavit Hartenberg matrix

This function calculates the 4x4 homogeneous transformation matrix,
representing the Denavit Hartenberg matrix of a robot arm link, given
link parameters of joint angle, length, joint offset and twist.
```
Usage: T = dhtrans(theta, offset, length, twist)

Arguments:  theta - joint angle (rotation about local z)
           offset - offset (translation along z)
           length - translation along link x axis
            twist - rotation about link x axis

Returns:        T - 4x4 Homogeneous transformation matrix.
```
See also: trans(), rotx(), roty(), rotz(), invht()
"""
function dhtrans(theta::Real, offset::Real, length::Real, twist::Real)

T = [ cos(theta) -sin(theta)*cos(twist)  sin(theta)*sin(twist) length*cos(theta)
      sin(theta)  cos(theta)*cos(twist) -cos(theta)*sin(twist) length*sin(theta)
          0.0           sin(twist)             cos(twist)         offset
          0.0              0.0                    0.0               1.0          ]
end

#---------------------------------------------------------------------
"""
homotrans - Homogeneous transformation of points/lines

Function to perform a transformation on 2D or 3D homogeneous or
inhomogeneous coordinates. The resulting coordinates are normalised to
have a homogeneous scale of 1.

```
Usage:     t = homotrans(P, v; checkinf = false)

Arguments:
          P  - 3 x 3 or 4 x 4 homogeneous transformation matrix
          v  - 3 x n or 4 x n matrix of homogeneous coordinates or ...
               2 x n or 3 x n matrix of inhomogeneous coordinates
Keyword Argument:
    checkinf - If true code will check if the homogeneous coordinates
               of any transformed point is at infinity (scale parameter
               is zero). If so, no normalisation of the point is attempted.

Returns   t  - Transformed coordinates. Homogeneous or inhomogeneous
               to match input data v.
```

The input v is assumed inhomogeneous if the number of rows is equal to
size(P,2) - 1

"""
function homotrans(P::Array, v::Array; checkinf = false)
    # ? Should normalisation of scale test for points at infinity ?

    (dim, nPts) = size(v,1,2)

    if size(P) == (dim+1, dim+1)   # Assume data is inhomogeneous
        t = P[1:dim, 1:dim]*v
        scale = P[dim+1,1:dim]*v

        # Add the translation component of the transform to the scale
        # then add the translation component to t and divide by the
        # scale to normalise the result.
        if checkinf   # Only normalise if the scale parameter is non-zero
            for j = 1:nPts
                scale[j] += P[dim+1,dim+1]
                if abs(scale[j]) > eps()
                    for i = 1:dim
                        t[i,j] = (t[i,j] + P[i,dim+1])/scale[j]
                    end
                end
            end

        else   # Do not check for points at infinity
            for j = 1:nPts
                scale[j] += P[dim+1,dim+1]

                for i = 1:dim
                    t[i,j] = (t[i,j] + P[i,dim+1])/scale[j]
                end
            end
        end

    elseif size(P) == (dim, dim)  # Homogeneous case
        t = P*v

        # Normalise to scale of 1
        if checkinf   # Only normalise if the scale parameter is non-zero
            for i = 1:nPts
                if abs(t[dim,i]) > eps()
                    for r = 1:dim-1
                        t[r,i] /= t[dim,i]
                    end
                    t[dim,i] = 1.0
                end
            end

        else   # Do not check for points at infinity
            for i = 1:nPts
                for r = 1:dim-1
                    t[r,i] /= t[dim,i]
                end
                t[dim,i] = 1.0
            end
        end

    else
        error("Transformation matrix and point data dimensions do not match")
    end

    return t
end

#---------------------------------------------------------------------
"""
inveuler - Inverse of Euler transform
```
Usage:  (euler1, euler2) = inveuler(T)

Argument:  T - 4x4 Homogeneous transformation matrix or 3x3 rotation matrix

Returns: euler1 = [phi1, theta1, psi1] - the 1st solution and,
         euler2 = [phi2, theta2, psi2] - the 2nd solution

     rotz(phi1)   * roty(theta1)    * rotz(psi1)      = T
  rotz(euler1[1]) * roty(euler1[2]) * rotz(euler1[3]) = T
```
See also: invrpy(), invht(), rotx(), roty(), rotz()
"""
function inveuler(T::Array{T1,2}) where T1 <: Real
    #=
    Reference: Richard P. Paul
    Robot Manipulators: Mathematics, Programming and Control.
    MIT Press 1981. Page 68
    =#

    phi1 = atan(T[2,3], T[1,3])
    phi2 = phi1 + pi

    theta1 = atan(cos(phi1)*T[1,3] + sin(phi1)*T[2,3], T[3,3])
    theta2 = atan(cos(phi2)*T[1,3] + sin(phi2)*T[2,3], T[3,3])

    psi1 = atan(-sin(phi1)*T[1,1] + cos(phi1)*T[2,1],
                 -sin(phi1)*T[1,2] + cos(phi1)*T[2,2])
    psi2 = atan(-sin(phi2)*T[1,1] + cos(phi2)*T[2,1],
                 -sin(phi2)*T[1,2] + cos(phi2)*T[2,2])

    euler1 = [phi1, theta1, psi1]
    euler2 = [phi2, theta2, psi2]
    return (euler1, euler2)
end

#---------------------------------------------------------------------
"""
invrpy - Inverse of Roll Pitch Yaw transform
```
Usage:  (rpy1, rpy2) = invrpy(RPY)

Argument:  RPY - 4x4 Homogeneous transformation matrix or 3x3 rotation matrix

Returns:  rpy1 = [phi1, theta1, psi1] - the 1st solution and
          rpy2 = [phi2, theta2, psi2] - the 2nd solution

   rotz(phi1)  * roty(theta1)  * rotx(psi1)    = RPY
 rotz(rpy1[1]) * roty(rpy1[2]) * rotx(rpy1[3]) = RPY
```
See also: inveuler(), invht(), rotx(), roty(), rotz()
"""
function invrpy(RPY::Array{T,2}) where T <: Real
    #=
    Reference: Richard P. Paul
    Robot Manipulators: Mathematics, Programming and Control.
    MIT Press 1981. Page 70
    =#
    # Z rotation
    phi1 = atan(RPY[2,1], RPY[1,1])
    phi2 = phi1 + pi

    # Y rotation
    theta1 = atan(-RPY[3,1], cos(phi1)*RPY[1,1] + sin(phi1)*RPY[2,1])
    theta2 = atan(-RPY[3,1], cos(phi2)*RPY[1,1] + sin(phi2)*RPY[2,1])

    # X rotation
    psi1 = atan(sin(phi1)*RPY[1,3] - cos(phi1)*RPY[2,3],
                -sin(phi1)*RPY[1,2] + cos(phi1)*RPY[2,2]);
    psi2 = atan(sin(phi2)*RPY[1,3] - cos(phi2)*RPY[2,3],
                -sin(phi2)*RPY[1,2] + cos(phi2)*RPY[2,2]);

    rpy1 = [phi1, theta1, psi1]
    rpy2 = [phi2, theta2, psi2]
    return (rpy1, rpy2)
end

#---------------------------------------------------------------------
"""
matrix2angleandaxis - Decompose homogeneous matrix to angle and axis
```
Usage: (ang, axis) = matrix2angleandaxis(T)

Argument:   T - 4x4 Homogeneous transformation matrix, or 3x3 rotation matrix.
Returns:  ang - Rotation angle in radians.
         axis - Unit 3-vector defining the rotation axis.
```
Note that only the top left 3x3 rotation component of T is used, any
translation component in T is ignored.

See also: angleaxis2matrix(), angleaxisrotate(), angleaxis(), normaliseangleaxis()
"""
function matrix2angleandaxis(T::Array)

    # This code follows the implementation suggested by Hartley and Zisserman
    R = T[1:3, 1:3]   # Extract rotation part of T

    # Trap case where rotation is very small.  ( See angleaxis2matrix() )
    Reye = R - I
    if norm(Reye) < 1e-8
        t = [T[3,2], T[1,3], T[2,1]]
        return norm(Reye), t
    end

    # Otherwise find rotation axis as the eigenvector having unit eigenvalue
    # Solve (R-I)v = 0
#    (eigval, eigvec) = eig(Reye)
    F = eigen(Reye)

    # Find index of eigenvalue with smallest magnitude, the
    # corresponding eigenvector will be the axis.
    absval = abs.(F.values)
    ind = sortperm(absval)          # Find index of smallest one
    if absval[ind[1]] > 0.001       # Hopefully it is close to 0
        @warn("Rotation matrix is dubious")
    end

    axis = real(F.vectors[:,ind[1]])    # Extract appropriate eigenvector

    if abs(norm(axis) - 1.0) > .0001    # Debug
        @warn("Non unit rotation axis")
    end

    # Now determine the rotation angle
    twocostheta = tr(R)-1.0
    twosinthetav = [R[3,2]-R[2,3], R[1,3]-R[3,1], R[2,1]-R[1,2]]
    twosintheta = dot(axis, twosinthetav)

    # Note use of twosinetheta to convert it from a 1x1 matrix to a scalar
    theta = atan(twosintheta, twocostheta)

    return theta, axis
end

#---------------------------------------------------------------------
"""
matrix2angleaxis - Homogeneous matrix to angle-axis description
```
Usage: t = matrix2angleaxis(T)

Argument:   T - 4x4 Homogeneous transformation matrix, or 3x3 rotation matrix.
Returns:    t - 3x1 column vector giving rotation axis with magnitude equal
                to the rotation angle in radians.
```
Note that only the top left 3x3 rotation component of T is used, any
translation component in T is ignored.

See also: angleaxis2matrix(),  angleaxisrotate(), angleaxis(), normaliseangleaxis()

"""
function matrix2angleaxis(T::Array)
    (theta, axis) = matrix2angleandaxis(T)
    return theta*axis
end

#---------------------------------------------------------------------
"""
matrix2quaternion - Homogeneous matrix to quaternion

Converts 4x4 homogeneous rotation matrix to quaternion
```
Usage: Q = matrix2quaternion(T)

Argument:   T - 4x4 Homogeneous transformation matrix
Returns:    Q - 4-vector quaternion in the form [w, xi, yj, zk]
```
See also: quaternion2matrix()
"""
function matrix2quaternion(T::Array)
    (theta, axis) = matrix2angleandaxis(T)
    return [cos(theta/2); axis*sin(theta/2)]
end

#---------------------------------------------------------------------
"""
angleaxis - Constructs angle-axis descriptor
```
Usage: t = angleaxis(theta, axis)

Arguments: theta - angle of rotation.
           axis  - 3x1 vector defining axis of rotation.
Returns:   t     - 3-vector giving rotation axis with magnitude equal to the
                   rotation angle in radians.
```
See also: matrix2angleaxis(), angleaxisrotate(), angleaxis2matrix(), normaliseangleaxis()
"""
function angleaxis(theta::Real, axis::Vector)

    if length(axis) != 3
        error("Axis must be a 3 vector or array")
    end

    ax = normalize(axis)  # Ensure unit magnitude

    # Normalise theta to lie in the range -pi to pi to ensure one-to-one mapping
    # between angle-axis descriptor and resulting rotation.
    theta = rem(theta, 2*pi)  # Remove multiples of 2pi

    if theta > pi
        theta = theta - 2*pi
    elseif  theta < -pi
        theta = theta + 2*pi
    end

    return theta*ax
end

#---------------------------------------------------------------------
"""
quaternion  - Construct quaternion

```
Usage:  Q = quaternion(theta, axis)

Arguments: theta - angle of rotation
           axis  - 3-vector defining axis of rotation.
Returns:   Q     - a 4-vector quaternion in the form [w, xi, yj, zk]
```
See also:  quaternion2matrix(), matrix2quaternion(), quaternionrotate()
"""
function quaternion(theta::Real, axis::Vector)

    if length(axis) != 3
        error("Axis must be a 3 vector or array")
    end

    Q = zeros(4)
    Q[1] = cos(theta/2)
    Q[2:4] = sin(theta/2)*normalize(axis)
    return Q
end

#---------------------------------------------------------------------
"""
normaliseangleaxis - Normalises angle-axis descriptor.

Function normalises theta so that it lies in the range -pi to pi to ensure
one-to-one mapping between angle-axis descriptor and resulting rotation.
```
Usage: t2 = normaliseangleaxis(t)

Argument:   t  - 3-vector giving rotation axis with magnitude equal to the
                 rotation angle in radians.
Returns:    t2 - Normalised angle-axis descriptor
```
See also: matrix2angleaxis(), angleaxis(), angleaxis2matrix(), angleaxisrotate()
"""
function normaliseangleaxis(t::Vector)

    if length(t) != 3
        error("Axis must be a 3-vector")
    end

    theta = norm(t)
    axis = t/theta

    theta = rem(theta, 2pi)  # Remove multiples of 2pi

    if theta > pi              # Note theta cannot be -ve
        theta = theta - 2pi;
    end

    return theta*axis
end

#---------------------------------------------------------------------
"""
quaternion2matrix - Quaternion to a 4x4 homogeneous transformation matrix
```
Usage:  T = quaternion2matrix(Q)

Argument: Q - a quaternion in the form [w, xi, yj, zk].
Returns:  T - 4x4 Homogeneous rotation matrix.
```
See also: matrix2quaternion(), quaternion(), quaternionrotate()
"""
function quaternion2matrix(Qu::Vector)

    if length(Qu) != 4
        error("Quaternion must be a 4 vector")
    end

    Q = normalize(Qu)    # Ensure Q has unit norm

    # Set up convenience variables
    w = Q[1]; x = Q[2]; y = Q[3]; z = Q[4]
    w2 = w^2; x2 = x^2; y2 = y^2; z2 = z^2
    xy = x*y; xz = x*z; yz = y*z
    wx = w*x; wy = w*y; wz = w*z

    T = [w2+x2-y2-z2  2*(xy - wz)  2*(wy + xz)   0.0
         2*(wz + xy)  w2-x2+y2-z2  2*(yz - wx)   0.0
         2*(xz - wy)  2*(wx + yz)  w2-x2-y2+z2   0.0
             0.0           0.0          0.0      1.0]

    return T
end

#---------------------------------------------------------------------
"""
quaternionconjugate - Conjugate of a quaternion
```
Usage: Qconj = quaternionconjugate(Q)

Argument: Q     - Quaternions in the form  Q = [Qw, Qi, Qj, Qk].
Returns:  Qconj - Conjugate.
```
See also: quaternion(), quaternionrotate(), quaternionproduct()
"""
function quaternionconjugate(Q::Vector)

    if length(Q) != 4
        error("Quaternion must be a 4 vector")
    end

    return [Q[1], -Q[2], -Q[3], -Q[4]]
end

#---------------------------------------------------------------------
"""
quaternionproduct - Computes product of two quaternions
```
Usage: Q = quaternionproduct(A, B)

Arguments: A, B - Quaternions assumed to be 4-vectors in the
                  form  A = [Aw, Ai, Aj, Ak].
Returns:   Q    - Quaternion product.
```
See also: quaternion(), quaternionrotate(), quaternionconjugate()
"""
function  quaternionproduct(A::Vector, B::Vector)

    if length(A) != 4 || length(B) != 4
        error("Quaternions must be 4 vectors")
    end

    Q = zeros(4,1)
    Q[1]  =  A[1]*B[1]  -  A[2]*B[2]  -  A[3]*B[3]  -  A[4]*B[4]
    Q[2]  =  A[1]*B[2]  +  A[2]*B[1]  +  A[3]*B[4]  -  A[4]*B[3]
    Q[3]  =  A[1]*B[3]  -  A[2]*B[4]  +  A[3]*B[1]  +  A[4]*B[2]
    Q[4]  =  A[1]*B[4]  +  A[2]*B[3]  -  A[3]*B[2]  +  A[4]*B[1]

    return Q
end

#---------------------------------------------------------------------
"""
quaternionrotate - Rotates a 3D vector by a quaternion
```
Usage:   vnew = quaternionrotate(Q, v)

Arguments: Q - a quaternion in the form [w, xi, yj, zk]
           v - a vector to rotate, either an inhomogeneous 3-vector or a
               homogeneous 4-vector.
Returns:   vnew - rotated vector.
```
See also: matrix2quaternion(), quaternion2matrix(), quaternion()
"""
function quaternionrotate(Q::Vector, v::Vector)
    #=
    % Code forms the equivalent 3x3 rotation matrix from the quaternion and
    % applies it to a vector
    %
    % Note that Qw^2 + Qi^2 + Qj^2 + Qk^2 = 1
    % So the top-left entry of the rotation matrix of
    %   Qw^2 + Qi^2 - Qj^2 - Qk^2
    % can be rewritten as
    %   Qw^2 + Qi^2 + Qj^2 + Qk^2 - 2Qj^2 - 2Qk^2
    % = 1 - 2Qj^2 - 2Qk^2
    %
    % Similar optimization applies to the other diagonal elements
    =#

    if length(Q) != 4
        error("Quaternion must be a 4 vector")
    end

    if length(v) != 4 && length(v) != 3
        error("Input vector must be of length 3 or 4")
    end

    # Copy v to vnew to allocate space.  If v is a 4 element
    # homogeneous vector this also sets the homogeneous scale factor
    # of vnew. Ensure vnew is float so that we avoid Inexact errors
    vnew = float(copy(v))

    Qw = Q[1];  Qi = Q[2];  Qj = Q[3];  Qk = Q[4]

    t2 =   Qw*Qi
    t3 =   Qw*Qj
    t4 =   Qw*Qk
    t5 =  -Qi*Qi
    t6 =   Qi*Qj
    t7 =   Qi*Qk
    t8 =  -Qj*Qj
    t9 =   Qj*Qk
    t10 = -Qk*Qk
    vnew[1] = 2*( (t8 + t10)*v[1] + (t6 -  t4)*v[2] + (t3 + t7)*v[3] ) + v[1]
    vnew[2] = 2*( (t4 +  t6)*v[1] + (t5 + t10)*v[2] + (t9 - t2)*v[3] ) + v[2]
    vnew[3] = 2*( (t7 -  t3)*v[1] + (t2 +  t9)*v[2] + (t5 + t8)*v[3] ) + v[3]

    return vnew
end

#---------------------------------------------------------------------
"""
vector2quaternion - Embeds 3-vector in a quaternion representation
```
Usage: Q = vector2quaternion(v)

Argument:  v - 3-vector.
Returns:   Q - Quaternion given by [0; v]
```
See also: quaternion(), quaternionrotate(), quaternionproduct(), quaternionconjugate()
"""
function vector2quaternion(v::Vector)

  if length(v) != 3
      error("v must be a 3-vector")
  end

  return [0.0, v[1], v[2], v[3]]
end


text:
using Test

println("Starting tests")

@time include("QuarticSolverTests.jl")

ext: jl

max_stars_repo_name: EP-Guy/QuarticSolver

max_issues_repo_name: EP-Guy/QuarticSolver

max_forks_repo_name: EP-Guy/QuarticSolver


text:
# include("../../src/SFFM.jl")
using LinearAlgebra, Plots, JLD2, GLM, SFFM

## define the model(s)
include("exampleModelDef.jl")

## load sims
@load pwd()*"/examples/paperNumerics/dump/sims.jld2" sims

## analytic X distribution for comparison
# construction
Ψₓ = SFFM.PsiFunX( approxModel)
ξₓ = SFFM.MakeXiX( approxModel, Ψₓ)
pₓ, πₓ, Πₓ, Kₓ = SFFM.StationaryDistributionX( approxModel, Ψₓ, ξₓ)

## section 4.5: error for approximation of stationary distribution of X
Basis = "lagrange"
Δs = [1.6; 0.8; 0.4]
NBasesRange = [1; 3; 5; 7]

πnorms = zeros(length(Δs), length(NBasesRange))
Ψnorms = zeros(length(Δs), length(NBasesRange))
times = Array{Any}(undef, length(Δs), length(NBasesRange))
gctimes = Array{Any}(undef, length(Δs), length(NBasesRange))
mems = Array{Any}(undef, length(Δs), length(NBasesRange))
alloc = Array{Any}(undef, length(Δs), length(NBasesRange))
approxSpec = Array{Any}(undef, length(Δs), length(NBasesRange))

πnormsme = zeros(length(Δs), length(NBasesRange))
Ψnormsme = zeros(length(Δs), length(NBasesRange))
timesme = Array{Any}(undef, length(Δs), length(NBasesRange))
gctimesme = Array{Any}(undef, length(Δs), length(NBasesRange))
memsme = Array{Any}(undef, length(Δs), length(NBasesRange))
allocme = Array{Any}(undef, length(Δs), length(NBasesRange))
approxSpecme = Array{Any}(undef, length(Δs), length(NBasesRange))

πnormsfv = zeros(length(Δs), length(NBasesRange))
Ψnormsfv = zeros(length(Δs), length(NBasesRange))
timesfv = Array{Any}(undef, length(Δs), length(NBasesRange))
gctimesfv = Array{Any}(undef, length(Δs), length(NBasesRange))
memsfv = Array{Any}(undef, length(Δs), length(NBasesRange))
allocfv = Array{Any}(undef, length(Δs), length(NBasesRange))
approxSpecfv = Array{Any}(undef, length(Δs), length(NBasesRange))

for d = 1:length(Δs), n = 1:length(NBasesRange)
    # define the mesh for each iteration
    nBases = NBasesRange[n]
    Δtemp = Δs[d]
    println("mesh details; h = ", Δtemp, " NBases = ", nBases)

    # collect time and memory stats
    ~, times[d, n], mems[d, n], gctimes[d, n], alloc[d, n] = @timed begin
        Nodes = collect(approxBounds[1, 1]:Δtemp:approxBounds[1, 2])
        mesh = SFFM.DGMesh(
            approxModel,
            Nodes,
            nBases,
            Basis =  "lagrange",
        )
        approxSpec[d, n] = (Δtemp, nBases, SFFM.TotalNBases(mesh) * SFFM.NPhases(approxModel))

        # compute the marginal via DG
        All = SFFM.MakeAll( approxModel, mesh, approxType = "interpolation")
        Ψ = SFFM.PsiFun( All.D)

        # the distribution of X when Y first returns to 0
        ξ = SFFM.MakeXi( All.B, Ψ)

        marginalX, p, K = SFFM.MakeLimitDistMatrices(
            All.B,
            All.D,
            All.R.RDict,
            Ψ,
            ξ,
            mesh,
            approxModel,
        )
    end

    # collect time and memory stats
    ~, timesme[d, n], memsme[d, n], gctimesme[d, n], allocme[d, n] = @timed begin
        Nodes = collect(approxBounds[1, 1]:Δtemp:approxBounds[1, 2])
        mesh = SFFM.FRAPMesh(
            approxModel,
            Nodes,
            nBases,
        )
        dgmesh = SFFM.DGMesh(
            approxModel,
            Nodes,
            nBases,
            Basis =  "lagrange",
        )
        approxSpecme[d, n] = (Δtemp, nBases, SFFM.TotalNBases(mesh) * SFFM.NPhases(approxModel))

        # compute the marginal via FRAP approx
        # me = SFFM.MakeME(SFFM.CMEParams[nBases], mean = SFFM.Δ(mesh)[1])
        B = SFFM.MakeFullGenerator( approxModel, mesh)
        R = SFFM.MakeR( approxModel, dgmesh, approxType = "interpolation")
        D = SFFM.MakeD( mesh, B, R)
        Ψme = SFFM.PsiFun(D)

        # the distribution of X when Y first returns to 0
        ξme = SFFM.MakeXi( B, Ψme)

        marginalXme, pme, Kme = SFFM.MakeLimitDistMatrices(
            B,
            D,
            R.RDict,
            Ψme,
            ξme,
            mesh,
            approxModel,
        )
    end

    # collect time and memory stats
    ~, timesfv[d, n], memsfv[d, n], gctimesfv[d, n], allocfv[d, n] = @timed begin
        Nodes = collect(approxBounds[1, 1]:Δtemp:approxBounds[1, 2])
        mesh = SFFM.FVMesh(
            approxModel,
            Nodes,
        )
        approxSpecfv[d, n] = (Δtemp, nBases, SFFM.TotalNBases(mesh) * SFFM.NPhases(approxModel))

        # compute the marginal via FRAP approx
        B = SFFM.MakeB(approxModel, mesh, nBases)
        R = SFFM.MakeR( approxModel, mesh, approxType = "interpolation")
        D = SFFM.MakeD( mesh, B, R)
        Ψfv = SFFM.PsiFun(D)

        # the distribution of X when Y first returns to 0
        ξfv = SFFM.MakeXi( B, Ψfv)

        marginalXfv, pfv, Kfv = SFFM.MakeLimitDistMatrices(
            B,
            D,
            R.RDict,
            Ψfv,
            ξfv,
            mesh,
            approxModel,
        )
    end
    dgmesh = SFFM.DGMesh(
            approxModel,
            Nodes,
            nBases,
            Basis =  "lagrange",
        )
    frapmesh = SFFM.FRAPMesh(
            approxModel,
            Nodes,
            nBases,
        )
    fvmesh = SFFM.FVMesh(
        approxModel,
        Nodes,
    )
    # convert marginalX to a distribution for analysis
    DGStationaryDist = SFFM.Coeffs2Dist(
        approxModel,
        dgmesh,
        marginalX,
        SFFM.SFFMProbability,
    )
    meStationaryDist = SFFM.Coeffs2Dist(
        approxModel,
        frapmesh,
        marginalXme,
        SFFM.SFFMProbability,
    )
    fvStationaryDist = SFFM.Coeffs2Dist(
        approxModel,
        fvmesh,
        marginalXfv,
        SFFM.SFFMProbability,
    )

    ## stationary distirbution stuff
    # evaluate the analytic result given the mesh
    analyticX = SFFM.SFFMProbability(
        [pₓ[:]; 0; 0],
        Πₓ(Matrix(mesh.Nodes[2:end]')) - Πₓ(Matrix(mesh.Nodes[1:end-1]')),
        mesh.Nodes[1:end-1] + SFFM.Δ(mesh) / 2,
    )

    # save them
    πnorms[d, n] = SFFM.starSeminorm( DGStationaryDist, analyticX)
    πnormsme[d, n] = SFFM.starSeminorm( meStationaryDist, analyticX)
    πnormsfv[d, n] = SFFM.starSeminorm( fvStationaryDist, analyticX)
end

let p = plot()
    types = [:solid,:dash,:dashdot,:dashdotdot]
    for n = 1:length(NBasesRange)
        p = plot!(
            Δs,
            πnorms[:, n],
            xaxis = :log,
            yaxis = :log,
            xlabel = "h",
            ylabel = "error",
            label = "DG: N_k = " * string(NBasesRange[n]),
            legend = :outertopright,
            linestyle = types[n],
            seriestype = :line,
            # markershape = :auto,
            title = "Error",
            grid = false,
            color = :red,
        )
        p = plot!(
            Δs,
            πnormsme[:, n],
            label = "ME: N_k = " * string(NBasesRange[n]),
            linestyle = types[n],
            seriestype = :line,
            # markershape = :auto,
            color = :blue,
        )
        p = plot!(
            Δs,
            πnormsfv[:, n],
            label = "FV: N = " * string(NBasesRange[n]),
            linestyle = types[n],
            seriestype = :line,
            # markershape = :auto,
            color = :black,
        )
    end
    display(p)
    # savefig(pwd()*"/examples/paperNumerics/dump/piErrorVsDelta.png")
end

let p = plot()
    types = [:solid,:dash,:dashdot,:dashdotdot,:solid,:dash]
    for d = 1:length(Δs)
        p = plot!(
            NBasesRange,
            πnorms[d, :],
            yaxis = :log,
            xlabel = "N_k",
            ylabel = "error",
            label = "DG: h = " * string(Δs[d]),
            legend = :outertopright,
            linestyle = types[d],
            seriestype = :line,
            # markershape = :auto,
            title = "Error",
            grid = false,
            color = :red,
        )
        p = plot!(
            NBasesRange,
            πnormsme[d, :],
            label = "ME: h = " * string(Δs[d]),
            linestyle = types[d],
            seriestype = :line,
            # markershape = :auto,
            colot = :blue,
        )
        p = plot!(
            NBasesRange,
            πnormsfv[d, :],
            label = "FV: h = " * string(Δs[d]),
            linestyle = types[d],
            seriestype = :line,
            # markershape = :auto,
            color = :black,
        )
    end
    display(p)
    # savefig(pwd()*"/examples/paperNumerics/dump/piErrorVsNBases.png")
end

begin
    n = 1
    data = log.(πnorms[:,n])
    println("DG rate of convergence for ", NBasesRange[n], " basis functions is ")
    display(
        lm([ones(length(Δs)) log.(Δs)], data)
    )
    println("")

    data = log.(πnormsme[:,n])
    ind = .!isnan.(data)
    data = data[ind]
    println("ME rate of convergence for ", NBasesRange[n], " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    data = log.(πnormsfv[:,n])
    ind = .!isnan.(data)
    data = data[ind]
    println("FV rate of convergence for ", NBasesRange[n], " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    n = 2
    data = log.(πnorms[:,n])
    println("DG rate of convergence for ", NBasesRange[n], " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    data = log.(πnormsme[:,n])
    println("ME rate of convergence for ", NBasesRange[n], " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    data = log.(πnormsfv[:,n])
    println("FV rate of convergence for ", NBasesRange[n], " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    n = 3
    data = log.(πnorms[:,n])
    println("DG rate of convergence for ", NBasesRange[n], " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    data = log.(πnormsme[:,n])
    println("ME rate of convergence for ", n, " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    data = log.(πnormsfv[:,n])
    println("FV rate of convergence for ", n, " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")
end

for d = 1:length(Δs)
    data = log.(πnorms[d,:])
    ind = .!isnan.(data)
    data = data[ind]
    println("DG rate of convergence for h = ", Δs[d])
    display(
        lm([ones(length(NBasesRange[ind])) NBasesRange[ind]], data)
    )
    println("")
end

for d = 1:length(Δs)
    data = log.(πnormsme[d,:])
    ind = .!isnan.(data)
    data = data[ind]
    println("ME rate of convergence for h = ", Δs[d])
    display(
        lm([ones(length(NBasesRange[ind])) NBasesRange[ind]], data)
    )
    println("")
end

for d = 1:length(Δs)
    data = log.(πnormsfv[d,:])
    ind = .!isnan.(data)
    data = data[ind]
    println("FV rate of convergence for h = ", Δs[d])
    display(
        lm([ones(length(NBasesRange[ind])) NBasesRange[ind]], data)
    )
    println("")
end

mp = (log.(πnorms[1:end-1, :]) .- log.(πnorms[2:end, :])) ./
    (log.(Δs[1:end-1]) - log.(Δs[2:end]))
mpme = (log.(πnormsme[1:end-1, :]) .- log.(πnormsme[2:end, :])) ./
    (log.(Δs[1:end-1]) - log.(Δs[2:end]))
mpfv = (log.(πnormsfv[1:end-1, :]) .- log.(πnormsfv[2:end, :])) ./
    (log.(Δs[1:end-1]) - log.(Δs[2:end]))
np = (log.(πnorms[:, 1:end-1]) .- log.(πnorms[:, 2:end])) ./
    ((NBasesRange[1:end-1]) - (NBasesRange[2:end]))'
npme = (log.(πnormsme[:, 1:end-1]) .- log.(πnormsme[:, 2:end])) ./
    ((NBasesRange[1:end-1]) - (NBasesRange[2:end]))'
npme = (log.(πnormsfv[:, 1:end-1]) .- log.(πnormsfv[:, 2:end])) ./
    ((NBasesRange[1:end-1]) - (NBasesRange[2:end]))'

stats = Array{Any}(undef,length(Δs)+1,length(NBasesRange)+1)
statsme = Array{Any}(undef,length(Δs)+1,length(NBasesRange)+1)
statsfv = Array{Any}(undef,length(Δs)+1,length(NBasesRange)+1)
for n in 1:length(NBasesRange), d in 1:length(Δs)
        stats[d+1,n+1] = (πnorms[d,n], times[d,n], mems[d,n]/(1024*1024), approxSpec[d,n][3])
        statsme[d+1,n+1] = (πnormsme[d,n], timesme[d,n], memsme[d,n]/(1024*1024), approxSpecme[d,n][3])
        statsfv[d+1,n+1] = (πnormsfv[d,n], timesfv[d,n], memsfv[d,n]/(1024*1024), approxSpecfv[d,n][3])
end
stats[1,2:end] = NBasesRange
stats[2:end,1] = Δs
statsme[1,2:end] = NBasesRange
statsme[2:end,1] = Δs
statsfv[1,2:end] = NBasesRange
statsfv[2:end,1] = Δs

using NumericIO
innards = stats[2:end,2:end]
innardsme = statsme[2:end,2:end]
innardsfv = statsfv[2:end,2:end]
asciiexponentfmt = NumericIO.IOFormattingExpNum(
	"\\times10^{", false, '+', '-', NumericIO.ASCII_SUPERSCRIPT_NUMERALS
)
fmt = NumericIO.IOFormattingReal(asciiexponentfmt,
	ndigits=3, decpos=0, decfloating=true, eng=false, minus='-', inf="Inf"
)

begin
    println("\\begin{tabular}{ | c | l | l | l | l | l | }")
    println("\\hline")
    println("\\multirow{2}{*}{\\(\\Delta\\)} & \\multirow{2}{*}{}& \\multicolumn{4}{c|}{Number of basis functions, \\(N_k\\)} \\\\ \\cline{3-6}
          &      & 1            & 3           & 5           & 7           \\\\\\hline"
        )
    for d in 1:length(Δs)
        println(Δs[d], " &")
        println("
        \\(\\begin{array}{l}
            \\mbox{ }\\\\
            \\mbox{error}
            \\\\ \\mbox{time (sec)}
            \\\\ \\mbox{memory (MB)}
            \\\\ n_\\varphi
        \\end{array}\\)"
        )
        for n in 1:length(NBasesRange)
            println("&")
                println("
                \\(\\begin{array}{l}
                    \\mbox{DG}\\\\
                    ", formatted(innards[d,n][1],fmt),"}
                    \\\\",round(innards[d,n][2],sigdigits=3),"
                    \\\\",round(innards[d,n][3],sigdigits=3),"
                    \\\\",convert(Int,round(innards[d,n][4],digits=0)),"
                \\end{array}\\)\\(\\begin{array}{l}
                \\mbox{ME}\\\\
                ", formatted(innardsme[d,n][1],fmt),"}
                \\\\",round(innardsme[d,n][2],sigdigits=3),"
                \\\\",round(innardsme[d,n][3],sigdigits=3),"
                \\\\",convert(Int,round(innardsme[d,n][4],digits=0)),"
            \\end{array}\\)\\(\\begin{array}{l}
            \\mbox{FV}\\\\
            ", formatted(innardsfv[d,n][1],fmt),"}
            \\\\",round(innardsfv[d,n][2],sigdigits=3),"
            \\\\",round(innardsfv[d,n][3],sigdigits=3),"
            \\\\",convert(Int,round(innardsfv[d,n][4],digits=0)),"
        \\end{array}\\)"
                )
        end
        println("\\\\ \\hline")
    end
    println("\\end{tabular}")
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
# the idea of this is to find a basis of the
# kernel of matrix. This is done using the eigenvalues.
function kernel_basis(A, tol = 1e-5)
	B = Array{Float64,1}[]
	kernel = nullspace(A);
	kernel_size = size(kernel)[2]
	for n = 1:kernel_size
		push!(B,kernel[:,n]);
	end
	return B
end


text:
"""
    label_propagation(g, maxiter=1000)

Community detection using the label propagation algorithm.
Return two vectors: the first is the label number assigned to each node, and
the second is the convergence history for each node. Will return after
`maxiter` iterations if convergence has not completed.

### References
- [Raghavan et al.](http://arxiv.org/abs/0709.2938)
"""
function label_propagation(g::AbstractGraph, maxiter=1000)
    T = eltype(g)
    n = nv(g)
    label = collect(one(T):n)
    active_vs = IntSet(vertices(g))
    c = NeighComm(collect(one(T):n), fill(-1,n), one(T))
    convergence_hist = Vector{Int}()
    random_order = Vector{T}(n)
    i = 0
    while !isempty(active_vs) && i < maxiter
        num_active = length(active_vs)
        push!(convergence_hist, num_active)
        i += 1

        # processing vertices in random order
        for (j,node) in enumerate(active_vs)
            random_order[j] = node
        end
        range_shuffle!(1:num_active, random_order)
        @inbounds for j=1:num_active
            u = random_order[j]
            old_comm = label[u]
            label[u] = vote!(g, label, c, u)
            if old_comm != label[u]
                for v in out_neighbors(g, u)
                    push!(active_vs, v)
                end
            else
                delete!(active_vs, u)
            end
        end
    end
    fill!(c.neigh_cnt, 0)
    renumber_labels!(label, c.neigh_cnt)
    label, convergence_hist
end

"""
    NeighComm{T}

Type to record neighbor labels and their counts.
"""
mutable struct NeighComm{T<:Integer}
    neigh_pos::Vector{T}
    neigh_cnt::Vector{Int}
    neigh_last::T
end

"""
    range_shuffle!(r, a)

Fast shuffle Array `a` in UnitRange `r`.
"""
function range_shuffle!(r::UnitRange, a::AbstractVector)
    (r.start > 0 && r.stop <= length(a)) || error("out of bounds")
    @inbounds for i=length(r):-1:2
        j = rand(1:i)
        ii = i + r.start - 1
        jj = j + r.start - 1
        a[ii],a[jj] = a[jj],a[ii]
    end
end

"""
    vote!(g, m, c, u)

Return the label with greatest frequency.
"""
function vote!(g::AbstractGraph, m::Vector, c::NeighComm, u::Integer)
    @inbounds for i=1:c.neigh_last-1
        c.neigh_cnt[c.neigh_pos[i]] = -1
    end
    c.neigh_last = 1
    c.neigh_pos[1] = m[u]
    c.neigh_cnt[c.neigh_pos[1]] = 0
    c.neigh_last = 2
    max_cnt = 0
    for neigh in out_neighbors(g,u)
        neigh_comm = m[neigh]
        if c.neigh_cnt[neigh_comm] < 0
            c.neigh_cnt[neigh_comm] = 0
            c.neigh_pos[c.neigh_last] = neigh_comm
            c.neigh_last += 1
        end
        c.neigh_cnt[neigh_comm] += 1
        if c.neigh_cnt[neigh_comm] > max_cnt
            max_cnt = c.neigh_cnt[neigh_comm]
        end
    end
    # ties breaking randomly
    range_shuffle!(1:c.neigh_last-1, c.neigh_pos)
    for lbl in c.neigh_pos
        if c.neigh_cnt[lbl] == max_cnt
            return lbl
        end
    end
end

function renumber_labels!(membership::Vector, label_counters::Vector{Int})
    N = length(membership)
    (maximum(membership) > N || minimum(membership) < 1) && error("Label must between 1 and |V|")
    j = 1
    @inbounds for i=1:length(membership)
        k = membership[i]
        if k >= 1
            if label_counters[k] == 0
                # We have seen this label for the first time
                label_counters[k] = j
                k = j
                j += 1
            else
                k = label_counters[k]
            end
        end
        membership[i] = k
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using ADCME
using MAT
using PyPlot
np = pyimport("numpy")

# Domain information
NT = 10
Δt = 1/NT
n = 20
m = n
h = 1.0/n
bdnode = Int64[]
for i = 1:m+1
    for j = 1:n+1
        if i==1 || i==m+1 || j==1|| j==n+1
            push!(bdnode, (j-1)*(m+1)+i)
        end
    end
end

is_training = false
b = 1.0
# Physical parameters
if !is_training
    E = 1.0
    ν = 0.35
    global H = E/(1-ν^2)*[1.0 ν 0
        ν 1 0
        0 0 1-ν]|>constant
else
    H = diagm(0=>ones(3))
    global H = Variable(H)
end
Q = SparseTensor(compute_fvm_tpfa_matrix(m, n, h))
K = compute_fem_stiffness_matrix(H, m, n, h)
L = SparseTensor(compute_interaction_matrix(m, n, h))
M = SparseTensor(compute_fvm_mass_matrix(m, n, h))
A = [K -b*L'
b*L/Δt 1/Δt*M-Q]
A, Abd = fem_impose_coupled_Dirichlet_boundary_condition(A, bdnode, m, n, h)
# error()
U = zeros(m*n+2(m+1)*(n+1), NT+1)
x = Float64[]; y = Float64[]
for j = 1:n+1
    for i = 1:m+1
        push!(x, (i-1)*h)
        push!(y, (j-1)*h)
    end
end

injection = (div(n,2)-1)*m + 3
production = (div(n,2)-1)*m + m-3

function condition(i, ta_u)
    i<=NT
end

function body(i, ta_u)
    u = read(ta_u, i)
    bdval = zeros(2*length(bdnode))
    rhs1 = vector([bdnode; bdnode.+ (m+1)*(n+1)], bdval, 2(m+1)*(n+1))
    rhs2 = zeros(m*n)
    rhs2[injection] += 1.0
    rhs2[production] -= 1.0
    rhs2 += b*L*u[1:2(m+1)*(n+1)]/Δt +
            M * u[2(m+1)*(n+1)+1:end]/Δt
    rhs = [rhs1;rhs2]
    rhs -= Abd * bdval
    o = A\rhs
    ta_u = write(ta_u, i+1, o)
    i+1, ta_u
end

i = constant(1, dtype=Int32)
ta_u = TensorArray(NT+1)
ta_u = write(ta_u, 1, constant(zeros(2(m+1)*(n+1)+m*n)))
_, u_out = while_loop(condition, body, [i, ta_u]; parallel_iterations=1)
u_out = stack(u_out)
u_out.set_shape((NT+1, size(u_out,2)))

upper_idx = Int64[]
for i = 1:m+1
    push!(upper_idx, (div(n,3)-1)*(m+1)+i)
    push!(upper_idx, (div(n,3)-1)*(m+1)+i + (m+1)*(n+1))
end
for i = 1:m
    push!(upper_idx, (div(n,3)-1)*m+i+2(m+1)*(n+1))
end

if is_training
    Ue = matread("U.mat")["U"]
    loss = sum((u_out[:, upper_idx] - Ue)^2)
    sess = Session(); init(sess)
    loss_ = BFGS!(sess, loss, 500)

    close("all")
    semilogy(loss_)
    xlabel("Iterations")
    ylabel("Loss")
    savefig("loss.png")
else
    sess = Session(); init(sess)
    U = run(sess, u_out)
    matwrite("U.mat", Dict("U"=> U[:, upper_idx]))
    visualize_displacement(U'|>Array, m, n, h, name="_tf")
    visualize_pressure(U'|>Array, m, n, h, name="_tf")

end


text:
using Test
using StaticArrays
using VariableTemplates

@testset "NTuple type" begin

    struct SingleType{FT} end
    struct NTupleType{FT} end
    Base.@kwdef struct Container{FT, N}
        ntuple_type::NTuple{N, NTupleType{FT}} =
            ntuple(i -> NTupleType{FT}(), N)
        single_type::SingleType{FT} = SingleType{FT}()
    end

    vars_state(m::NTupleType, FT) = @vars(a::FT, b::SVector{3, FT})
    vars_state(m::NTuple{N, NTupleType}, FT) where {N} =
        Tuple{ntuple(i -> vars_state(m[i], FT), N)...}
    vars_state(::SingleType, FT) = @vars(a::FT)

    function vars_state(m::Container, FT)
        @vars(
            single_type::vars_state(m.single_type, FT),
            ntuple_type::vars_state(m.ntuple_type, FT)
        )
    end

    FT = Float32
    N = 5
    container = Container{FT, N}()
    st = vars_state(container, FT)
    num_state = varsize(vars_state(container, FT))

    # Vars
    local_state = MArray{Tuple{num_state}, FT}(undef)
    fill!(local_state, 0)
    v = Vars{st}(local_state)
    v.single_type.a = 2
    @test v.single_type.a == 2
    v.ntuple_type[1].a = 2
    @test v.ntuple_type[1].a == 2
    v.ntuple_type[2].b = SVector(FT(1), FT(2), FT(3))
    @test v.ntuple_type[2].b == SVector(FT(1), FT(2), FT(3))
    @test v.ntuple_type[1].b == SVector(FT(0), FT(0), FT(0))
    v.ntuple_type[3].a = 3
    @test v.ntuple_type[3].a == 3
    v.ntuple_type[N].a = 3
    @test v.ntuple_type[N].a == 3

    # Grad
    local_state = MArray{Tuple{3, num_state}, FT}(undef)
    fill!(local_state, 0)
    v = Grad{st}(local_state)
    v.single_type.a = SVector(2, 2, 2)
    @test v.single_type.a == SVector(2, 2, 2)
    v.ntuple_type[1].a = SVector(1, 2, 3)
    @test v.ntuple_type[1].a == SVector(1, 2, 3)

    v.ntuple_type[2].b = SArray{Tuple{3, 3}, FT}([1, 2, 3, 4, 5, 6, 7, 8, 9])
    @test v.ntuple_type[2].b ==
          SArray{Tuple{3, 3}, FT}([1, 2, 3, 4, 5, 6, 7, 8, 9])
    @test v.ntuple_type[3].b ==
          SArray{Tuple{3, 3}, FT}([0, 0, 0, 0, 0, 0, 0, 0, 0])
    @test v.ntuple_type[1].b ==
          SArray{Tuple{3, 3}, FT}([0, 0, 0, 0, 0, 0, 0, 0, 0])

    # flattenednames
    fn = flattenednames(st)
    @test fn == [
        "single_type.a",
        "ntuple_type[1].a",
        "ntuple_type[1].b[1]",
        "ntuple_type[1].b[2]",
        "ntuple_type[1].b[3]",
        "ntuple_type[2].a",
        "ntuple_type[2].b[1]",
        "ntuple_type[2].b[2]",
        "ntuple_type[2].b[3]",
        "ntuple_type[3].a",
        "ntuple_type[3].b[1]",
        "ntuple_type[3].b[2]",
        "ntuple_type[3].b[3]",
        "ntuple_type[4].a",
        "ntuple_type[4].b[1]",
        "ntuple_type[4].b[2]",
        "ntuple_type[4].b[3]",
        "ntuple_type[5].a",
        "ntuple_type[5].b[1]",
        "ntuple_type[5].b[2]",
        "ntuple_type[5].b[3]",
    ]

end


text:
abstract type
    Reserve <: Service
end

"""
ProportionalReserve(name::String,
        contributingdevices::Device,
        timeframe::Float64,
        requirement::Dict{Any,Dict{Int,TimeSeries.TimeArray}})

Data Structure for a proportional reserve product for system simulations.
The data structure can be called calling all the fields directly or using named fields.
name - description
contributingdevices - devices from which the product can be procured
timeframe - the relative saturation timeframe
requirement - the required quantity of the product

# Examples

```jldoctest

"""
struct ProportionalReserve <: Reserve
    name::String
    contributingdevices::Array{Device}
    timeframe::Float64
    requirement::TimeSeries.TimeArray
end

function ProportionalReserve(name::String,
                            contributingdevices::Array{G},
                            timeframe::Float64,
                            requirement::Float64,
                            loads::Array{T}) where {G <: Device, T <: ElectricLoad}

    totalload = zeros(0)
    for i in TimeSeries.timestamp(loads[1].scalingfactor)
        t = zeros(0)
        for load in loads
            push!(t,(load.maxactivepower*values(load.scalingfactor[i]))[1])
        end
        push!(totalload,sum(t))
    end
    requirement = TimeSeries.TimeArray(TimeSeries.timestamp(loads[1].scalingfactor),totalload*requirement)

    ProportionalReserve(name, contributingdevices, timeframe, requirement)
end

ProportionalReserve(;name = "init",
        contributingdevices = [ThermalDispatch()],
        timeframe = 0.0,
        requirement = 0.03,
        loads = [ PowerLoad()]) = ProportionalReserve(name, contributingdevices, timeframe, requirement, loads)

"""
StaticReserve(name::String,
        contributingdevices::Device,
        timeframe::Float64,
        requirement::Float64})

Data Structure for the procurement products for system simulations.
The data structure can be called calling all the fields directly or using named fields.
name - description
contributingdevices - devices from which the product can be procured
timeframe - the relative saturation timeframe
requirement - the required quantity of the product

# Examples

```jldoctest

"""
struct StaticReserve <: Reserve
    name::String
    contributingdevices::Array{Device}
    timeframe::Float64
    requirement::Float64

    function StaticReserve(name::String,
                                contributingdevices::Array{Q},
                                timeframe::Float64,
                                generators::Array{G}) where {Q <: Device, G <: TechThermal}

        requirement = maximum([gen.activepowerlimits[:max] for gen in generators])

        new(name, contributingdevices, timeframe, requirement)
    end
end

StaticReserve(;name = "init",
        contributingdevices = [ThermalDispatch()],
        timeframe = 0.0,
        generators = [TechThermal()]) = StaticReserve(name, contributingdevices, timeframe, generators)


text:
module BridgeLandmarks

using Bridge
using CSV
using DataFrames
using DelimitedFiles
using Distributions
using ForwardDiff
using GaussianDistributions
using JLD2
using LinearAlgebra
using Base.Iterators
using PDMats
using Parameters
using RecursiveArrayTools
using SparseArrays
using StaticArrays
using StatsFuns

dir() = joinpath(@__DIR__, "..")

const d = 2   # dimension of landmarks
const sk = 1  # entries to skip for likelihood evaluation
const itostrat =  true # include Ito-Stratonovich correction in AHS-model
const γconstant = true # if true, no parameter updating for γ

import Bridge: kernelr3!, R3!, target, auxiliary, constdiff, llikelihood, _b!, B!, σ!, b!

export Point, PointF, Unc, UncF, State, deepvec

export Landmarks, LandmarksAux, MarslandShardlow, MarslandShardlowAux, Pars_ms, Pars_ahs, FlatPrior, show_updates
export landmarksforward, itostrat, construct_nfs, lm_mcmc, gramkernel, landmarksmatching, template_estimation#, plotlandmarksmatching, plottemplate_estimation
export Lmplotbounds, extractcomp, tc
export d, sk, to

#plotlandmarkpositions = Ref{Any}((args...) -> nothing )

include("nstate.jl")
include("state.jl")

include("backwardsfiltering.jl")
include("guidedproposal.jl")
include("landmarksmatching.jl")
include("lmguid.jl")  # contains main routines for mcmc
include("models.jl")

include("obsinfo.jl") # set observation info
include("pars.jl")  # set tuning pars
include("patches.jl")
include("postprocessing.jl")

include("template_estimation.jl")
include("updatematching.jl")

# include("basic_outplots.jl")
# include("plotting.jl")

end # module


text:
import Base.show, Base.rand
import POMDPs:
        isterminal,
        pdf,
        observation,
        generate_s,
        generate_o,
        initialstate_distribution,
        discount,
        reward

import LPDM: isterminal, max_belief_clusters
using StatsFuns
using Random

# Typealias appropriately
const LD1State  = Float64
const LD1Action = Float64
const LD1Obs    = Float64
const LD1Belief = LPDMBelief

Base.show(io::IO, x::Float64) = print(io,"$(@sprintf("%.4f", x))")

abstract type AbstractLD1 <: POMDP{Float64, Float64, Float64} end

"""
1-Dimensional light-dark problem

    min_noise::Float64
        minimum standard deviation of observation noise

    min_noise_loc::Float64
        x location of minimum observation noise

    Q::Matrix{Float64}
    R::Matrix{Float64}
        cost function matrices: reward = -s'Qs -a'Ra

    term_radius::Float64
        problem will terminate if the state is within this radius of zero

    init_dist::Any
        initial state distribution

    discount::Float64
        discount factor
"""
@with_kw mutable struct LightDark1D <: LightDarkLPDM.AbstractLD1
    min_noise::Float64      = 0.0
    min_noise_loc::Float64  = 5.0
    Q::Float64              = 0.5
    R::Float64              = 0.5
    term_radius::Float64    = 1e-5
    # init_dist::Any          = Normal(2.0, 0.5)
    #init_dist::Any          = Normal(-3.0, 0.5)
    discount::Float64       = 1.0
    count::Int              = 0
end

# POMDPs.jl API functions:
# POMDPs.generate_s(p::AbstractLD1, s::Float64, a::Float64, rng::AbstractRNG) = generate_s(p, s, a, rng, det = false)
POMDPs.generate_s(p::AbstractLD1, s::Float64, a::Float64, rng::AbstractRNG) = generate_s(p, s, a, rng, false)
POMDPs.generate_o(p::AbstractLD1, s::Float64, a::Float64, sp::Float64, rng::AbstractRNG) = generate_o(p, sp, rng)
POMDPs.observation(p::AbstractLD1, sp::Float64) = Normal1D(sp,obs_std(p,sp))
POMDPs.observation(p::AbstractLD1, a::Float64, sp::Float64) = observation(p, sp)
POMDPs.observation(p::AbstractLD1, s::Float64, a::Float64, sp::Float64) = observation(p,sp)

function generate_o(p::AbstractLD1, sp::Float64, rng::AbstractRNG)
    d = observation(p, sp)
    o = rand(rng, Distributions.Normal(d.mean, d.std))
    if p.obs_mode == :discrete
        o_disc = p.bin_centers[encode(p.lindisc,o)]
        return o_disc
    elseif p.obs_mode == :continuous
    #     o=clamp(o, -p.max_x, p.max_x) #TODO: EB, 02/11/2020: disable clamping for now, see if it stops clustering outliers into wrong clusters
        return o
    else
        error("Invalid obs_mode = $(p.obs_mode)")
    end
    # return obs_index(p,o_disc) # return a single combined obs index
end

# generate_o(p::AbstractLD1, sp::Float64, rng::AbstractRNG) = rand(rng, observation(p, sp))

# POMDPs.reward(p::AbstractLD1, s::Float64, a::Float64)              = -1.0

POMDPs.reward(p::AbstractLD1, s::Float64, a::Float64) =
                        (p.reward_mode == :quadratic) ? -(p.Q*s^2 + p.R*a^2) : -1

POMDPs.reward(p::AbstractLD1, s::Float64, a::Float64, sp::Float64) = POMDPs.reward(p,s,a)
POMDPs.discount(p::AbstractLD1) = p.discount

# Replaces the default call
function LPDM.isterminal(pomdp::AbstractLD1, particles::Vector{LPDMParticle{LD1State}})
    expected_state = 0.0
    weight_sum = 0.0
    # LPDM.normalize!(particles)
    for p in particles
        expected_state += p.state*p.weight
        weight_sum += p.weight
    end
    expected_state /= weight_sum
    return isterminal(pomdp,expected_state)
end

POMDPs.isterminal(p::AbstractLD1, s::Float64) = (abs(s) <= p.term_radius)

struct Normal1D
    mean::Float64
    std::Float64
end

# Base.rand(p::AbstractLD1, s::Float64, rng::LPDM.RNGVector) = rand(rng, Normal(s, std(p.init_dist)))

# function Base.rand(rng::LPDM.RNGVector,
#                    d::Normal)
#
#     # error("IN HERE!!!")
#     # a random number selected from normal distribution
#     r1 = norminvcdf(mean(d), std(d), rand(rng))
#     return r1
# end

# rand(rng::AbstractRNG, d::Normal) = d.mean + d.std*Float64(randn(rng, 2))
# POMDPs.pdf(d::Normal, o::Float64) = exp(-0.5*sum((o-d.mean).^2)/d.std^2)/(2*pi*d.std^2); println("observing!")
# function POMDPs.pdf(d::Normal, o::Float64) #DEBUG version
#     # error("In my pdf")
#     return exp(-0.5*sum((o-d.mean).^2)/d.std^2)/(2*pi*d.std^2)
# end
#
# mean(d::Normal) = d.mean
# mode(d::Normal) = d.mean
# Base.eltype(::Type{Normal}) = Float64

# chose this on 2/6/17 because I like the bowtie particle patterns it produces
# unclear which one was actually used in the paper
# Masters Thesis by Pas assumes the sqrt version
# obs_std(p::AbstractLD1, x::Float64) = sqrt(0.5*(p.min_noise_loc-x)^2 + p.min_noise)

# obs_std(p::AbstractLD1, x::Float64) = 0.5*(p.min_noise_loc-x)^2 + p.min_noise

# EB 10/20/18, implementing the version from Platt et al:
obs_std(p::AbstractLD1, x::Float64) = sqrt(0.5*(p.min_noise_loc-x)^2 + p.min_noise)

# NOTE: the whole point of defining Normal1D is to avoid the so-called type piracy when defining this POMDPs.pdf method
POMDPs.pdf(d::Normal1D, o::Float64) = Distributions.pdf(Distributions.Normal(d.mean, d.std),o)

# function POMDPs.pdf(d::Normal1D, o::Float64)
#     # pdf=Distributions.pdf(d,o)
#     println("In pdf: d=$d, o=$o")
#     println("pdf=$(Distributions.pdf(Distributions.Normal(d.mean, d.std),o))")
#     return nothing
# end

# function generate_s(p::AbstractLD1, s::Float64, a::Float64, rng::RNGVector; det::Bool = false)
function generate_s(p::AbstractLD1, s::Float64, a::Float64, rng::RNGVector, det::Bool = false)
    p.count += 1
    if p.action_std == 0.0 || det
        return s + a
    else
        return s + LPDM.rand(rng, Distributions.Normal(a, p.action_std))
    end
end

# function generate_s_det(p::AbstractLD1, s::Float64, a::Float64)
#     p.count += 1
#     return s+a
# end

# generate_o(p::AbstractLD1, sp::Float64, rng::AbstractRNG) = sp #DEBUG: trying no noise at all

# function generate_o(p::AbstractLD1, sp::Float64, rng::AbstractRNG)
#     println("sp=$sp")
#     return sp #DEBUG: trying no noise at all
# end

# Initial distribution corresponds to how the observations would look
function state_distribution(pomdp::AbstractLD1, s0::LD1State, config::LPDMConfig, rng::RNGVector)
    states = Vector{LPDMParticle{Float64}}();
    weight = 1/(config.n_particles^2) # weight of each individual particle
    particle = LPDMParticle{LD1State}(0.0, 1, weight)
    d = observation(pomdp,s0)
    for i = 1:config.n_particles^2 #TODO: Inefficient, possibly improve. Maybe too many particles
        particle = LPDMParticle{Float64}(LPDM.rand(rng, Distributions.Normal(d.mean, d.std)), i, weight)
        push!(states, particle)
    end
    # println("n states: $(length(states))")
    return states
end

function POMDPs.generate_sor(p::AbstractLD1, s::Float64, a::Float64, rng::RNGVector, det::Bool = false)
# function POMDPs.generate_sor(p::AbstractLD1, s::Float64, a::Float64, rng::RNGVector; det::Bool = false)

    s = generate_s(p, s, a, rng, det)
    o = generate_o(p, s, rng)
    r = reward(p, s, a)
    # println("o! $o")
    return s, o, r
end

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
abstract type Contact{T} <: Bound{T} end

### Constraints and derivatives
## Discrete-time position wrappers (for dynamics)
@inline g(contact::Contact, state::State, Δt) = g(contact, posargsnext(state, Δt)...)

## Position level constraints (for dynamics)
# @inline g(contact::Contact{T}) where {T} = szeros(T,6)
@inline g(contact::Contact, x::AbstractVector, q::UnitQuaternion) = contact.ainv3 * (x + vrotate(contact.p,q) - contact.offset)

## Derivatives NOT accounting for quaternion specialness
@inline function ∂g∂pos(contact::Contact, x::AbstractVector, q::UnitQuaternion)
    p = contact.p
    X = contact.ainv3
    Q = contact.ainv3 * (VLmat(q) * Lmat(UnitQuaternion(p)) * Tmat() + VRᵀmat(q) * Rmat(UnitQuaternion(p)))
    return X, Q
end

## Discrete-time position derivatives (for dynamics)
# Wrappers 1
@inline function ∂g∂ʳpos(contact::Contact, body::Body)
    return ∂g∂ʳpos(contact, body.state)
end

#Wrappers 2
@inline ∂g∂ʳpos(contact::Contact{T}) where T = szeros(T, 6)
∂g∂ʳpos(contact::Contact, state::State) = ∂g∂ʳpos(contact, posargsk(state)...)

#Derivatives accounting for quaternion specialness
@inline function ∂g∂ʳpos(contact::Contact, x::AbstractVector, q::UnitQuaternion)
    X, Q = ∂g∂pos(contact, x, q)
    Q = Q * LVᵀmat(q)

    return [X Q]
end

## Discrete-time velocity derivatives (for dynamics)
# Wrappers 1
@inline function ∂g∂ʳvel(contact::Contact, body::Body, Δt)
    return ∂g∂ʳvel(contact, body.state, Δt)
end

#Wrappers 2
@inline ∂g∂ʳvel(contact::Contact{T}) where {T} = szeros(T, 6)
∂g∂ʳvel(contact::Contact, state::State, Δt) = ∂g∂ʳvel(contact, posargsnext(state, Δt)..., fullargssol(state)..., Δt)

#Derivatives accounting for quaternion specialness
@inline function ∂g∂ʳvel(contact::Contact,x2::AbstractVector, q2::UnitQuaternion, x1::AbstractVector,v1::AbstractVector, q1::UnitQuaternion,w1::AbstractVector, Δt)
    X, Q = ∂g∂pos(contact, x2, q2)
    V = X * Δt
    Ω = Q * Lmat(q1) * derivωbar(w1, Δt) * Δt / 2
    return [V Ω]
end

##Schurf und SchurD
#Schurf
@inline function schurf(ineqc, contact::Contact, i, body::Body, μ, Δt)
    return schurf(contact, body.state, ineqc.γsol[2][i], ineqc.ssol[2][i], μ, Δt)[SA[1; 2; 3; 4; 5; 6]]
end
@inline function schurf(contact::Contact, state::State, γ, s, μ, Δt)
    φ = g(contact, posargsnext(state, Δt)...)
    return ∂g∂ʳpos(contact, state)' * (γ / s * φ - μ / s)
end

#SchurD
@inline function schurD(ineqc, contact::Contact, i, body::Body, Δt)
    return schurD(contact, body.state, ineqc.γsol[2][i], ineqc.ssol[2][i], Δt)
end
@inline function schurD(contact::Contact,state::State, γ, s, Δt)
    return ∂g∂ʳpos(contact, state)' * γ/s * ∂g∂ʳvel(contact, state, Δt)
end

## Additional force for friction
@inline additionalforce(bound::Contact, state::State) = additionalforce(bound, posargsk(state)...)
@inline additionalforce(bound::Contact, body::Body) = additionalforce(bound, body.state)


text:
using DelimitedFiles
using Plots

EVDdata = DelimitedFiles.readdlm("../Assets/wikipediaEVDdatesconverted.csv", ',')  # don't forget the delimiter!
epidays = EVDdata[:, 1]
allcases = EVDdata[:, 2]

gr() # specify the backend

plot(epidays, allcases,
     title = "West African EVD epidemic",
     xlabel = "Days since 22 March 2014",
     ylabel = "Total cases to date",
     marker = (:diamond, 8),
     line = (:path, "gray"),
     legend = false,
     grid = false)

savefig("../Assets/EVD.png")


text:
"""
Linear Functions
"""
module LinearTransforms

using SparseArrays, LinearAlgebra

export AbstractLinearTransform, LinearTransform, opMask,
    opMaskFactored,
    drandn, rrandn, dottest

import Base: *, -, +
import Base: show, size, adjoint, eltype, isreal, Matrix

include("abstractlinop.jl")
include("utilities.jl")
include("masking.jl")
include("operations.jl")

end # module


text:
using DistanceTransforms
using Test

@testset ExtendedTestSet "mean_hausdorff" begin
    @testset ExtendedTestSet "mean_hausdorff" begin
    x = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
    y = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]

    @test mean_hausdorff(x, y) ≈ 0
    end
end

@testset ExtendedTestSet "mean_dice" begin
    @testset ExtendedTestSet "mean_dice" begin
    x = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
    y = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]

    @test mean_dice(x, y) ≈ 1
    end
end


text:
"""
QQQQ
"""
function test_fit_example1(;sim_runs=10^2)
    Λ₄, λ45, λ54, Λ₅ = 15.0, 5.0, 7.0, 16.0
    μ41, μ42, μ43, μ51, μ52, μ53 = 4.0, 3.0, 3.0, 1.0, 7.0, 1.0
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]
    initial_dist = [0.5,0.5]

    maph = MAPHDist(initial_dist',T_example, T0_example)

    @show maph

    absorbtion_probs_of_true_dist = absorption_probs(maph)
    mean_of_true_dist = mean(maph)
    scvs_of_true_dist = scv(maph)

    times = zeros(sim_runs)
    absorbing_states = zeros(Int,sim_runs)

    @showprogress "Simulating data" for i in 1:sim_runs
        time, state = rand(maph, full_trace = false)
        times[i] = time
        absorbing_states[i] = state
    end

    estimated_dist = fit_maph(times, absorbing_states, 4, max_iter = 0)

    absorbtion_probs_of_estimated_dist = absorption_probs(estimated_dist)
    mean_of_estimated_dist = mean(estimated_dist)
    scvs_of_estimated_dist = scv(estimated_dist)

    @show absorbtion_probs_of_true_dist,absorbtion_probs_of_estimated_dist
    @show mean_of_estimated_dist, mean_of_true_dist
    @show scvs_of_estimated_dist, scvs_of_true_dist

    return true
end


text:
include(joinpath(JULIA_HOME, Base.DATAROOTDIR, "julia", "build_sysimg.jl"))

println("Building JuliaBox Julia system image...")
build_sysimg("/opt/julia_packages/jimg/stable/sys", "native", "/home/juser/jimg.jl", force=true)

max_stars_repo_name: shashi/JuliaBox

max_issues_repo_name: wsshin/JuliaBox

max_forks_repo_name: wsshin/JuliaBox


text:
# code inside VegaGraphs.jl
module VegaGraphs
using VegaLite
export sum_values
include("graph_functions.jl")
end

max_stars_repo_name: jpbrasile/VegaGraphs

max_issues_repo_name: jpbrasile/VegaGraphs

max_forks_repo_name: jpbrasile/VegaGraphs


text:
import SeisIO.Formats: formats, FmtVer, FormatDesc, HistVec
import SeisIO: BUF,
  InstrumentResponse,
  KW,
  SeedBlk,
  SeisIOBuf,
  TimeSpec,
  buf_to_double,
  buf_to_i16,
  buf_to_int,
  channel_match,
  check_for_gap!,
  checkbuf!,
  checkbuf_8!,
  endtime,
  fillx_u32_be!,
  fillx_u32_le!,
  fix_units,
  is_u8_digit,
  mk_t,
  mktime,
  stream_int,
  sμ,
  track_hdr!,
  trunc_x!,
  y2μs,
  μs


text:
using Bokeh
using Compat
if VERSION < v"0.4-"
    using Dates
else
    using Base.Dates
end

start = now()
x = map(d -> start + Dates.Hour(d), 1:(30 * 24))
println(typeof(x))
# x: Array{DateTime,1}
day = 1000 * 3600 * 24
@compat y = 0.5sin(map(Int, x) * pi / day) + cos(map(Int, x) * pi / (day * 7))
# note x needs to be an array of type DateTime or Date
plot(x, y)
showplot("dates_xy.html")

start = Date(1980, 1, 1)
x = map(d -> start + Dates.Day(d), 1:(365 * 30))
println(typeof(x))
# x: Array{Date,1}
@compat y = 0.5sin(map(Int, x) * pi / 365) + cos(map(Int, x) * pi / (365 * 7))

plot(y, x)
showplot("dates_yx.html")


text:
#=

Utility functions for simulations for tribe

Ignacio Quintero Mächler

August 29 2017

t(-_-t)

=#

"""
    simulate_tribe(X_initial::Float64,
                   nareas   ::Int64,
                   tree_file::String;
                   ωx       = 0.0,
                   σ²       = 0.5,
                   λ1       = 0.5,
                   λ0       = 0.2,
                   ω1       = 0.0,
                   ω0       = 0.0,
                   const_δt = 1e-4)

Simulate tribe model.
"""
function simulate_tribe(X_initial::Float64,
                        nareas   ::Int64,
                        tree_file::String;
                        ωx       = 0.0,
                        σ²       = 0.5,
                        λ1       = 0.5,
                        λ0       = 0.2,
                        ω1       = 0.0,
                        ω0       = 0.0,
                        const_δt = 1e-4)

  # bounds checks for parameters
  0.0 >= σ² && error("σ² has to be > 0.0")
  0.0 >= λ1 && error("λ1 has to be > 0.0")
  0.0 >= λ0 && error("λ0 has to be > 0.0")

  Y_initial =
    [rand() < λ1/(λ1 + λ0) ? 1 : 0 for i in Base.OneTo(nareas)]::Array{Int64,1}

  while iszero(sum(Y_initial))
    Y_initial =
      [rand() < λ1/(λ1 + λ0) ? 1 : 0 for i in Base.OneTo(nareas)]::Array{Int64,1}
  end

  tree, bts = read_tree(tree_file)

  br = branching_times(tree)

  # sort according to branching times
  brs = sortslices(br, dims = 1, by = x -> x[5], rev = true)

  # sort branching times
  sort!(bts, rev = true)

  # add present
  push!(bts, 0.0)

  # number of speciation events
  nbt = lastindex(bts) - 1

  # calculate speciation waiting times
  swt = Array{Float64,1}(undef,nbt)
  for i in Base.OneTo(nbt)
    swt[i] = bts[i] - bts[i+1]
  end

  # initial values
  Xt = fill(X_initial, 2)

  Y_initial = reshape(vec(Y_initial), 1, :)
  Yt = vcat(Y_initial, Y_initial)

  # start of alive
  alive = sortslices(br, dims = 1, by = x -> x[1])[1:2,2]

  nalive = lastindex(alive)

  # is ωx positive? (for lineage averages)
  isωxP = ωx >= 0.0

  rate = sqrt(const_δt*σ²)

  # loop through waiting times
  for j in Base.OneTo(nbt)

    nreps = reps_per_period(swt[j], const_δt)

    # simulate during the speciation waiting time
    nconst_sim!(Xt, Yt, nreps, const_δt, ωx, ω1, ω0, rate, λ1, λ0, isωxP)

    if j == nbt
      break
    end

    # which lineage speciates
    wsp = brs[j,2]

    # where to insert
    wti = findfirst(isequal(wsp), alive)

    # index to insert
    idx = sort(push!(collect(1:nalive), wti))

    # insert in Xt & Yt
    Xt = Xt[idx]
    Yt = Yt[idx,:]

    # update alive
    chs = brs[wsp .== brs[:,1],2]

    alive[wti] = chs[1]
    insert!(alive, wti+1, chs[2])

    nalive = length(alive)

  end

  tip_traits =
    Dict(convert(Int64, alive[i]) => Xt[i]   for i = Base.OneTo(nalive))
  tip_areas  =
    Dict(convert(Int64, alive[i]) => Yt[i,:] for i = Base.OneTo(nalive))

  pop!(bts)

  return tip_traits, tip_areas, tree, bts

end

"""
    reps_per_period(br_length::Float64, const_δt::Float64)

Estimate number of repetitions for each speciation waiting time.
"""
reps_per_period(br_length::Float64, const_δt::Float64) =
  round(Int64,cld(br_length,const_δt))

"""
    nconst_sim!(Xt   ::Array{Float64,1},
                Yt   ::Array{Int64,2},
                nreps::Int64,
                δt   ::Float64,
                rate ::Float64,
                ωx   ::Float64,
                λ1   ::Float64,
                λ0   ::Float64,
                ω1   ::Float64,
                ω0   ::Float64,
                isωxP::Bool)

Simulate biogeographic and trait evolution along a
speciation waiting time.
"""
function nconst_sim!(Xt   ::Array{Float64,1},
                     Yt   ::Array{Int64,2},
                     nreps::Int64,
                     δt   ::Float64,
                     ωx   ::Float64,
                     ω1   ::Float64,
                     ω0   ::Float64,
                     rate ::Float64,
                     λ1   ::Float64,
                     λ0   ::Float64,
                     isωxP::Bool)

  # n species and narea areas
  n, narea = size(Yt)
  Ytp      = similar(Yt)
  nch      = zeros(Int64, n)

  # allocate memory for lineage averages and differences
  δX = zeros(n, n)      # X pairwise differences
  δY = zeros(n, n)      # Y pairwise differences
  la = zeros(n)         # lineage averages
  ld = zeros(n, narea)  # area specific lineage rates

  for i in Base.OneTo(nreps)

    # estimate area and lineage averages and area occupancy
    δXY_la_ld!(δX, δY, la, ld, Xt, Yt, n, narea, isωxP)

    # trait step
    traitsam_1step!(Xt, la, δt, rate, ωx, n)

    # biogeographic step
    copyto!(Ytp, Yt)
    biogeosam_1step!(ω1, ω0, λ1, λ0, ld, Ytp, nch, δt, n, narea)

    while check_sam(Ytp, nch, n, narea)
      copyto!(Ytp, Yt)
      biogeosam_1step!(ω1, ω0, λ1, λ0, ld, Ytp, nch, δt, n, narea)
    end

    copyto!(Yt,Ytp)
  end

  return nothing
end

"""
    δXY_la_ld!(δX   ::Array{Float64,2},
               δY   ::Array{Float64,2},
               la   ::Array{Float64,1},
               ld   ::Array{Float64,2},
               Xt   ::Array{Float64,1},
               Yt   ::Array{Int64,2},
               n    ::Int64,
               narea::Int64)

Estimate area and lineage specific averages given sympatry configuration.
"""
function δXY_la_ld!(δX   ::Array{Float64,2},
                    δY   ::Array{Float64,2},
                    la   ::Array{Float64,1},
                    ld   ::Array{Float64,2},
                    Xt   ::Array{Float64,1},
                    Yt   ::Array{Int64,2},
                    n    ::Int64,
                    narea::Int64,
                    isωxP::Bool)
  @inbounds begin

    # estimate pairwise distances
    for j = Base.OneTo(n), i = Base.OneTo(n)
      i == j && continue
      # δX differences
      δX[i,j] = Xt[i] - Xt[j]
      # δY differences
      sk      = 0.0
      δY[i,j] = 0.0
      @simd for k = Base.OneTo(narea)
        if Yt[j,k] == 1
            sk += 1.0
            δY[i,j] += Float64(Yt[i,k])
        end
      end
      δY[i,j] /= sk
    end

    # estimate lineage averages
    la[:] .= 0.0
    if isωxP
      for j = Base.OneTo(n), i = Base.OneTo(n)
          i == j && continue
          la[i] += δX[j,i] * δY[j,i]
      end
    else
      for j = Base.OneTo(n), i = Base.OneTo(n)
          i == j && continue
          la[i] += sign(δX[j,i]) * δY[j,i] * exp(-abs(δX[j,i]))
      end
    end

    # estimate lineage sum of distances
    ld[:] .= NaN
    for k = Base.OneTo(narea), i = Base.OneTo(n)
      xmin = 1.0e20
      for j = Base.OneTo(n)
        j == i && continue
        if Yt[j,k] == 1
          x = abs(δX[j,i])
          iszero(x) && continue
          xmin = x < xmin ? x : xmin
        end
      end
      ld[i,k] = xmin == 1.0e20 ? 0.0 : xmin
    end
  end

  return nothing
end

"""
    traitsam_1step!(Xt  ::Array{Float64,1},
                    la  ::Array{Float64,1},
                    δt  ::Float64,
                    rate::Float64,
                    ωx  ::Float64,
                    n   ::Int64)

Sample one step for trait evolution history: `X(t + δt)`.
"""
function traitsam_1step!(Xt  ::Array{Float64,1},
                         la  ::Array{Float64,1},
                         δt  ::Float64,
                         rate::Float64,
                         ωx  ::Float64,
                         n   ::Int64)

  @inbounds begin

    for i in Base.OneTo(n)
      Xt[i] += Eδx(la[i], ωx, δt) + randn()*rate
    end

  end
end

"""
    biogeosam_1step!(ω1   ::Float64,
                     ω0   ::Float64,
                     λ1   ::Float64,
                     λ0   ::Float64,
                     ld   ::Array{Float64,2},
                     Ytp  ::Array{Int64,2},
                     nch  ::Array{Int64,1},
                     δt   ::Float64,
                     n    ::Int64,
                     narea::Int64)

Sample one step for biogeographic history: `Y(t + δt)`.
"""
function biogeosam_1step!(ω1   ::Float64,
                          ω0   ::Float64,
                          λ1   ::Float64,
                          λ0   ::Float64,
                          ld   ::Array{Float64,2},
                          Ytp  ::Array{Int64,2},
                          nch  ::Array{Int64,1},
                          δt   ::Float64,
                          n    ::Int64,
                          narea::Int64)

  @inbounds begin

    nch[:] .= 0
    for k = Base.OneTo(narea), i = Base.OneTo(n)
      if iszero(Ytp[i,k])
        if rand() < f_λ(λ1,ω1,ld[i,k])*δt
          setindex!(Ytp,1,i,k)
          nch[i] += 1
        end
      else
        if rand() < f_λ(λ0,ω0,ld[i,k])*δt
          setindex!(Ytp,0,i,k)
          nch[i] += 1
        end
      end
    end

  end

  return nothing
end

"""
    check_sam(Ytc::Array{Int64,2}, nch::Array{Int64,1}, n::Int64, k::Int64)

Returns false if biogeographic step consists of *only one* change
or if the species does *not* go globally extinct.
"""
function check_sam(Ytp   ::Array{Int64,2},
                   nch   ::Array{Int64,1},
                   n     ::Int64,
                   nareas::Int64)

  @inbounds begin

    # check lineage did not go extinct
    for i in Base.OneTo(n)
      s = 0
      @simd for k in Base.OneTo(nareas)
        s += Ytp[i,k]
      end
      if iszero(s)
        return true
      end
    end

    # check that, at most, only one event happened per lineage
    for i in nch
      if i > 1
        return true
      end
    end

    return false
  end
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
module Curves

import ..CryptoGroups: <|, tobits, validate, a, b, cofactor, order

include("elliptic_curves.jl")
include("ecpoint.jl")

end


text:
module Neural

using Random: AbstractRNG, GLOBAL_RNG
using ArgCheck, Flux

import Flux, Zygote, Tracker, Distributions, BSON

### Tracker extensions

trackerparams(m) = m |> Flux.params |> Tracker.Params

Flux.gradient(f, ps::Tracker.Params; kwargs...) = Tracker.gradient(f, ps; kwargs...)
Flux.data(xs::Tracker.Params) = Tracker.data(xs)
Flux.data(xs::Tracker.TrackedArray) = Tracker.data(xs)

# https://github.com/FluxML/Flux.jl/blob/bdeb9c6d584668c7cef1ce71caf659d611c86d65/src/optimise/train.jl#L9-L18
function Flux.Optimise.update!(opt, xs::Tracker.Params, Δs)
    for x in xs
        Δs[x] == nothing && continue
        Tracker.update!(x, -Flux.Optimise.apply!(opt, Tracker.data(x), Tracker.data(Δs[x])))
    end
end

function Tracker.gradient(f, xs::Tracker.Params; once=true)
    l = f()
    Tracker.losscheck(l)
    Tracker.@interrupts Tracker.back!(l; once=once)
    gs = Tracker.Grads()
    for x in xs
        gs[Tracker.tracker(x)] = Tracker.extract_grad!(x)
    end
    return gs
end

track_arr(x) = x isa AbstractArray ? Tracker.TrackedArray(x) : x
track(m) = Flux.fmap(track_arr, m)

# Avoid tracking some fields
function Flux.fmap(f::typeof(track_arr), m::BatchNorm; cache=IdDict())
    haskey(cache, m) && return cache[x]
    cache[m] = BatchNorm(m.λ, f(m.β), f(m.γ), m.μ, m.σ², m.ϵ, m.momentum)
end

export trackerparams, track

### Flux extensions

using Humanize: Humanize

nparams(m) = sum(prod.(size.(Flux.params(m))))
function nparams(T, m)
    @argcheck T isa AbstractString
    n = nparams(m)
    return ishuman ? Humanize.digitsep(n) : n
end

Zygote.@nograd nparams

export nparams

### Zygote extensions

function Zygote.pullback(f, p1::Flux.Params, p2::Flux.Params, prest::Flux.Params...)
    ps = (p1, p2, prest...)
    n = length(ps)
    cx = Zygote.Context()
    ys, back = Zygote._pullback(cx, f)
    @assert length(ys) == n
    caches = [copy(cx.cache) for i in 1:n]
    function back_i(Δ, i)
        cx.cache = caches[i]
        for p in ps[i]
          Zygote.cache(cx)[p] = nothing
        end
        back(Δ)
        Zygote.Grads(cx.cache)
    end
    ys, tuple((Δ -> back_i(Δ, i) for i in 1:n)...)
end

function Zygote.gradient(f, p1::Flux.Params, p2::Flux.Params, prest::Flux.Params...)
    ys, backs = Zygote.pullback(f, p1, p2, prest...)
    n = length(ys)
    makesens(i) = tuple(map(j -> i == j ? Zygote.sensitivity(ys[i]) : zero(ys[i]), 1:n)...)
    return tuple([backs[i](makesens(i)) for i in 1:n]...)
end

Zygote.@nograd Flux.gpu

# TODO: enable below if saving and loading params is buggy
# Flux.trainable(bn::BatchNorm) = (bn.β, bn.γ, bn.μ, bn.σ²)

###

abstract type AbstractNeuralModel end

include("trainable.jl")
export train!, saveparams, loadparams!

NeuralNet(Sin::Int, Sout::Int, args...) = DenseNet(Sin, Sout, args...)
NeuralNet(Sin::NTuple{3, Int}, Sout::Int, args...) = ConvNet(Sin, Sout, args...)
NeuralNet(Sin::Int, Sout::NTuple{3, Int}, args...) = ConvNet(Sin, Sout, args...)
include("architecture.jl")
export NeuralNet, DenseNet, ConvNet

export Neural

end # module


text:
# Autogenerated wrapper script for basiclu_jll for powerpc64le-linux-gnu
export libbasiclu

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("basiclu")
JLLWrappers.@declare_library_product(libbasiclu, "libbasiclu.so.2")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libbasiclu,
        "lib/libbasiclu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Interpolation routines for 4th order Weighted Essentially Non-Oscillatory (WENO)
schemes.
"""
module WENO4

export interpolate_weno4

"""
    function interpolate_weno4(
        xs::AbstractVector{<:Number},
        xp::AbstractVector{<:Number},
        fp::AbstractVector{<:Number};
        extrapolate=false
    )

Performs 1D interpolation using the fourth-order Weighted Essentially
Non-Oscillatory (WENO) scheme from [Janett et al (2019)](https://ui.adsabs.harvard.edu/abs/2019A%26A...624A.104J/abstract).
Based on https://github.com/Goobley/Weno4Interpolation by Chris Osborne.

# Arguments

- `xs::AbstractVector{<:Number}`: The x-coordinates at which to evaluate the
  interpolated values
- `xp::AbstractVector{<:Number}`: The x-coordinates of the data points. Must
  have at least four elements, be increasing, and have no NaNs.
- `fp::AbstractVector{<:Number}`: The y-coordinates of the data points. Must
  have same length as `xp`.

# Keywords

- `extrapolate::Bool`: whether or not to use extrapolation. If `true`,
  will perform quadratic extrapolation from the edge points. When
  `false`, it will not throw an error if `xs` is outside `xp`, but return
  the closest `fp` value.

"""
function interpolate_weno4(
    xs::AbstractVector{<:Number},
    xp::AbstractVector{<:Number},
    fp::AbstractVector{<:Number};
    extrapolate=false,
)
    Ngrid = length(xp)
    result = similar(xs)
    ε = 1.0f-6
    prevβ = -1
    β2 = zero(eltype(fp))
    β3 = zero(eltype(fp))
    @inbounds for (idx, x) in enumerate(xs)
        if x < xp[1]
            if !extrapolate
                result[idx] = fp[1]
                continue
            end
            # quadratic extrapolation from first three points
            i = 1
        elseif x > xp[end]
            if !extrapolate
                result[idx] = fp[end]
                continue
            end
            # quadratic extrapolation from last three points
            i = Ngrid
        else
            i = binary_search(xp, x)
        end
        if i == Ngrid
            i -= 1
        end
        # set stencil, pad with zeros when not relevant
        xi = xp[i]
        xip = xp[i+1]
        yi = fp[i]
        yip = fp[i+1]
        if i == 1
            xim = zero(eltype(xp))
            xipp = xp[i+2]
            yim = zero(eltype(fp))
            yipp = fp[i+2]
        elseif i == Ngrid - 1
            xim = xp[i-1]
            xipp = zero(eltype(xp))
            yim = fp[i-1]
            yipp = zero(eltype(fp))
        else
            xim = xp[i-1]
            xipp = xp[i+2]
            yim = fp[i-1]
            yipp = fp[i+2]
        end
        q2, q3 = weno4_q(x, xim, xi, xip, xipp, yim, yi, yip, yipp)
        if i == 1
            result[idx] = q3
        elseif i == Ngrid - 1
            result[idx] = q2
        else
            if i != prevβ  # reuse previously computed when possible
                β2, β3 = weno4_β(xim, xi, xip, xipp, yim, yi, yip, yipp)
                prevβ = i
            end
            γ2 = -(x - xipp) / (xipp - xim)
            γ3 = (x - xim) / (xipp - xim)
            α2 = γ2 / (ε + β2)
            α3 = γ3 / (ε + β3)
            ω2 = α2 / (α2 + α3)
            ω3 = α3 / (α2 + α3)
            result[idx] = ω2 * q2 + ω3 * q3
        end
    end
    return result
end

"""
    binary_search(a::AbstractVector{<:Number}, x)

Return the index of the last value in `a` less than or equal to `x`.
`a` is assumed to be sorted.
"""
function binary_search(a::AbstractVector{<:Number}, x)
    n = length(a)
    ileft = 2
    i = n - 1
    while ileft <= i
        mid = (ileft + i) >> 1
        if @inbounds x < a[mid]
            i = mid - 1
        else
            ileft = mid + 1
        end
    end
    return i
end

"""
Compute the smoothness indicators β2 and β3 for 4th order WENO interpolation.
Arguments are the x and y values for the 4-point stencil.
"""
function weno4_β(xim, xi, xip, xipp, yim, yi, yip, yipp)
    him = xi - xim
    hi = xip - xi
    hip = xipp - xip
    H = him + hi + hip
    # Derivatives
    yyim = -((2 * him + hi) * H + him * (him + hi)) / (him * (him + hi) * H) * yim
    yyim += ((him + hi) * H) / (him * hi * (hi + hip)) * yi
    yyim -= (him * H) / ((him + hi) * hi * hip) * yip
    yyim += (him * (him + hi)) / ((hi + hip) * hip * H) * yipp
    yyi = -(hi * (hi + hip)) / (him * (him + hi) * H) * yim
    yyi += (hi * (hi + hip) - him * (2 * hi + hip)) / (him * hi * (hi + hip)) * yi
    yyi += (him * (hi + hip)) / ((him + hi) * hi * hip) * yip
    yyi -= (him * hi) / ((hi + hip) * hip * H) * yipp
    yyip = (hi * hip) / (him * (him + hi) * H) * yim
    yyip -= (hip * (him + hi)) / (him * hi * (hi + hip)) * yi
    yyip += ((him + 2 * hi) * hip - (him + hi) * hi) / ((him + hi) * hi * hip) * yip
    yyip += ((him + hi) * hi) / ((hi + hip) * hip * H) * yipp
    yyipp = -((hi + hip) * hip) / (him * (him + hi) * H) * yim
    yyipp += (hip * H) / (him * hi * (hi + hip)) * yi
    yyipp -= ((hi + hip) * H) / ((him + hi) * hi * hip) * yip
    yyipp += ((2 * hip + hi) * H + hip * (hi + hip)) / ((hi + hip) * hip * H) * yipp
    # Smoothness indicators
    β2 = (hi + hip)^2 * (abs(yyip - yyi) / hi - abs(yyi - yyim) / him)^2
    β3 = (him + hi)^2 * (abs(yyipp - yyip) / hip - abs(yyip - yyi) / hi)^2
    return (β2, β3)
end

"""
Compute the q2 and q3 Lagrange interpolation factors for 4th order WENO interpolation.
Arguments are the value to interpolate plus the x and y values for the 4-point stencil.
"""
function weno4_q(x, xim, xi, xip, xipp, yim, yi, yip, yipp)
    him = xi - xim
    hi = xip - xi
    hip = xipp - xip
    q2 = yim * ((x - xi) * (x - xip)) / (him * (him + hi))
    q2 -= yi * ((x - xim) * (x - xip)) / (him * hi)
    q2 += yip * ((x - xim) * (x - xi)) / ((him + hi) * hi)
    q3 = yi * ((x - xip) * (x - xipp)) / (hi * (hi + hip))
    q3 -= yip * ((x - xi) * (x - xipp)) / (hi * hip)
    q3 += yipp * ((x - xi) * (x - xip)) / ((hi + hip) * hip)
    return (q2, q3)
end

end


text:
using Latexify
using LaTeXStrings
using Test

test_array = ["x/y * d" :x ; :( (t_sub_sub - x)^(2*p) ) 3//4 ]

#### Test the setting of default keyword arguments.

@test latexify("x * y") ==
raw"$x \cdot y$"

set_default(cdot = false)

@test latexify("x * y") ==
raw"$x y$"

@test get_default() == Dict{Symbol,Any}(:cdot => false)

set_default(cdot = true, transpose = true)

@test get_default() == Dict{Symbol,Any}(:cdot => true,:transpose => true)
@test get_default(:cdot) == true
@test get_default(:cdot, :transpose) == (true, true)
@test get_default([:cdot, :transpose]) == Bool[1, 1]

reset_default()
@test get_default() == Dict{Symbol,Any}()

@test latexify("x * y") ==
raw"$x \cdot y$"


text:
"""
Sampling parameters to use with a texture.
"""
Base.@kwdef struct Sampling
    magnification::Vk.Filter = Vk.FILTER_CUBIC_IMG
    minification::Vk.Filter = Vk.FILTER_CUBIC_IMG
    mipmap_mode::Vk.SamplerMipmapMode = Vk.SAMPLER_MIPMAP_MODE_LINEAR
    address_modes::NTuple{3,Vk.SamplerAddressMode} = ntuple(Returns(Vk.SAMPLER_ADDRESS_MODE_CLAMP_TO_BORDER), 3)
    mip_lod_bias::Float32 = 0.
    anistropy_enable::Bool = true
    max_anisotropy::Float32 = 16.
    compare_enable::Bool = false
    compare_op::Vk.CompareOp = Vk.COMPARE_OP_GREATER_OR_EQUAL
    lod_bounds::NTuple{2,Float32} = (0., 1000.)
    border_color::Vk.BorderColor = Vk.BORDER_COLOR_INT_OPAQUE_BLACK
    unnormalized_coordinates::Bool = false
end

function Vk.Sampler(device, sampling::Sampling)
    create_info = Vk.SamplerCreateInfo(
        sampling.magnification,
        sampling.minification,
        sampling.mipmap_mode,
        sampling.address_modes...,
        sampling.mip_lod_bias,
        sampling.anistropy_enable,
        sampling.max_anisotropy,
        sampling.compare_enable,
        sampling.compare_op,
        sampling.lod_bounds...,
        sampling.border_color,
        sampling.unnormalized_coordinates,
    )
    unwrap(Vk.create_sampler(device, create_info))
end

const DEFAULT_SAMPLING = Sampling()

"""
Texture identified by name with sampling parameters.

This texture is to be transformed into a texture index (to index into an array of sampled images or combined image-samplers depending on whether sampling parameters are provided) to be included as a material parameter in push constant data.
"""
struct Texture
    name::Symbol
    sampling::Union{Nothing,Sampling}
end

Texture(name::Symbol) = Texture(name, nothing)


text:
@testset "lq_approx.jl" begin

    @testset "QuantEcon LQ Lecture Example" begin
        # Set parameters
        α, β, ρ1, ρ2, σ = 10.0, 0.95, 0.9, 0.0, 1.0

        R = 1 / β
        A = [1.0 0.0 0.0;
               α  ρ1  ρ2;
            0.0 1.0 0.0]
        A12 = zeros(3,1)
        ALQ_l = hcat(A, A12)
        ALQ_r = [0 -R 0 R]
        ALQ = vcat(ALQ_l, ALQ_r)

        RLQ = [0.0 0.0 0.0 0.0;
               0.0 0.0 0.0 0.0;
               0.0 0.0 0.0 0.0;
               0.0 0.0 0.0 1e-9]

        QLQ = 1.0
        BLQ = [0.0; 0.0; 0.0; R]
        β_LQ = β

        f(s, x) = -([1, s...]' * RLQ * [1, s...] + x' * QLQ * x);
        g(s, x) = ALQ * [1, s...] + BLQ * x;

        s_star, x_star = [1.0, 10.9, 1.0], 0.0;

        f_star = f(s_star, x_star)
        Df_star = [0.0, 0.0, -2.0e-9, 0.0]
        D²f_star =  [0.0 0.0 0.0 0.0;
                     0.0 0.0 0.0 0.0;
                     0.0 0.0 -2.0e-9 0.0;
                     0.0 0.0 0.0 -2.0]
        g_star = g(s_star, x_star)[2:end]
        Dg_star = [0.9 0.0 0.0 0.0;
                   1.0 0.0 0.0 0.0;
                   -1.052631578947368 0.0 1.052631578947368 1.052631578947368]

       lq = approx_lq(s_star, x_star, f_star, Df_star, D²f_star, g_star,
                      Dg_star, β)

       @test isapprox(lq.R, RLQ)
       @test isapprox(lq.Q, [QLQ])
       @test isapprox(lq.A, ALQ)
       @test isapprox(lq.B, BLQ)
    end

    @testset "One State and One Action Variable Example" begin
        A = [1.0 0.0;
             0.5 0.9];

        R = [5.0 0.0;
             0.0 0.3]

        N = [0.05 0.1]

        Q = [0.1];

        B = [0.0; 1.5];

        f(s, x) = -([1, s...]' * R * [1, s...] .+ x' * Q * x
                    .+ 2 * x' * N * [1, s...]);
        g(s, x) = A * [1, s...] + B * x;

        s_star = 5.0;
        x_star = 0.0;

        g_star = g(s_star, x_star)[2];
        f_star = f(s_star, x_star)[1];

        Df_star = [-3.0;
                   -1.1];
        D²f_star = [-0.6 -0.2;
                    -0.2 -0.2]
        Dg_star = [0.9 1.5];

        lq = approx_lq(s_star, x_star, f_star, Df_star, D²f_star, g_star,
                       Dg_star, 1.0)

        @test isapprox(lq.R, R)
        @test isapprox(lq.Q, Q)
        @test isapprox(lq.A, A)
        @test isapprox(lq.B, B)
        @test isapprox(lq.N, N)
    end

    @testset "One State and Two Action Variables Example" begin
        A = [1.0 0.0;
               0.5 0.9];

        R = [5.0 0.0;
               0.0 0.3]

        N = [0.05 0.1;
                0.01 0.02]

        Q = [0.1 -0.05;
                -0.05 0.3];

        B = [0.0 0.0; 1.5 -0.9];

        f(s, x) = -([1, s...]' * R * [1, s...] .+ x' * Q * x
                    .+ 2 * x' * N * [1, s...]);
        g(s, x) = A * [1, s...] + B * x;

        x_star = [0.0, 0.0];
        s_star = 5.0;

        g_star = g(s_star, x_star)[2];
        f_star = f(s_star, x_star)[1];

        Df_star = [-3.0;
                   -1.1;
                   -0.22];
        D²f_star = [-0.6 -0.2 -0.04;
                    -0.2 -0.2 0.1;
                    -0.04 0.1 -0.6]
        Dg_star = [0.9 1.5 -0.9];

        lq = approx_lq(s_star, x_star, f_star, Df_star, D²f_star, g_star,
                       Dg_star, 1.0)

        @test isapprox(lq.R, R)
        @test isapprox(lq.Q, Q)
        @test isapprox(lq.A, A)
        @test isapprox(lq.B, B)
        @test isapprox(lq.N, N)

    end

    @testset "Three States and Two Action Variables Example" begin
        A = [1.0 0.0 0.0 0.0;
             -0.2 0.5 0.01 0.0;
             0.05 -0.5 0.8 0.3;
            -0.1 -0.2 0.3 0.6];

        R = [5.0 0.00 0.00 0.00;
             0.0 0.15 0.21 -0.05;
             0.0 0.21 0.5 -0.01;
             0.0 -0.05 -0.01 0.8]

        N = [-0.1 0.005 -0.15 -0.2;
             0.005 -0.15 -0.2 0.2]

        Q = [1.0 -0.5;
             -0.5 0.9];

        B = [0.0 0.0;
             0.3 0.5;
             0.0 0.2;
             0.1 0.3];

       f(s, x) = ([1, s...]' * R * [1, s...] + x' * Q * x
                  + 2 * x' * N * [1, s...]);
       g(s, x) = A * [1, s...] + B * x;

       s_star = [-0.7916666666666763, -19.583333333333226, -14.541666666666579];
       x_star = [0.0, 0.0];

       g_star = g(s_star, x_star)[2:end];  # Drop constant
       f_star = f(s_star, x_star);

       Df_star = [-7.008333333332452;
                  -19.62499999999752;
                  -22.79583333333047;
                  11.483749999998526;
                  2.264166666666381];

       D²f_star = [0.3 0.42 -0.1 0.01 -0.3;
                   0.42 1.0 -0.02 -0.3 -0.4;
                   -0.1 -0.02 1.6 -0.4 0.4;
                   0.01 -0.3 -0.4 2.0 -1.0;
                   -0.3 -0.4 0.4 -1.0 1.8];

       Dg_star = [0.5 0.01 0.0 0.3 0.5;
                  -0.5 0.8 0.3 0.0 0.2;
                  -0.2 0.3 0.6 0.1 0.3];

       lq = approx_lq(s_star, x_star, f_star, Df_star, D²f_star, g_star,
                      Dg_star, 1.0);

        @test isapprox(-lq.R, R)
        @test isapprox(-lq.Q, Q)
        @test isapprox(lq.A, A)
        @test isapprox(lq.B, B)
        @test isapprox(-lq.N, N)

    end

    @testset "Monetary Policy Example" begin  # See 9.4.6
        vlq =  1.0e+02 * [-1.516771891935916  # Obtained from Fackler's code
                          -1.421873224803633
                          -1.253495277719767
                          -1.054476825045721
                          -0.919716533913779
                          -0.849214404323941
                          -0.842970436276206
                          -0.900984629770574
                          -0.975359737401322
                          -1.023256984807047
                          -1.350214126639643
                          -1.257350153002016
                          -1.093041592907463
                          -0.900127220717386
                          -0.771471010069413
                          -0.707072960963543
                          -0.706933073399777
                          -0.771051347378114
                          -0.849495841998174
                          -0.899427782898555
                          -1.056517641369012
                          -0.967723054720698
                          -0.811553268604770
                          -0.630847057382630
                          -0.514399007702594
                          -0.462209119564662
                          -0.474277392968833
                          -0.550603827915108
                          -0.637187096513793
                          -0.691188424403486
                          -0.714520526767850
                          -0.631830020603504
                          -0.487868395455513
                          -0.325474425685280
                          -0.227338617457150
                          -0.193460970771124
                          -0.223841485627201
                          -0.318480162025382
                          -0.417271591592004
                          -0.477376999965666
                          -0.490780548132429
                          -0.414194122452051
                          -0.282440658271998
                          -0.138358929953671
                          -0.058535363177447
                          -0.042969957943327
                          -0.091662714251310
                          -0.204613632101397
                          -0.315613222635957
                          -0.381822711493587
                          -0.385297705462748
                          -0.314815360266339
                          -0.195270057054223
                          -0.069500570187802
                          -0.007989244863485
                          -0.010736081081270
                          -0.077741078841160
                          -0.209004238143153
                          -0.332211989645650
                          -0.404525558987249
                          -0.398071998758808
                          -0.333693734046368
                          -0.226356591802190
                          -0.118899346387675
                          -0.075700262515263
                          -0.096759340184955
                          -0.182076579396751
                          -0.331651980150650
                          -0.467067892621084
                          -0.545485542446652
                          -0.529103428020609
                          -0.470829243792138
                          -0.375700262515897
                          -0.286555258553288
                          -0.261668416132782
                          -0.301039735254381
                          -0.404669215918082
                          -0.572556858123887
                          -0.720180931562259
                          -0.804702661871796
                          -0.682156123064999
                          -0.627951325825841
                          -0.540961118528224
                          -0.464024275533553
                          -0.451345594080985
                          -0.502925074170520
                          -0.618762715802159
                          -0.798858518975902
                          -0.954621366392899
                          -1.043212483691748
                          -0.778391993248151
                          -0.726221889503649
                          -0.643301069195345
                          -0.572468306684642
                          -0.565893705716043
                          -0.623577266289547
                          -0.745518988405155
                          -0.931718872062866
                          -1.091551106469175
                          -1.182176917262681]

        bet = [0.8 0.5; 0.2 0.6]
        gamm = [-0.8, 0.0]
        Omega = [0.3 0.0; 0.0 1.0]
        s_target = [0., 1.]
        alpha = [0.9, 0.4]
        Sigma = 0.04 * Matrix(I, 2, 2)
        discount = 0.9;

        k = [3, 3];
        m = [10, 10];
        breaks = m - (k.-1);
        s_min = [-15, -10];
        s_max = [15, 10];
        basis = Basis(map(SplineParams, breaks, s_min, s_max, k)...);

        s_star = [0., 1.];  # steady-state states
        x_star = (s_star[1] - alpha[1] - dot(bet[1,:],  s_star)) / gamm[1];
        e_star = [0.0, 0.0];

        f(s, x) = -(s - s_target)' * Omega * (s - s_target) / 2
        g(s, x, e) = alpha + bet * s + gamm * x + e

        f_star = f(s_star, x_star)
        Df_star = [0.0, 0.0, 0.0]
        D²f_star = [-0.3 0.0 0.0;
                    0.0 -1.0 0.0;
                    0.0 0.0 0.0]
        g_star = g(s_star, x_star, e_star)
        Dg_star = [0.8 0.5 -0.8;
                   0.2 0.6 0.0]

        lq = approx_lq(s_star, x_star, f_star, Df_star, D²f_star, g_star,
                       Dg_star, discount);

        P, F, d = stationary_values(lq)

        states = [ones(100) nodes(basis)[1]]
        vlq_approx = -[states[i, 1:3]' * P * states[i, 1:3] for i ∈ 1:100]

        @test isapprox(vlq_approx, vlq)
    end

end


text:
include("dqns/dqns.jl")
include("policy_gradient/policy_gradient.jl")


text:
# Julia wrapper for header: /usr/include/openssl/hmac.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c None HMAC_CTX_init (Ptr{HMAC_CTX},) libcrypto
@c None HMAC_CTX_cleanup (Ptr{HMAC_CTX},) libcrypto
@c Int32 HMAC_Init (Ptr{HMAC_CTX}, Ptr{None}, Int32, Ptr{EVP_MD}) libcrypto
@c Int32 HMAC_Init_ex (Ptr{HMAC_CTX}, Ptr{None}, Int32, Ptr{EVP_MD}, Ptr{ENGINE}) libcrypto
@c Int32 HMAC_Update (Ptr{HMAC_CTX}, Ptr{Uint8}, size_t) libcrypto
@c Int32 HMAC_Final (Ptr{HMAC_CTX}, Ptr{Uint8}, Ptr{Uint32}) libcrypto
@c Ptr{Uint8} HMAC (Ptr{EVP_MD}, Ptr{None}, Int32, Ptr{Uint8}, size_t, Ptr{Uint8}, Ptr{Uint32}) libcrypto
@c Int32 HMAC_CTX_copy (Ptr{HMAC_CTX}, Ptr{HMAC_CTX}) libcrypto
@c None HMAC_CTX_set_flags (Ptr{HMAC_CTX}, Uint32) libcrypto