text:
#############
# Rotations #
#############

export Rx
"""
Rotation matrix, for a rotation about the x-axis.

Arguments:
- `angle::Real`: Counter-clockwise angle of rotation as viewed looking back along the postive direction of the rotation axis.
- `use_degrees:Bool`: If `true` interpret input as being in degrees.

Returns:
- `r::Array{<:Real, 2}`: Rotation matrix

References:
1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.27.
"""
function Rx(angle::Real ; use_degrees::Bool=false)

    if use_degrees
        angle *= pi/180.0
    end

    c = cos(angle)
    s = sin(angle)

    return [ +1.0  0.0  0.0;
              0.0 +c   +s;
              0.0 -s   +c]
end

export Ry
"""
Rotation matrix, for a rotation about the y-axis.

Arguments:
- `angle::Real`: Counter-clockwise angle of rotation as viewed looking back along the postive direction of the rotation axis.
- `use_degrees:Bool`: If `true` interpret input as being in degrees.

Returns:
- `r::Array{<:Real, 2}`: Rotation matrix

References:
1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.27.
"""
function Ry(angle::Real ; use_degrees::Bool=false)

    if use_degrees
        angle *= pi/180.0
    end

    c = cos(angle)
    s = sin(angle)

    return [ +c    0.0 -s;
              0.0 +1.0  0.0;
             +s    0.0 +c]
end

export Rz
"""
Rotation matrix, for a rotation about the z-axis.

Arguments:
- `angle::Real`: Counter-clockwise angle of rotation as viewed looking back along the postive direction of the rotation axis.
- `use_degrees:Bool`: If `true` interpret input as being in degrees.

Returns:
- `r::Array{<:Real, 2}`: Rotation matrix

References:
1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.27.
"""
function Rz(angle::Real ; use_degrees::Bool=false)

    if use_degrees
        angle *= pi/180.0
    end

    c = cos(angle)
    s = sin(angle)

    return [ +c    +s   0.0;
             -s    +c   0.0;
              0.0  0.0 +1.0]
end

###################
# Data Structures #
###################

export Quaternion
"""
The `Quaternion` type defines a _scalar-first_ quaternion for representation
of attitude.

Data members:
- `q0::Float64`: Scalar part of quaternion
- `q1::Float64`: First vector component of quaternion
- `q2::Float64`: Second vector component of quaternion
- `q3::Float64`: Third vector component of quaternion

References:
1. J. Diebel, _Representing attitude: Euler angles, unit quaternions, and rotation vectors._ Matrix 58(15-16) (2006).
"""
mutable struct Quaternion
    q0::Float64
    q1::Float64
    q2::Float64
    q3::Float64

    function Quaternion(q0::Real, q1::Real, q2::Real, q3::Real)
        n = sqrt(q0^2 + q1^2 + q2^2 + q3^2)

        return new(q0/n, q1/n, q2/n, q3/n)
    end
end

export EulerAngle
"""
The `EulerAngle` type provides a represenation of EulerAngles for storing attitude
information.

Valid sequences are: `121, 123, 131, 132, 212, 213, 231, 232, 312, 313, 321, 323`.

Data members:
- `seq::Integer`: Order of application of angles with respect to body axis.
- `phi::Float64`: First Euler angle
- `theta::Float64`: Second Euler angle
- `psi::Float64`: Third Euler angle

References:
1. J. Diebel, _Representing attitude: Euler angles, unit quaternions, and rotation vectors._ Matrix 58(15-16) (2006).
"""
mutable struct EulerAngle
    seq::Int
    phi::Float64
    theta::Float64
    psi::Float64

    function EulerAngle(seq::Integer, phi::Real, theta::Real, psi::Real)
        if !(seq in [121, 123, 131, 132, 212, 213, 231, 232, 312, 313, 321, 323])
            throw(ArgumentError("Invalid EulerAngle sequence: $seq"))
        end

        return new(seq, phi, theta, psi)
    end
end

export EulerAxis
"""
The `EulerAxis` type provides a representation of the Euler angle-and-axis attitude
representation.

Data members:
- `theta::Float64`: Angle of rotation
- `vec::Array{Float64, 1}`: Axis of rotation

References:
1. J. Diebel, _Representing attitude: Euler angles, unit quaternions, and rotation vectors._ Matrix 58(15-16) (2006).
"""
mutable struct EulerAxis
    angle::Float64
    axis::Array{Float64, 1}

    function EulerAxis(angle::Real, axis::Array{<:Real, 1})
        if length(axis) != 3
            throw(ArgumentError("Invalid array for EulerAxis initialization. Input size: $(size(axis)), Required size: (3,)"))
        end

        return new(angle, axis)
    end
end

##############
# Quaternion #
##############

# Quaternion Constructors
function Quaternion(vec::Array{<:Real, 1})
    if length(vec) != 4
        throw(ArgumentError("Invalid array for Quaternion initialization. Input length: $(length(vec)), Required length: 4"))
    end

    Quaternion(vec...)
end

function Quaternion(mat::Array{<:Real, 2})
    if size(mat) == (1,4)
        # Actually vector initialization. so it and return early
        return Quaternion(mat...)
    elseif size(mat) != (3,3)
        throw(ArgumentError("Invalid array for Quaternion initialization. Input size: $(size(mat)), Required size: (3,3)"))
    end

    temp = zeros(Float64, 4)
    temp[1] = 1 + mat[1, 1] + mat[2, 2] + mat[3, 3]
    temp[2] = 1 + mat[1, 1] - mat[2, 2] - mat[3, 3]
    temp[3] = 1 - mat[1, 1] + mat[2, 2] - mat[3, 3]
    temp[4] = 1 - mat[1, 1] - mat[2, 2] + mat[3, 3]

    # Get the maximum value and its index
    den, ind = findmax(temp)
    den      = sqrt(den) # Short-cut equivalence to reuse test information

    # Select optimal inverse mapping
    q0, q1, q2, q3 = 0.0, 0.0, 0.0, 0.0
    if ind == 1
        q0 = 0.5 * den
        q1 = 0.5 * (mat[2, 3] - mat[3, 2]) / den
        q2 = 0.5 * (mat[3, 1] - mat[1, 3]) / den
        q3 = 0.5 * (mat[1, 2] - mat[2, 1]) / den
    elseif ind == 2
        q0 = 0.5 * (mat[2, 3] - mat[3, 2]) / den
        q1 = 0.5 * den
        q2 = 0.5 * (mat[1, 2] + mat[2, 1]) / den
        q3 = 0.5 * (mat[3, 1] + mat[1, 3]) / den
    elseif ind == 3
        q0 = 0.5 * (mat[3, 1] - mat[1, 3]) / den
        q1 = 0.5 * (mat[1, 2] + mat[2, 1]) / den
        q2 = 0.5 * den
        q3 = 0.5 * (mat[2, 3] + mat[3, 2]) / den
    elseif ind == 4
        q0 = 0.5 * (mat[1, 2] - mat[2, 1]) / den
        q1 = 0.5 * (mat[3, 1] + mat[1, 3]) / den
        q2 = 0.5 * (mat[2, 3] + mat[3, 2]) / den
        q3 = 0.5 * den
    end

    return Quaternion(q0, q1, q2, q3)
end

function Quaternion(e::EulerAngle)
    # Extract Quaternion components
    q0, q1, q2, q3 = 0.0, 0.0, 0.0, 0.0

    # Compute sine and cosine values for angles

    # Get Quaternion components depending on Euler Angle sequence
    # reduce number of trig calls
    c1 = cos(e.phi/2.0)
    c2 = cos(e.theta/2.0)
    c3 = cos(e.psi/2.0)
    s1 = sin(e.phi/2.0)
    s2 = sin(e.theta/2.0)
    s3 = sin(e.psi/2.0)

    # Populate the quaternion
    if e.seq == 121
        q0 = c1*c2*c3 - s1*c2*s3
        q1 = c1*c2*s3 + c2*c3*s1
        q2 = c1*c3*s2 + s1*s2*s3
        q3 = c1*s2*s3 - s1*c3*s2

    elseif e.seq == 123
        q0 =  c1*c2*c3 + s1*s2*s3
        q1 = -c1*s2*s3 + c2*c3*s1
        q2 =  c1*c3*s2 + s1*c2*s3
        q3 =  c1*c2*s3 - s1*c3*s2

    elseif e.seq == 131
        q0 =  c1*c2*c3 - s1*c2*s3
        q1 =  c1*c2*s3 + c2*c3*s1
        q2 = -c1*s2*s3 + s1*c3*s2
        q3 =  c1*c3*s2 + s1*s2*s3

    elseif e.seq == 132
        q0 =  c1*c2*c3 - s1*s2*s3
        q1 =  c1*s2*s3 + c2*c3*s1
        q2 =  c1*c2*s3 + s1*c3*s2
        q3 =  c1*c3*s2 - s1*c2*s3

    elseif e.seq == 212
        q0 =  c1*c2*c3 - s1*c2*s3
        q1 =  c1*c3*s2 + s1*s2*s3
        q2 =  c1*c2*s3 + c2*c3*s1
        q3 = -c1*s2*s3 + s1*c3*s2

    elseif e.seq == 213
        q0 =  c1*c2*c3 - s1*s2*s3
        q1 =  c1*c3*s2 - s1*c2*s3
        q2 =  c1*s2*s3 + c2*c3*s1
        q3 =  c1*c2*s3 + s1*c3*s2

    elseif e.seq == 231
        q0 =  c1*c2*c3 + s1*s2*s3
        q1 =  c1*c2*s3 - s1*c3*s2
        q2 = -c1*s2*s3 + c2*c3*s1
        q3 =  c1*c3*s2 + s1*c2*s3

    elseif e.seq == 232
        q0 =  c1*c2*c3 - s1*c2*s3
        q1 =  c1*s2*s3 - s1*c3*s2
        q2 =  c1*c2*s3 + c2*c3*s1
        q3 =  c1*c3*s2 + s1*s2*s3

    elseif e.seq == 312
        q0 =  c1*c2*c3 + s1*s2*s3
        q1 =  c1*c3*s2 + s1*c2*s3
        q2 =  c1*c2*s3 - s1*c3*s2
        q3 = -c1*s2*s3 + c2*c3*s1

    elseif e.seq == 313
        q0 =  c1*c2*c3 - s1*c2*s3
        q1 =  c1*c3*s2 + s1*s2*s3
        q2 =  c1*s2*s3 - s1*c3*s2
        q3 =  c1*c2*s3 + c2*c3*s1

    elseif e.seq == 321
        q0 =  c1*c2*c3 - s1*s2*s3
        q1 =  c1*c2*s3 + s1*c3*s2
        q2 =  c1*c3*s2 - s1*c2*s3
        q3 =  c1*s2*s3 + c2*c3*s1

    elseif e.seq == 323
        q0 =  c1*c2*c3 - s1*c2*s3
        q1 = -c1*s2*s3 + s1*c3*s2
        q2 =  c1*c3*s2 + s1*s2*s3
        q3 =  c1*c2*s3 + c2*c3*s1
    else
        # Should get an invalid sequence, but it is possible if a user
        # Directly sets the sequence number
        throw(ArgumentError("Invalid EulerAngle sequence: $e.seq"))
    end

    return Quaternion(q0, q1, q2, q3)
end

function Quaternion(e::EulerAxis)
    # Extract Quaternion components
    q0 = cos(e.angle/2.0)
    q1 = e.axis[1]*sin(e.angle/2.0)
    q2 = e.axis[2]*sin(e.angle/2.0)
    q3 = e.axis[3]*sin(e.angle/2.0)

    return Quaternion(q0, q1, q2, q3)
end

#########################
# Quaternion Operations #
#########################

function Base.getindex(q::Quaternion, I::UnitRange{<:Integer})
    # Allocate vector once
    vec = as_vector(q)

    # Return selected index or range
    return [vec[i] for i in I]
end

function Base.getindex(q::Quaternion, I::Integer)
    if I == 1
        return q.q0
    elseif I == 2
        return q.q1
    elseif I == 3
        return q.q2
    elseif I == 4
        return q.q3
    else
        throw(BoundsError())
    end
end

Base.getindex(q::Quaternion, ::Colon) = [q.q0, q.q1, q.q2, q.q3]

# Return quaternion as a vector
export as_vector
"""
Return quaternion as a vector.

Equivalent to q[:]

Arguments:
- `q::Quaternion`: Quaternion

Returns:
- `vec::Array{Float64, 1}`: Quaternion as a (4,) vector
"""
function as_vector(q::Quaternion)
    return q[:]
end

# Return quaternion as a matrix
export as_matrix
"""
Return the rotation matrix representation of a Quaternion.

Arguments:
- `q::Quaternion`: Quaternion

Returns:
- `mat::Array{Float64, 2}`: Rotation Matrix on SO(3).
"""
function as_matrix(q::Quaternion)
    # initialize Empty Matrix
    mat = zeros(Float64, 3, 3)

    # Construct matrix from Quaternion
    mat[1, 1] = q.q0*q.q0 + q.q1*q.q1 - q.q2*q.q2 - q.q3*q.q3
    mat[1, 2] = 2*q.q1*q.q2 + 2*q.q0*q.q3
    mat[1, 3] = 2*q.q1*q.q3 - 2*q.q0*q.q2
    mat[2, 1] = 2*q.q1*q.q2 - 2*q.q0*q.q3
    mat[2, 2] = q.q0*q.q0 - q.q1*q.q1 + q.q2*q.q2 - q.q3*q.q3
    mat[2, 3] = 2*q.q2*q.q3 + 2*q.q0*q.q1
    mat[3, 1] = 2*q.q1*q.q3 + 2*q.q0*q.q2
    mat[3, 2] = 2*q.q2*q.q3 - 2*q.q0*q.q1
    mat[3, 3] = q.q0*q.q0 - q.q1*q.q1 - q.q2*q.q2 + q.q3*q.q3

    return mat
end

function Base.copy(q::Quaternion)
    return Quaternion(q.q0, q.q1, q.q2, q.q3)
end

function Base.deepcopy(q::Quaternion)
    return Quaternion(q.q0, q.q1, q.q2, q.q3)
end

"""
Compute the norm of a Quaternion.

Equivalent to `sqrt(q0^2 + q1^2 + q2^2 + q3^2)`

Arguments:
- `q::Quaternion`: Quaternion

Returns:
- `q_norm::Float64`: Norm of quaternion.
"""
function LinearAlgebra.norm(q::Quaternion)
    return sqrt(q.q0^2 + q.q1^2 + q.q2^2 + q.q3^2)
end

"""
Normalize a Quaternion in-place.

Equivalent to q=q/norm(q)

Arguments:
- `q::Quaternion`: Quaternion

Returns:
- `q_norm::Float64`: Norm of quaternion.
"""
function LinearAlgebra.normalize(q::Quaternion)
    # Get Quaternion norm
    q_norm = norm(q)

    # Normalize q in-place
    q.q0 = q.q0/q_norm
    q.q1 = q.q1/q_norm
    q.q2 = q.q2/q_norm
    q.q3 = q.q3/q_norm

    # Ensure return value is nothing
    nothing
end

"""
Get conjugate Quaternion.

Arguments:
- `q::Quaternion`: Input Quaternion

Returns:
- `q_conj::Quaternion`: Conjugate Quaternion of input
"""
function Base.conj(q::Quaternion)
    return Quaternion(q.q0, -q.q1, -q.q2, -q.q3)
end

"""
Get Quaternion inverse.

Arguments:
- `q::Quaternion`: Input Quaternion

Returns:
- `q_inv::Quaternion`: Inverse Quaternion of input
"""
function Base.inv(q::Quaternion)
    # Same as Quaternion conjugate since all quaternions are normalized to have
    # unit norm on construction
    return conj(q)
end

function Base.:-(q::Quaternion)
    return Quaternion(-q.q0, -q.q1, -q.q2, -q.q3)
end

function Base.:-(qa::Quaternion, qb::Quaternion)
    return [qa.q0 - qb.q0
            qa.q1 - qb.q1
            qa.q2 - qb.q2
            qa.q3 - qb.q3]
end

function Base.:+(qa::Quaternion, qb::Quaternion)
    return [qa.q0 + qb.q0
            qa.q1 + qb.q1
            qa.q2 + qb.q2
            qa.q3 + qb.q3]
end

function Base.:+(q::Quaternion, n::Real)
    return [q.q0 + n
            q.q1 + n
            q.q2 + n
            q.q3 + n]
end

function Base.:+(n::Real, q::Quaternion)
    return q+n
end

function Base.:*(qa::Quaternion, qb::Quaternion)
    # # Quaternion Multiplication
    # qcos = self.data[0]*other.data[0] - np.dot(self.data[1:4], other.data[1:4])
    # qvec = self.data[0]*other.data[1:4] + other.data[0]*self.data[1:4] + np.cross(self.data[1:4], other.data[1:4])

    # Quaternion Multiplication
    qcos = qa.q0*qb.q0 - dot(qa[2:4], qb[2:4])
    qvec = qa.q0*qb[2:4] + qb.q0*qa[2:4] + cross(qa[2:4], qb[2:4])

    return Quaternion(qcos, qvec...)
end

function Base.:*(q::Quaternion, n::Real)
    return q[:]*n
end

function Base.:*(n::Real, q::Quaternion)
    return q*n
end

export slerp
"""
Perform spherical linear interpolation (SLERP) on two quaternions. Interpolatles
from quaternion, `q1`, to quaternion, `q2`, at normalized interpolation time, `t`.

Interpolation time must be in the range `[0, 1]` a value of `0` will return `q1`,
while a value of `1` will return `q2`.

Arguments:
- `q1::Quaternion`: Starting Quaternion
- `q2::Quaternion`: Ending Quaternion
- `t::Real`: Normalized interpolation time. [0, 1]

Returns:
- `q:Quaternion`: Quaternion attitude interpolation from q1 toward q2 at time t.
"""
function slerp(q0::Quaternion, q1::Quaternion, t::Real)
    # Check Range on t
    if t < 0.0 || t > 1.0
        throw(ArgumentError("Invalid interpolation time $t. t must be in the range [0, 1]."))
    end

    # Extract vectors and normalize
    q0 = copy(q0)[:]
    q1 = copy(q1)[:]

    # Compute cosine of the angle between the two vectors
    dp = dot(q0, q1)

    # If the dot product is negative, the quaternions have opposite handed-ness
    # and slerp won't take the shortest path. Fix by reversing one quaternion.
    if dp < 0.0
        q1  = -q1
        dp = -dp
    end

    # If the inputs are too close we use linear interpolation instead
    if dp > 0.9995
        return Quaternion(q0 + (q1 - q0)*t)
    end

    theta0 = acos(dp) # Angle between input vectors
    theta  = theta0*t  # Angle between q0 and result quaternion

    s0 = cos(theta) - dp*sin(theta)/sin(theta0)
    s1 = sin(theta) / sin(theta0)

    return Quaternion((s0 * q0) + (s1 * q1))
end

##############
# EulerAngle #
##############

function EulerAngle(seq::Integer, vec::Array{<:Real, 1})
    if length(vec) != 3
        throw(ArgumentError("Invalid array for EulerAngle initialization. Input length: $(length(vec)), Required length: 3"))
    end

    EulerAngle(seq, vec...)
end

function EulerAngle(seq::Integer, mat::Array{<:Real, 2})
    if size(mat) != (3,3)
        throw(ArgumentError("Invalid array for Quaternion initialization. Input size: $(size(mat)), Required size: (3,3)"))
    end

    # Extract elements out of rotation matrix
    r11 = mat[1, 1]
    r12 = mat[1, 2]
    r13 = mat[1, 3]
    r21 = mat[2, 1]
    r22 = mat[2, 2]
    r23 = mat[2, 3]
    r31 = mat[3, 1]
    r32 = mat[3, 2]
    r33 = mat[3, 3]

    # Select euler angle sequence
    phi, theta, psi = 0.0, 0.0, 0.0
    if seq == 121
        phi   = atan(r21, r31)
        theta = acos(r11)
        psi   = atan(r12, -r13)
    elseif seq == 123
        phi   = atan(r23, r33)
        theta = -asin(r13)
        psi   = atan(r12, r11)
    elseif seq == 131
        phi   = atan(r31, -r21)
        theta = acos(r11)
        psi   = atan(r13, r12)
    elseif seq == 132
        phi   = atan(-r32, r22)
        theta = asin(r12)
        psi   = atan(-r13, r11)
    elseif seq == 212
        phi   = atan(r12, -r32)
        theta = acos(r22)
        psi   = atan(r21, r23)
    elseif seq == 213
        phi   = atan(-r13, r33)
        theta = asin(r23)
        psi   = atan(-r21, r22)
    elseif seq == 231
        phi   = atan(r31, r11)
        theta = -asin(r21)
        psi   = atan(r23, r22)
    elseif seq == 232
        phi   = atan(r32, r12)
        theta = acos(r22)
        psi   = atan(r23, -r21)
    elseif seq == 312
        phi   = atan(r12, r22)
        theta = -asin(r32)
        psi   = atan(r31, r33)
    elseif seq == 313
        phi   = atan(r13, r23)
        theta = acos(r33)
        psi   = atan(r31, -r32)
    elseif seq == 321
        phi   = atan(-r21, r11)
        theta = asin(r31)
        psi   = atan(-r32, r33)
    elseif seq == 323
        phi   = atan(r23, -r13)
        theta = acos(r33)
        psi   = atan(r32, r31)
    else
        throw(ArgumentError("Invalid EulerAngle sequence: $seq"))
    end

    return EulerAngle(seq, phi, theta, psi)
end

function EulerAngle(seq::Integer, q::Quaternion)
    # Construct angle from Quaternion by going through a rotation matrix
    return EulerAngle(seq::Integer, as_matrix(q))
end

function EulerAngle(seq::Integer, ea::EulerAxis)
    # Construct angle from EulerAxis by going through a rotation matrix
    return EulerAngle(seq::Integer, as_matrix(ea))
end

# Access Operators
function Base.getindex(e::EulerAngle, I::UnitRange{<:Integer})
    # Allocate vector once
    vec = as_vector(e)

    # Return selected index or range
    return [vec[i] for i in I]
end

function Base.getindex(e::EulerAngle, I::Integer)
    if I == 1
        return e.phi
    elseif I == 2
        return e.theta
    elseif I == 3
        return e.psi
    else
        throw(BoundsError())
    end
end

Base.getindex(e::EulerAngle, ::Colon) = [e.phi, e.theta, e.psi]

"""
Return Euler angles as a vector.

Equivalent to: `[e.phi, e.theta, e.psi]` for `EulerAngle` `e`

Arguments:
- `e::EulerAngle` Euler Angle

Returns:
- `evec::Array{Float64, 1}` Euler angles components in vector form.
"""
function as_vector(e::EulerAngle)
    return [e.phi, e.theta, e.psi]
end

function as_matrix(e::EulerAngle)
    # Get EulerAngle as matrix by going through Quaternions
    return as_matrix(Quaternion(e))
end

function Base.copy(e::EulerAngle)
    return EulerAngle(e.seq, e.phi, e.theta, e.psi)
end

function Base.deepcopy(e::EulerAngle)
    return EulerAngle(e.seq, e.phi, e.theta, e.psi)
end

#############
# EulerAxis #
#############

function EulerAxis(angle::Real, v1::Real, v2::Real, v3::Real)
    return EulerAxis(angle, [v1, v2, v3])
end

function EulerAxis(vec::Array{<:Real, 1})
    if length(vec) != 4
        throw(ArgumentError("Invalid array for EulerAxis initialization. Input size: $(size(vec)), Required size: (4,)"))
    end

    return EulerAxis(vec[1], [vec[2], vec[3], vec[4]])
end

function EulerAxis(mat::Array{<:Real, 2})
    if size(mat) != (3,3)
        throw(ArgumentError("Invalid array for EulerAxis initialization. Input size: $(size(mat)), Required size: (3,3)"))
    end

    return EulerAxis(Quaternion(mat))
end

function EulerAxis(q::Quaternion)
    # Extract quaternion vector and normalize
    qv = as_vector(q)
    q  = qv/norm(qv)

    # Ensure first element is positive
    if q[1] < 0
        q = -q
    end

    # Compute Euler Angle
    angle     = 2*acos(q[1])
    qvec_norm = norm(q[2:4])
    vec       = [0.0, 0.0, 0.0]

    if qvec_norm > 1e-15
        vec = q[2:4]/qvec_norm
    end

    return EulerAxis(angle, vec)
end

function EulerAxis(e::EulerAngle)
    # If input is an EulerAngle first compute a quaternion then get the
    # EulerAxis form
    EulerAxis(Quaternion(e))
end

# Access operators

function Base.getindex(e::EulerAxis, I::UnitRange{<:Integer})
    # Allocate vector once
    vec = as_vector(e)

    # Return selected index or range
    return [vec[i] for i in I]
end

function Base.getindex(e::EulerAxis, I::Integer)
    if I == 1
        return e.angle
    elseif I == 2
        return e.axis[1]
    elseif I == 3
        return e.axis[2]
    elseif I == 4
        return e.axis[3]
    else
        throw(BoundsError())
    end
end

Base.getindex(e::EulerAxis, ::Colon) = [e.angle, e.axis[1], e.axis[2], e.axis[3]]

function as_vector(e::EulerAxis)
    return e[:]
end

function as_matrix(e::EulerAxis)
    # Get matrix form from Quaternion for ease
    return as_matrix(Quaternion(e))
end

function Base.copy(e::EulerAxis)
    return EulerAxis(e.angle, e.axis)
end

function Base.deepcopy(e::EulerAxis)
    return EulerAxis(e.angle, e.axis)
end

####################
# Type Conversions #
####################


text:
function optionsResource(ctx, subscription_id::String, resource_group_name::String, uri::String, accountkey::Union{Nothing,AccountKey}=nothing, headers::Dict{Symbol,String}=Dict{Symbol,String}(); kwargs...)
    resource_op(ctx, "OPTIONS", subscription_id, resource_group_name, uri, accountkey, headers; kwargs...)
end

function putResource(ctx, subscription_id::String, resource_group_name::String, uri::String, accountkey::Union{Nothing,AccountKey}=nothing, headers::Dict{Symbol,String}=Dict{Symbol,String}(), data::Union{AbstractString,Vector{UInt8},IO,Nothing}=nothing; kwargs...)
    resource_op(ctx, "PUT", subscription_id, resource_group_name, uri, accountkey, headers, data; kwargs...)
end

function getResource(ctx, subscription_id::String, resource_group_name::String, uri::String, accountkey::Union{Nothing,AccountKey}=nothing, headers::Dict{Symbol,String}=Dict{Symbol,String}(); kwargs...)
    resource_op(ctx, "GET", subscription_id, resource_group_name, uri, accountkey, headers; kwargs...)
end

function deleteResource(ctx, subscription_id::String, resource_group_name::String, uri::String, accountkey::Union{Nothing,AccountKey}=nothing, headers::Dict{Symbol,String}=Dict{Symbol,String}(); kwargs...)
    resource_op(ctx, "DELETE", subscription_id, resource_group_name, uri, accountkey, headers; kwargs...)
end

function resource_op(ctx, op::String, subscription_id::String, resource_group_name::String, uri::String, accountkey::Union{Nothing,AccountKey}, headers::Dict{Symbol,String}, data::Union{AbstractString,Vector{UInt8},IO,Nothing}=nothing; kwargs...)
    if accountkey === nothing
        accountkey = extract_account_and_key(ctx, subscription_id, resource_group_name, uri)
    end
    req = Azure.REST.ServiceRequest(accountkey.account, op, uri, data; headers...)
    Azure.REST.execute(req, accountkey.key; kwargs...)
end

function extract_account_and_key(ctx, subscription_id::String, resource_group_name::String, uri::String)
    # find the storage account
    account = String(split(split(uri, "/")[3], '.'; limit=2)[1])

    # get keys for it
    allkeys = storageAccountsListKeys(Azure.api(ctx, StorageAccountsApi), resource_group_name, account, apiver(StorageAccountsApi), subscription_id)
    key = (allkeys.keys[1]).value

    AccountKey(account, key)
end

function append_params(uristr::String; kwargs...)
    haskwargs = false
    for (k,v) in kwargs
        if v !== nothing
            haskwargs = true
            break
        end
    end

    if haskwargs
        uri = URIs.parse_uri(uristr)
        query = URIs.queryparams(uri)
        for (k,v) in kwargs
            (v !== nothing) && (query[string(k)] = string(v))
        end
        uri = merge(uri; query=query)
        uristr = string(uri)
    end

    uristr
end

ext: jl


text:
import IRTools
import IRTools: IR, @dynamo, self, insertafter!

# Some abbreviations used in this file
# * sid : SSA ID, ID of a variable in SSA form of IR
# * tid : Tape ID, ID of a variable on a Tape
# * fn : function being called
# * args : arguments to a function
# * fargs : array of [fn, args...]
# * farg_defs : SSA definitions of fargs, i.e. IRTools.Variable or objects
# * res / ret - result or return value

################################################################################
#                                Frame                                         #
################################################################################

"""Frame of a call stack"""
mutable struct Frame
    # SSA var ID to Tape var ID. Note that SSA ID refers to the _original_ IR,
    # not the transformed one. In fact, any unique name of SSA instructions
    # would fit, using SSA IDs is just convenient
    ssa2tape::Dict{Int, Int}
    # result ID - tape ID corresponding to latest return value
    # from this call frame
    resultid::Int
end

function Base.show(io::IO, fr::Frame)
    map_str = join(["$k=>$v" for (k, v) in fr.ssa2tape], ",")
    print(io, "Frame($map_str, $(fr.resultid))")
end

################################################################################
#                        IRTracer (defs and utils)                             #
################################################################################

mutable struct IRTracer
    primitives::Set{Any}
    tape::Tape
    frames::Vector{Frame}
end

function IRTracer(;primitives=PRIMITIVES)
    tape = Tape()
    return IRTracer(primitives, tape, [])
end

Base.show(io::IO, t::IRTracer) = print(io, "IRTracer($(length(t.tape)))")

promote_const_value(x::QuoteNode) = x.value
promote_const_value(x::GlobalRef) = getproperty(x.mod, x.name)
promote_const_value(x) = x

function ssa_args_to_tape_vars!(t::IRTracer, arg_defs::Union{Vector, Tuple})
    result = Vector{Int}(undef, length(arg_defs))
    for (i, arg) in enumerate(arg_defs)
        if arg isa IRTools.Variable
            result[i] = t.frames[end].ssa2tape[arg.id]
        else
            val = promote_const_value(arg)
            arg_var = record!(t.tape, Constant, val)
            result[i] = arg_var
        end
    end
    return result
end

"""Push a new call frame to tracer, setting function params accordingly"""
function push_frame!(t::IRTracer, arg_defs...)
    tape_ids = ssa_args_to_tape_vars!(t, arg_defs)
    frame = Frame(Dict(i => tape_id for (i, tape_id) in enumerate(tape_ids)), -1)
    push!(t.frames, frame)
end

"""Pop call frame from tracer"""
function pop_frame!(t::IRTracer, res_sid::Int)
    frame = pop!(t.frames)
    # create mapping from the current SSA ID to the last instruction on the tape
    t.frames[end].ssa2tape[res_sid] =
        (frame.resultid == -1 ? length(t.tape) : frame.resultid)
end

"""Set target branch parameters to variables corresponding to SSA args"""
function set_branch_params!(t::IRTracer, ssa_args, target_params)
    tape_vars = ssa_args_to_tape_vars!(t, ssa_args)
    ssa2tape = t.frames[end].ssa2tape
    for (v, p) in zip(tape_vars, target_params)
        ssa2tape[p] = v
    end
end

"""Set return variable for the current frame"""
function set_return!(t::IRTracer, arg_sid_ref)
    # global STATE = (t, arg_sid_ref)
    tape_var = ssa_args_to_tape_vars!(t, [arg_sid_ref[]])[1]
    t.frames[end].resultid = tape_var
end

################################################################################
#                        IRTracer (body) + irtrace()                           #
################################################################################

function rewrite_special_cases!(ir::IR)
    for (v, st) in ir
        if Meta.isexpr(st.expr, :new)
            ir[v] = Expr(:call, __new__, st.expr.args...)
        end
    end
end

function record_or_recurse!(t::IRTracer, res_sid::Int, farg_defs, fargs...)
    fn, args = fargs[1], fargs[2:end]
    if fn in t.primitives || (fn isa Type && fn <: NamedTuple)
        res = fn(args...)
        tape_ids = ssa_args_to_tape_vars!(t, farg_defs[2:end])
        # record corresponding op to the tape
        res_tid = record!(t.tape, Call, res, fn, tape_ids)
        # update mapping from SSA var to tape var
        # note that in functions with loops this mapping may change over time
        t.frames[end].ssa2tape[res_sid] = res_tid
    else
        push_frame!(t, farg_defs...)
        res = t(fn, args...)
        pop_frame!(t, res_sid)
    end
    return res
end

function record_const!(t::IRTracer, res_sid, val)
    val = val isa QuoteNode ? val.value : val
    res_tid = record!(t.tape, Constant, val)
    t.frames[end].ssa2tape[res_sid] = res_tid
    return val
end

function trace_branches!(ir::IR)
    # if a block ends with a branch, we map its parameters to tape IDs
    # which currently correspond to argument SSA IDs
    for block in IRTools.blocks(ir)
        for branch in IRTools.branches(block)
            if IRTools.isreturn(branch)
                ret_v = branch.args[1]
                push!(ir, Expr(:call, set_return!, self, Ref(ret_v)))
            else
                ssa_args = branch.args
                target_params = [v.id for v in ir.blocks[branch.block].args]
                push!(block, Expr(:call, set_branch_params!, self, ssa_args, target_params))
            end
        end
    end
end

@dynamo function (t::IRTracer)(fargs...)
    ir = IR(fargs...)
    ir == nothing && return   # intrinsic functions
    rewrite_special_cases!(ir)
    for (v, st) in ir
        ex = st.expr
        # note the difference:
        # * `ex.args` is an array and thus will be passed to a function as is,
        # including definitions of IRTools.Variable;
        # * `ex.args...` is top-level to this expression and thus all Variable's
        # will be replaced with actual values during runtime
        if Meta.isexpr(ex, :call)
            ir[v] = Expr(:call, record_or_recurse!, self, v.id, ex.args, ex.args...)
        else
            # elseif ex isa GlobalRef
            ir[v] = Expr(:call, record_const!, self, v.id, ex)
        end
    end
    trace_branches!(ir)
    return ir
end

function irtrace(f, args...; primitives=PRIMITIVES, optimize=true)
    t = IRTracer(; primitives=primitives)
    for arg in args
        record!(t.tape, Input, arg)
    end
    push!(t.frames, Frame(Dict(i + 1 => i for i in 1:length(args)), -1))
    val = t(f, args...)
    t.tape.resultid = t.frames[1].resultid
    tape = t.tape
    if optimize
        tape = simplify(tape)
    end
    return val, tape
end


text:
using Test, LinearMaps, LinearAlgebra

@testset "function maps" begin
    N = 100
    function myft(v::AbstractVector)
        # not so fast fourier transform
        N = length(v)
        w = zeros(complex(eltype(v)), N)
        for k = 1:N
            kappa = (2*(k-1)/N)*pi
            for n = 1:N
                w[k] += v[n]*exp(kappa*(n-1)*im)
            end
        end
        return w
    end
    MyFT = @inferred LinearMap{ComplexF64}(myft, N) / sqrt(N)
    U = Matrix(MyFT) # will be a unitary matrix
    @test @inferred U'U ≈ Matrix{eltype(U)}(I, N, N)

    CS = @inferred LinearMap(cumsum, 2)
    @test size(CS) == (2, 2)
    @test @inferred !issymmetric(CS)
    @test @inferred !ishermitian(CS)
    @test @inferred !isposdef(CS)
    @test @inferred !(LinearMaps.ismutating(CS))
    @test @inferred Matrix(CS) == [1. 0.; 1. 1.]
    @test @inferred Array(CS) == [1. 0.; 1. 1.]
    CS = @inferred LinearMap(cumsum, 10; ismutating=false)
    v = rand(10)
    cv = cumsum(v)
    @test CS * v == cv
    @test *(CS, v) == cv
    @test_throws ErrorException CS' * v
    CS = @inferred LinearMap(cumsum, x -> reverse(cumsum(reverse(x))), 10; ismutating=false)
    cv = cumsum(v)
    @test @inferred CS * v == cv
    @test @inferred *(CS, v) == cv
    @test @inferred CS' * v == reverse!(cumsum(reverse(v)))
    @test @inferred mul!(similar(v), transpose(CS), v) == reverse!(cumsum(reverse(v)))

    CS! = @inferred LinearMap(cumsum!, 10; ismutating=true)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test_throws ErrorException CS!'v
    @test_throws ErrorException transpose(CS!) * v

    CS! = @inferred LinearMap{ComplexF64}(cumsum!, 10; ismutating=true)
    v = rand(ComplexF64, 10)
    cv = cumsum(v)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test_throws ErrorException CS!'v
    @test_throws ErrorException adjoint(CS!) * v
    @test_throws ErrorException mul!(similar(v), CS!', v)
    @test_throws ErrorException mul!(similar(v), transpose(CS!), v)
    CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10; ismutating=true)
    @inferred adjoint(CS!)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test @inferred CS' * v == reverse!(cumsum(reverse(v)))
    @test @inferred mul!(similar(v), transpose(CS), v) == reverse!(cumsum(reverse(v)))
    @test @inferred mul!(similar(v), adjoint(CS), v) == reverse!(cumsum(reverse(v)))

    # Test fallback methods:
    L = @inferred LinearMap(x -> x, x -> x, 10)
    v = randn(10)
    @test @inferred (2 * L)' * v ≈ 2 * v
    @test @inferred transpose(2 * L) * v ≈ 2 * v
    L = @inferred LinearMap{ComplexF64}(x -> x, x -> x, 10)
    v = rand(ComplexF64, 10)
    w = similar(v)
    @test @inferred (2 * L)' * v ≈ 2 * v
    @test @inferred transpose(2 * L) * v ≈ 2 * v

    A = rand(ComplexF64, 10, 10)
    L = LinearMap{ComplexF64}(x -> A*x, 10)
    @test L * v == A * v == mul!(w, L, v)
    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), 10)
    @test L * v == A * v == mul!(w, L, v)
    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), (y, x) -> mul!(y, A', x), 10)
    @test L * v == A * v == mul!(w, L, v)
    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)

    A = Symmetric(rand(ComplexF64, 10, 10))
    L = LinearMap{ComplexF64}(x -> A*x, 10; issymmetric=true)
    @test L * v == A * v == mul!(w, L, v)
    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), 10; issymmetric=true)
    @test L * v == A * v == mul!(w, L, v)
    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)

    A = Hermitian(rand(ComplexF64, 10, 10))
    L = LinearMap{ComplexF64}(x -> A*x, 10; ishermitian=true)
    @test L * v == A * v == mul!(w, L, v)
    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
    L = LinearMap{ComplexF64}((y, x) -> mul!(y, A, x), 10; ishermitian=true)
    @test L * v == A * v == mul!(w, L, v)
    @test adjoint(L) * v ≈ A'v ≈ mul!(w, L', v)
    @test transpose(L) * v ≈ transpose(A)*v ≈ mul!(w, transpose(L), v)
end


text:
"Returns 1xN equally spaced reds going from light to dark for plotting"
function reds(N::Int)
    colors = N==1 ? [RGB(0.8,0.3,0.3)] : range(RGB(1.0,0.4,0.4),stop=RGB(0.6,0.2,0.2),length=N)
    return reshape(colors,1,:)
end

"Returns 1xN equally spaced greens going from light to dark for plotting"
function greens(N::Int)
    colors = N==1 ? [RGB(0.3,0.8,0.3)] : range(RGB(0.4,1.0,0.4),stop=RGB(0.2,0.6,0.2),length=N)
    return reshape(colors,1,:)
end

"Returns 1xN equally spaced blues going from light to dark for plotting"
function blues(N::Int)
    colors = N==1 ? [RGB(0.3,0.3,0.8)] : range(RGB(0.4,0.4,1.0),stop=RGB(0.2,0.2,0.6),length=N)
    return reshape(colors,1,:)
end


text:
// passing integers to printDouble().

int main() {
	printDouble(1);
	return 0 ;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
test_gage_principal_axes(
    Array[[0., 1., 1.],
     [0., 0., 0.]],
    [0.  1.  0.;
     sqrt(2)/2. 0. -sqrt(2)/2.;
     sqrt(2)/2. 0. sqrt(2)/2. ], "test using array of arrays a input")


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
mutable struct TypeA
    n::Int
end

mutable struct TypeB
    x::Float64
end

"""
`fun!(a::TypeA, b::TypeB, c::AbstractVector{Float64})` only mutates `b`.
"""
function fun!(a::TypeA, b::TypeB, c::AbstractVector{Float64})
    b.x = sum(c) / a.n
end

a, b, c = TypeA(10), TypeB(0.0), collect(1.0:10.0)
fun!(a, b, c)
@show a b c;

# %%
@doc fun!

# %%
?fun!

# %%

max_stars_repo_path: 0016/doc string example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/doc string example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/doc string example.jl

max_forks_repo_name: genkuroki/public


text:
cl = ComplexLine(100rand() + im * 100rand())
s = Spherical(cl)
c = Cartesian(cl)
g = Geographic(cl)
r3 = normalize(ℝ³(rand(3)))

@test isapprox(cl, ComplexLine(s))
@test isapprox(cl, ComplexLine(c))
@test isapprox(cl, ComplexLine(g))
@test isapprox(s, Spherical(c))
@test isapprox(s, Spherical(g))
@test isapprox(c, Cartesian(s))
@test isapprox(c, Cartesian(g))
@test isapprox(g, Geographic(s))
@test isapprox(g, Geographic(c))
@test isapprox(g, Geographic(Cartesian(ℝ³(Cartesian(g)))))

@test isapprox(ℝ³(cl), ℝ³(c))
@test isapprox(ComplexLine(r3), ComplexLine(Cartesian(r3)))

number = rand(5:10)
a = [Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2,) for i in 1:number]
b = [Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]

@test isapprox(a, a, atol = TOLERANCE)
@test isapprox(a, b, atol = TOLERANCE) == false

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
#Get mantissa and exponent parts of a number
function splitnum(x::Real)
    parts::Array{String,1} = split(string(float64(x)),"e")
    if length(parts) == 1
        if x != 0
            exponent::Int16 = floor(log10(abs(parsefloat(parts[1]))))
            mantissa::Float64 = parsefloat("$(parts[1])e$(-exponent)")
        else
            mantissa = 0.0
            exponent = 0
        end
    else
        mantissa = parsefloat(parts[1])
        exponent = parseint(parts[2])
    end
    return mantissa, exponent
end

#Format to n significant figures
function to_sf(x::Real,num_sf::Int64;expformat::String="e%",force=false)
    #Usage:
    #   s = to_sf(x,nsf)
    #   s = to_sf(..., force=force)
    #   s = to_sf(..., expformat=expformat)
    #       x           real number to be converted
    #       nsf         number of significant figures to be presented in
    #       force       true if scientific notation is to be enforced for small
    #                   exponents in [-1,3] (default false to allow casual notation
    #                   sometimes)
    #       expformat   string for exponent to be displayed, occurrence of '%' denotes
    #                   position of exponent value
    #       s           output string of formatted number

    #Check input
    if num_sf < 1
        error("Number of significant figures must be an integer greater or equal to 1.")
    end

    #Account for zeros
    if x == 0; return "0"; end

    #Split number into mantissa and exponent parts
    mantissa::Float64, exponent::Int16 = splitnum(x)
    #Round to correct number of significant figures
    digits::String = string(int64(round(parsefloat("$(mantissa)e$(num_sf-1)"))))

    if !force
        #Allow standard notation for...
        if exponent == -1;
            #... numbers in [0.1,1), and...
            return "0."*digits
        elseif num_sf-1 <= exponent <= 3
            #... integers in [1,9999]
            return rpad(digits,exponent+1,"0")
        end
    end

    #Use scientific notation
    if length(digits) == 1
        #no decimal point for single digit
        numstr::String = digits
    else
        #Include decimal point
        numstr = "$(digits[1])."*digits[2:end]
    end
    if exponent != 0
        #append exponent
        expsymb = replace(expformat,"%",string(exponent))
        numstr *= expsymb
    end

    return numstr
end

#Format to n decimal places
function to_dp(x::Real,num_dp::Int64)
    #Usage:
    #   s = to_dp(x,ndp)
    #       x       real number to be converted
    #       ndp     number of decimal places to be presented in. Can be negative
    #               to round to places left of decimal point
    #       s       output string of formatted number

    if num_dp > 0
        return eval(parse("@sprintf(\"%.$(num_dp)f\",$(float64(x)))"))
    else
        #Get significant digits
        digits::String = string(int64(round(parsefloat("$(float64(x))e$(num_dp)"))))
        if digits == "0"
            return digits
        else
            #pad with zeros
            return digits*"0"^abs(num_dp)
        end
    end
end


text:
# Top-level file

"""
Support for sparse tensors (i.e., multi-dimensional arrays).
Provides `AbstractSparseTensor` and subtypes.
"""
module SparseTensors

# AbstractArray interface
import Base: size, getindex, setindex!
import Base: IndexStyle

# Overloads for specializing outputs
import Base: show, summary, similar

# Overloads for improving efficiency
import Base: findall, iterate

# Exported types
export AbstractSparseTensor, SparseTensorCOO, SparseTensorDOK

# AbstractSparseTensor interface
export dropstored!, numstored, storedindices, storedvalues, storedpairs

# Exported functions
export indtype

include("abstract.jl")
include("coo.jl")
include("dok.jl")

end


text:
module BQL

import Base: show, showerror, >>, !, +, -, <, <=, >, >=, ==, !=, &, |, *, /, //

export
    Combinator,
    IllegalInputError,
    IllegalOperandsError,
    IllegalContextError,
    Const,
    This,
    Field,
    Fields,
    Var,
    Vars,
    Select,
    Given,
    UnaryOp,
    BinaryOp,
    Sieve,
    Count,
    Min,
    Max,
    First,
    vars

abstract Combinator

call(C::Combinator, x) = call(C, x, Dict())

call(C::Combinator, x, ctx) = error("$C: not implemented")

vars(C::Combinator) = Set()

immutable IllegalInputError <: Exception
    C::Combinator
    x
end

function showerror(io::IO, err::IllegalInputError)
    print(io, "IllegalInputError: cannot apply ", err.C, " to ")
    dump(io, err.x, 2)
end

immutable IllegalOperandsError <: Exception
    C::Combinator
    y
end

function showerror(io::IO, err::IllegalOperandsError)
    print(io, "IllegalOperandsError: cannot evaluate ", err.C, " on ")
    dump(io, err.y, 2)
end

immutable IllegalContextError <: Exception
    C::Combinator
end

function showerror(io::IO, err::IllegalContextError)
    print(io, "IllegalContextError: cannot evaluate ", err.C, " in the given context")
end

immutable Const <: Combinator
    val
end

call(C::Const, x, ctx) = x == nothing ? nothing : C.val

immutable This <: Combinator
end

call(C::This, x, ctx) = x

immutable Field <: Combinator
    name::AbstractString
end

Field(name::Symbol) = Field(string(name))

call(C::Field, x, ctx) =
    isa(x, Dict) && haskey(x, C.name) ? x[C.name] :
    x == nothing ? nothing :
    throw(IllegalInputError(C, x))

Fields(names...) = map(Field, names)

immutable Var <: Combinator
    name::AbstractString
end

Var(name::Symbol) = Var(string(name))

call(C::Var, x, ctx) =
    haskey(ctx, C.name) ? ctx[C.name] :
    throw(IllegalContextError(C))

vars(C::Var) = Set([C.name])

Vars(names...) = map(Var, names)

immutable Compose <: Combinator
    F::Combinator
    G::Combinator
end

>>(F::Combinator, G::Combinator) = Compose(F, G)

show(io::IO, C::Compose) = print(io, C.F, " >> ", C.G)

function call(C::Compose, x, ctx)
    y = C.F(x, ctx)
    if isa(y, Array)
        z = []
        for yi in y
            zi = C.G(yi, ctx)
            if isa(zi, Array)
                append!(z, zi)
            elseif zi != nothing
                push!(z, zi)
            end
        end
    else
        z = C.G(y)
    end
    z
end

vars(C::Compose) = union(vars(C.F), vars(C.G))

immutable Select <: Combinator
    name_to_F::Dict{AbstractString, Combinator}
end

show(io::IO, C::Select) = print(io, "Select(", join([string(name=>F) for (name, F) in C.name_to_F], ","), ")")

Select(ps::Pair...) = Select(Dict{AbstractString, Combinator}(map(p -> string(p.first) => p.second, ps)))

call(C::Select, x, ctx) = Dict([name => F(x, ctx) for (name, F) in C.name_to_F])

vars(C::Select) = union(map(vars, values(C.name_to_F))...)

immutable Given <: Combinator
    F::Combinator
    name_to_F::Dict{AbstractString, Combinator}
end

show(io::IO, C::Given) =
    print(io, "Given(", C.F, ",", join([string(name=>F) for (name, F) in C.name_to_F], ","), ")")

Given(C::Combinator, ps::Pair...) =
    Given(C, Dict{AbstractString, Combinator}(map(p -> string(p.first) => p.second, ps)))

function call(C::Given, x, ctx)
    ctx = copy(ctx)
    for (name, F) in C.name_to_F
        ctx[name] = F(x, ctx)
    end
    C.F(x, ctx)
end

vars(C::Given) = setdiff(vars(C.F), keys(C.name_to_F))

immutable UnaryOp <: Combinator
    op
    F::Combinator
end

show(io::IO, C::UnaryOp) = print(io, C.op, C.F)

function call(C::UnaryOp, x, ctx)
    y = C.F(x, ctx)
    return (
        isa(y, Array) ? map(C.op, y) :
        y != nothing ? C.op(y) :
        nothing)
end

vars(C::UnaryOp) = vars(C.F)

for op in [:(!), :(+), :(-)]
    @eval $op(F::Combinator) = UnaryOp($op, F)
end

immutable BinaryOp <: Combinator
    op
    F::Combinator
    G::Combinator
end

show(io::IO, C::BinaryOp) = print(io, C.F, C.op, C.G)

function call(C::BinaryOp, x, ctx)
    y1 = C.F(x, ctx)
    y2 = C.G(x, ctx)
    return _binop(C.op, y1, y2)
end

_binop(op, y1, y2) =
    y1 == nothing || y2 == nothing ? nothing :
    isa(y1, Array) && !isa(y2, Array) ? map(y1i -> _binop(op, y1i, y2), y1) :
    !isa(y1, Array) && isa(y2, Array) ? map(y2i -> _binop(op, y1, y2i), y2) :
    op == (+) && isa(y1, Dict) && isa(y2, Dict) ? merge(y1, y2) :
    op(y1, y2)

vars(C::BinaryOp) = union(vars(C.F), vars(C.G))

const Scalar = Union{Void, Bool, Number, AbstractString}
for op in [:(<), :(<=), :(>), :(>=), :(==), :(!=), :(&), :(|), :(+), :(-), :(*), :(/)]
    @eval begin
        $op(F::Combinator, G::Combinator) = BinaryOp($op, F, G)
        $op(F::Combinator, g::Scalar) = BinaryOp($op, F, Const(g))
        $op(f::Scalar, G::Combinator) = BinaryOp($op, Const(f), G)
    end
end

//(F::Combinator, G::Combinator) = BinaryOp(div, F, G)
//(F::Combinator, g::Scalar) = BinaryOp(div, F, Const(g))
//(f::Scalar, G::Combinator) = BinaryOp(div, Const(f), G)

immutable Sieve <: Combinator
    P::Combinator
end

function call(C::Sieve, x, ctx)
    p = C.P(x, ctx)
    return (
        p == true ? x :
        p == nothing || p == false ? nothing :
        throw(IllegalOperandsError(C, p)))
end

vars(C::Sieve) = vars(C.P)

immutable Count <: Combinator
    F::Combinator
end

function call(C::Count, x, ctx)
    y = C.F(x, ctx)
    return isa(y, Array) ? length(y) : throw(IllegalOperandsError(C, y))
end

vars(C::Count) = vars(C.F)

immutable Min <: Combinator
    F::Combinator
end

function call(C::Min, x, ctx)
    y = C.F(x, ctx)
    return (
        isa(y, Array) && length(y) > 0 ? minimum(y) :
        isa(y, Array) ? nothing :
        throw(IllegalOperandsError(C, y)))
end

vars(C::Min) = vars(C.F)

immutable Max <: Combinator
    F::Combinator
end

function call(C::Max, x, ctx)
    y = C.F(x, ctx)
    return (
        isa(y, Array) && length(y) > 0 ? maximum(y) :
        isa(y, Array) ? nothing :
        throw(IllegalOperandsError(C, y)))
end

vars(C::Max) = vars(C.F)

immutable First <: Combinator
    F::Combinator
    N
end

First(F) = First(F, nothing)

function call(C::First, x, ctx)
    y = C.F(x, ctx)
    if isa(C.N, Void)
        return (
            isa(y, Array) && length(y) > 0 ? y[1] :
            isa(y, Array) ? nothing :
            throw(IllegalOperandsError(C, y)))
    else
        N = isa(C.N, Combinator) ? C.N(x, ctx) : C.N
        return isa(y, Array) ? y[1:N] : throw(IllegalOperandsError(C, y))
    end
end

vars(C::First) = vars(C.F)

end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
A package for solving various marriage market matching models.

Two models are supported:

- `StaticMatch`: the static frictionless matching model of Choo & Siow (2006)
- `SearchMatch`: the search and matching model of Shimer & Smith (2000)

`SearchMatch` is generalized in a few ways beyond the model presented in the paper.
This includes a match-specific "love" shock to make the matching probabilistic.
"""
module MarriageMarkets

export StaticMatch, estimate_static_surplus, SearchMatch, SearchClosed, SearchInflow

include("static-match.jl")
include("search-match.jl")

end # module


text:
# compute moving average using p most recent values, including current value
function ma(x, p)
    m = similar(x)
    for i = p:size(x,1)
        m[i] = mean(x[i-p+1:i])
    end
    return m
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
module physics

include("num.jl")
using .num
using Elliptic

export dos3D, dosbethe, baredos, fermi, zeta

function dos3D(t::Float64,ω::Float64)
    z = range(-0.995,length=250,stop=0.995)
    z = convert(Array{Float64,1},z);
    ma = 1. .- ((ω / (2. * t) .+ z) ./ 2.).^2.
    km = zeros(Float64,length(ma))
    for (imax,nma) in enumerate(ma)
        if nma > 0. && nma < 1.
            km[imax] = Elliptic.K.(nma)
        end
    end
    intx = km ./ sqrt.(1. .- z.^2)
    return num.trapz(z,intx) / (π^3. * 2. * t)
end

function dosbethe(t::Float64,ω::Float64)
    if ω^2 < 4*t^2
        return (sqrt(4*t^2 - ω^2) / (2*π*t^2))
    else
        return 0.0
    end
end

function baredos(typedos::String,t::Float64,ω::Float64)
    if typedos == "cubic"
        return dos3D(t,ω)
    elseif typedos == "bethe"
        return dosbethe(t,ω)
    else
        println(error("current typedos is either cubic or bethe"))
    end
end

function fermi(ω::Float64,T::Float64)
    kB = 8.61733034e-5;
    if T==0
        if ω <= 0
            return 1.0
        else
            return 0.0
        end
    else
        β = 1. / kB / T
        return 1. ./ ( exp(ω*β) + 1.)
    end
end

zeta(ω,zeroplus) = ω + im * zeroplus

end

max_stars_repo_name: muhammadgaffar/SKRIPSI

max_issues_repo_name: muhammadgaffar/SKRIPSI

max_forks_repo_name: muhammadgaffar/SKRIPSI


text:
module Constants
using OffsetArrays

export get_offset_constants, get_offset_axial_constants, init!, ConstantContext
mutable struct ConstantContext
    nmax::Int64
    bcof::Array{Float64,2}
    fnr::Array{Float64,1}
    monen::Array{Int64,1}
    vwh_coef::Array{Float64,4}
    vcc::Array{Float64,3}
    fnm1::Array{Float64,2}
    fn::Array{Float64,2}
    fnp1::Array{Float64,2}
end

function get_offset_constants(ctx::ConstantContext)
    notd = ctx.nmax
    nbc = 6notd + 6
    monen = OffsetArray(ctx.monen, (-2notd):(2notd))
    bcof = OffsetArray(ctx.bcof, 0:nbc, 0:nbc)
    fnr = OffsetArray(ctx.fnr, 0:(2nbc))
    vwh_coef = OffsetArray(ctx.vwh_coef, (-notd):notd, 1:notd, -1:1, -1:1)
    return bcof, fnr, monen, vwh_coef
end

function get_offset_axial_constants(ctx::ConstantContext)
    nmax = size(ctx.fn)[2]
    vcc = OffsetArray(ctx.vcc, 1:nmax, 1:(nmax * (nmax + 2)), 0:(2nmax))
    fnm1 = OffsetArray(ctx.fnm1, (-nmax):nmax, 1:nmax)
    fn = OffsetArray(ctx.fn, (-nmax):nmax, 1:nmax)
    fnp1 = OffsetArray(ctx.fnp1, (-nmax):nmax, 1:nmax)
    return vcc, fnm1, fn, fnp1
end

function init()
    ctx = ConstantContext(
        0,
        Array{Float64,2}(undef, 0, 0),
        Array{Float64,1}(undef, 0),
        Array{Int64,1}(undef, 0),
        Array{Float64,4}(undef, 0, 0, 0, 0),
        Array{Float64,3}(undef, 0, 0, 0),
        Array{Float64,2}(undef, 0, 0),
        Array{Float64,2}(undef, 0, 0),
        Array{Float64,2}(undef, 0, 0),
    )
    return ctx
end

function init(notd)
    ctx = init()
    init!(ctx, notd)
    return ctx
end

function init!(ctx::ConstantContext, notd)
    if notd <= ctx.nmax
        return
    end
    ctx.nmax = notd
    nbc = 6notd + 6
    monen = OffsetArray(zeros(Int64, 4notd + 1), (-2notd):(2notd))
    bcof = OffsetArray(zeros(nbc + 1, nbc + 1), 0:nbc, 0:nbc)
    fnr = OffsetArray(zeros(2nbc + 1), 0:(2nbc))
    monen[-2notd] = 1
    for n in (-2notd + 1):(2notd)
        monen[n] = -monen[n - 1]
    end
    for n in 0:(2nbc)
        fnr[n] = √n
    end
    bcof[0, 0] = 1
    for n in 0:(nbc - 1)
        for l in (n + 1):nbc
            bcof[n, l] = fnr[n + l] * bcof[n, l - 1] / fnr[l]
            bcof[l, n] = bcof[n, l]
        end
        bcof[n + 1, n + 1] = fnr[2n + 2] * fnr[2n + 1] * bcof[n, n] / fnr[n + 1] / fnr[n + 1]
    end

    vwh_coef = OffsetArray(zeros(2notd + 1, notd, 3, 3), (-notd):notd, 1:notd, -1:1, -1:1)

    @simd for n in 1:notd
        nn1 = n * (n + 1)
        np1 = n + 1
        nm1 = n - 1
        fnorm1 = -0.5 / fnr[2n + 1] / fnr[n] / fnr[n + 1]
        fnorm2 = -0.5 * fnr[2n + 1] / fnr[n] / fnr[n + 1]
        m = -n
        mp1 = m + 1
        mm1 = m - 1
        vwh_coef[m, n, 1, 1] = -fnorm1 * n * fnr[np1 + m] * fnr[np1 + mp1]
        vwh_coef[m, n, 1, -1] = fnorm1 * np1 * fnr[n - m] * fnr[nm1 - m]
        vwh_coef[m, n, -1, 1] = fnorm1 * n * fnr[np1 - m] * fnr[np1 - mm1]
        vwh_coef[m, n, -1, -1] = 0
        vwh_coef[m, n, 0, 1] = fnorm1 * n * fnr[np1 + m] * fnr[np1 - m]
        vwh_coef[m, n, 0, -1] = 0
        vwh_coef[m, n, 1, 0] = -fnorm2 * fnr[n - m] * fnr[np1 + m]
        vwh_coef[m, n, -1, 0] = -0
        vwh_coef[m, n, 0, 0] = -fnorm2 * m
        @simd for m in (-n + 1):(n - 1)
            mp1 = m + 1
            mm1 = m - 1
            vwh_coef[m, n, 1, 1] = -fnorm1 * n * fnr[np1 + m] * fnr[np1 + mp1]
            vwh_coef[m, n, 1, -1] = fnorm1 * np1 * fnr[n - m] * fnr[nm1 - m]
            vwh_coef[m, n, -1, 1] = fnorm1 * n * fnr[np1 - m] * fnr[np1 - mm1]
            vwh_coef[m, n, -1, -1] = -fnorm1 * np1 * fnr[n + m] * fnr[nm1 + m]
            vwh_coef[m, n, 0, 1] = fnorm1 * n * fnr[np1 + m] * fnr[np1 - m]
            vwh_coef[m, n, 0, -1] = fnorm1 * np1 * fnr[n + m] * fnr[n - m]
            vwh_coef[m, n, 1, 0] = -fnorm2 * fnr[n - m] * fnr[np1 + m]
            vwh_coef[m, n, -1, 0] = -fnorm2 * fnr[n + m] * fnr[np1 - m]
            vwh_coef[m, n, 0, 0] = -fnorm2 * m
        end
        m = n
        mp1 = m + 1
        mm1 = m - 1
        vwh_coef[m, n, 1, 1] = -fnorm1 * n * fnr[np1 + m] * fnr[np1 + mp1]
        vwh_coef[m, n, 1, -1] = 0
        vwh_coef[m, n, -1, 1] = fnorm1 * n * fnr[np1 - m] * fnr[np1 - mm1]
        vwh_coef[m, n, -1, -1] = -fnorm1 * np1 * fnr[n + m] * fnr[nm1 + m]
        vwh_coef[m, n, 0, 1] = fnorm1 * n * fnr[np1 + m] * fnr[np1 - m]
        vwh_coef[m, n, 0, -1] = 0
        vwh_coef[m, n, 1, 0] = -0
        vwh_coef[m, n, -1, 0] = -fnorm2 * fnr[n + m] * fnr[np1 - m]
        vwh_coef[m, n, 0, 0] = -fnorm2 * m
    end

    ctx.monen = parent(monen)
    ctx.bcof = parent(bcof)
    ctx.fnr = parent(fnr)
    return ctx.vwh_coef = parent(vwh_coef)
end

end  # module Constants


text:
 # Full width cells in Jupyter notebook
full_width_notebook(width=100) = display(HTML("<style>.container { width:$width% !important; }</style>"))

# Display of action nodes.
function MCTS.node_tag(s::AST.ASTState)
    state_str::String = "0x"*string(s.hash, base=16)
    if s.terminal
        return "Terminal [$state_str]."
    else
        return state_str
    end
end

# Display of state nodes.
function MCTS.node_tag(a::AST.ASTSeedAction)
    if a == action # selected optimal action
        return "—[$(string(a))]—"
    else
        return "[$(string(a))]"
    end
end

"""
Visualize MCTS tree structure for AST MDP.
"""
function visualize(planner::MCTS.DPWPlanner)
    tree = search!(planner; return_tree=true)
    d3 = visualize(tree)
    return d3::D3Tree
end

function visualize(tree::MCTS.DPWTree)
    d3::D3Tree = D3Tree(tree, init_expand=1)
    return d3::D3Tree
end


text:
####
#### Implementation of the TikTak method.
####

export TikTak

struct TikTak
    quasirandom_N::Int
    initial_N::Int
    θ_min::Float64
    θ_max::Float64
    θ_pow::Float64
end

"""
$(SIGNATURES)

The “TikTak” multistart method, as described in *Arnoud, Guvenen, and Kleineberg (2019)*.

This implements the *multistart* part, can be called with arbitrary local methods, see
[`multistart_minimization`](@ref).

# Arguments

- `quasirandom_N`: the number of quasirandom points for the first pass (using a Sobol
  sequence).

# Keyword arguments

- `keep_ratio`: the fraction of best quasirandom points which are kept

- `θ_min` and `θ_max` clamp the weight parameter, `θ_pow` determines the power it is raised
  to.

The defaults are from the paper cited above.
"""
function TikTak(quasirandom_N; keep_ratio = 0.1, θ_min = 0.1, θ_max = 0.995, θ_pow = 0.5)
    @argcheck 0 < keep_ratio ≤ 1
    TikTak(quasirandom_N, ceil(keep_ratio * quasirandom_N), θ_min, θ_max, θ_pow)
end

function _weight_parameter(t::TikTak, i)
    @unpack initial_N, θ_min, θ_max, θ_pow = t
    clamp((i / initial_N)^θ_pow, θ_min, θ_max)
end

"""
$(SIGNATURES)

Evaluate and return points of an `N`-element Sobol sequence.

When `use_threads`, execution is parallelized using `Threads.@spawn`.
"""
function sobol_starting_points(minimization_problem::MinimizationProblem, N::Integer,
                               use_threads::Bool)
    @unpack objective, lower_bounds, upper_bounds = minimization_problem
    s = SobolSeq(lower_bounds, upper_bounds)
    skip(s, N)                  # better uniformity
    points = Iterators.take(s, N)
    _initial(x) = (location = x, value = objective(x))
    if use_threads
        map(fetch, map(x -> @spawn(_initial(x)), points))
    else
        map(_initial, points)
    end
end

"""
$(SIGNATURES)

Helper function to keep the `N` points with the lowest `value`.
"""
function _keep_lowest(xs, N)
    @argcheck 1 ≤ N ≤ length(xs)
    partialsort(xs, 1:N, by = p -> p.value)
end

"""
$(SIGNATURES)

Solve `minimization_problem` by using `local_method` within `multistart_method`.

When `use_threads`, initial point search is parallelized using `Threads.@spawn`.
"""
function multistart_minimization(multistart_method::TikTak, local_method,
                                 minimization_problem; use_threads = true)
    @unpack quasirandom_N, initial_N, θ_min, θ_max, θ_pow = multistart_method
    quasirandom_points = sobol_starting_points(minimization_problem, quasirandom_N,
                                               use_threads)
    initial_points = _keep_lowest(quasirandom_points, initial_N)
    function _step(visited_minimum, (i, initial_point))
        θ = _weight_parameter(multistart_method, i)
        x = @. (1 - θ) * initial_point.location + θ * visited_minimum.location
        local_minimum = local_minimization(local_method, minimization_problem, x)
        local_minimum ≡ nothing && return visited_minimum
        local_minimum.value < visited_minimum.value ? local_minimum : visited_minimum
    end
    foldl(_step, enumerate(initial_points); init = first(initial_points))
end


text:
module LyonPassiveEar

include("soscascade.jl")
include("sosfilters.jl")
include("agc.jl")
include("utils.jl")
include("lyon_passive_ear.jl")

export lyon_passive_ear

end # module


text:
function PlotSpanningTree(xi,yi,clustersign,Lx,Ly,treepointers,fighandle)

I = treepointers(:,1);
J = treepointers(:,2);

x1 = xi(I);
x2 = xi(J);
y1 = yi(I);
y2 = yi(J);

if clustersign == 1
    c = [1 0 0];
elseif clustersign == -1
    c = [0 0 1];
end

lw = 1.5;
%figure(fighandle)
hold on

for ii = 1:length(I)
    if (x2(ii) - x1(ii) < -Lx/2) && ( abs(y2(ii) - y1(ii)) < Ly/2)
        plot([x1(ii)-Lx,  x2(ii)],[ y1(ii), y2(ii)],'-','Linewidth',lw,'Color',c)
        plot([x1(ii),  x2(ii)+Lx],[ y1(ii) y2(ii)],'-','Linewidth',lw,'Color',c)
    elseif (x2(ii) - x1(ii) > Lx/2) && (abs(y2(ii) - y1(ii)) < Ly/2)
        plot([x2(ii)-Lx,  x1(ii)],[ y2(ii), y1(ii)],'-','Linewidth',lw,'Color',c)
        plot([x2(ii),  x1(ii)+Lx],[ y2(ii) y1(ii)],'-','Linewidth',lw,'Color',c)
    elseif (y2(ii) - y1(ii) < -Ly/2) && (abs(x2(ii) - x1(ii)) < Lx/2)
        plot([x1(ii) x2(ii)],[y1(ii) y2(ii)+Ly],'-','Linewidth',lw,'Color',c)
        plot([x2(ii) x1(ii)],[y2(ii) y1(ii)-Ly],'-','Linewidth',lw,'Color',c)
    elseif (y2(ii) - y1(ii) > Ly/2) && (abs(x2(ii) - x1(ii)) < Lx/2)
        plot([x1(ii) x2(ii)],[y1(ii) y2(ii)-Ly],'-','Linewidth',lw,'Color',c)
        plot([x2(ii) x1(ii)],[y2(ii) y1(ii)+Ly],'-','Linewidth',lw,'Color',c)
    elseif (y2(ii) - y1(ii) < -Ly/2) && (x2(ii) - x1(ii) < -Lx/2)
        plot([x1(ii)-Lx,  x2(ii)],[ y1(ii)-Ly, y2(ii)],'-','Linewidth',lw,'Color',c)
        plot([x1(ii),  x2(ii)+Lx],[ y1(ii) y2(ii)+Ly],'-','Linewidth',lw,'Color',c)
    elseif (x2(ii) - x1(ii) > Lx/2) && (y2(ii) - y1(ii) < -Ly/2)
        plot([x2(ii)-Lx,  x1(ii)],[ y2(ii)+Ly, y1(ii)],'-','Linewidth',lw,'Color',c)
        plot([x2(ii),  x1(ii)+Lx],[ y2(ii) y1(ii)-Ly],'-','Linewidth',lw,'Color',c)
    elseif (x2(ii) - x1(ii) < -Lx/2) && (y2(ii) - y1(ii) > Ly/2)
        plot([x2(ii)+Lx,  x1(ii)],[ y2(ii)-Ly, y1(ii)],'-','Linewidth',lw,'Color',c)
        plot([x2(ii),  x1(ii)-Lx],[ y2(ii) y1(ii)+Ly],'-','Linewidth',lw,'Color',c)
    elseif (x2(ii) - x1(ii) < -Lx/2) && (y2(ii)-y1(ii) <-Ly/2)
        plot([x2(ii)+Lx,  x1(ii)],[ y2(ii)+Ly, y1(ii)],'-','Linewidth',lw,'Color',c)
        plot([x2(ii),  x1(ii)-Lx],[ y2(ii) y1(ii)-Ly],'-','Linewidth',lw,'Color',c)
    elseif (x2(ii) - x1(ii) > Lx/2) && (y2(ii)-y1(ii) > Ly/2)
        plot([x2(ii),  x1(ii)+Lx],[ y2(ii), y1(ii)+Ly],'-','Linewidth',lw,'Color',c)
        plot([x2(ii)-Lx,  x1(ii)],[ y2(ii)-Ly y1(ii)],'-','Linewidth',lw,'Color',c)

    else
        plot([x1(ii) x2(ii)],[y1(ii) y2(ii)],'-','Linewidth',lw,'Color',c)
    end
end

return


text:
using Plots
using BenchmarkTools

#Problem 1(a)
function vfsolve(vnew, kgrid, tolerance, imax)
    α = 0.37
    δ = 0.1
    β = 0.95
    σ = 1.5

    v = vnew .+ 2*tolerance
    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1

    c = kgrid.^α + (1-δ)*kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid.^α + (1-δ)*kgrid .- kgrid' .< 0] .= -Inf

    while maximum(abs.(v - vnew)) > tolerance && i <= imax
        v = vnew;

        (vnew, cartesianindex) = findmax(u .+ β*v', dims = 2);
        i += 1;
    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex, M=i)
end

kupper = 5
klower = 0.05
n = 500
kgrid = collect(range(klower, stop = kupper, length = n))
(v, kprime, kprimeindex, M) = vfsolve(zeros(n), kgrid, 0.001, 1000);
M
#M=95 the output value function is close enough to the true value function
#1(b)
#define a new function that take imax as the only control for the number of iterations
function vfsolve1(vnew, kgrid, imax)
    α = 0.37
    δ = 0.1
    β = 0.95
    σ = 1.5

    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1
    c = kgrid.^α + (1-δ)*kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid.^α + (1-δ)*kgrid .- kgrid' .< 0] .= -Inf

    while i <= imax

        (vnew, cartesianindex) = findmax(u .+ β*vnew', dims = 2);
        i += 1;

    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex)
end

(v1, kprime1, kprimeindex1) = vfsolve1(zeros(n), kgrid, 1);
(v5, kprime5, kprimeindex5) = vfsolve1(zeros(n), kgrid, 5);
(v10, kprime10, kprimeindex10) = vfsolve1(zeros(n), kgrid, 10);
plot(kgrid, v1, label = "v1")
plot!(kgrid, v5, label = "v5")
plot!(kgrid, v10, label = "v10")
plot!(kgrid, v, label = "v")

#1(c)

#define a function to find the value in krgrid which is closest to the input x and return the index
function findindex(x)
    a=findmin(abs.(x .- kgrid), dims=1)
    return a[2]
end

#policy(k) returns the optimal k'
function policy(x)
    e=getindex(findall(kgrid .== x),1)
    return kprime[e]
end

# the function find and plot the time path of k
function findkpath(T, k_0)
    w = findindex(k_0)
    kpath = zeros(T)
    k = kgrid[w]
        for i in 1:T
            k = policy(k)
            kpath[i] = k
        end
    return kpath
end

# # the function find and plot the time path of k

function plotkpath(T, k_0)
w = findindex(k_0)
kpath = zeros(T)
k = kgrid[w]
    for i in 1:T
        k = policy(k)
        kpath[i] = k
    end
time = collect(range(1, T, length = T))
plot(time, kpath, label = "kpath")
end

findkpath(100, 0.5)
plotkpath(100, 0.5)

#the steady state level of capital is around 4
#1(d)
findkpath(100, 4.5)
plotkpath(100, 4.5)
#the steady state level of capital is also around 4
#(e)

#(e)(a)
function vfsolve_E(vnew, kgrid, tolerance, imax)
    α = 0.37
    δ = 0.1
    β = 0.95
    σ = 0.5

    v = vnew .+ 2*tolerance
    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1

    c = kgrid.^α + (1-δ)*kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid.^α + (1-δ)*kgrid .- kgrid' .< 0] .= -Inf

    while maximum(abs.(v - vnew)) > tolerance && i <= imax
        v = vnew;

        (vnew, cartesianindex) = findmax(u .+ β*v', dims = 2);
        i += 1;
    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex, n=i)
end
kupper = 5
klower = 0.05
n = 500
kgrid = collect(range(klower, stop = kupper, length = n))
(v_E, kprime_E, kprimeindex_E, n_E) = vfsolve_E(zeros(n), kgrid, 0.001, 1000)

#(e)(b)
function vfsolve1_E(vnew, kgrid, imax)
    α = 0.37
    δ = 0.1
    β = 0.95
    σ = 0.5

    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1
    c = kgrid.^α + (1-δ)*kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid.^α + (1-δ)*kgrid .- kgrid' .< 0] .= -Inf

    while i <= imax

        (vnew, cartesianindex) = findmax(u .+ β*vnew', dims = 2);
        i += 1;

    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex)
end

(v1_E, kprime1_E, kprimeindex1_E) = vfsolve1_E(zeros(n), kgrid, 1);
(v5_E, kprime5_E, kprimeindex5_E) = vfsolve1_E(zeros(n), kgrid, 5);
(v10_E, kprime10_E, kprimeindex10_E) = vfsolve1_E(zeros(n), kgrid, 10);
plot(kgrid, v1_E, label = "v1_E")
plot!(kgrid, v5_E, label = "v5_E")
plot!(kgrid, v10_E, label = "v10_E")
plot!(kgrid, v_E, label = "v_E")

#(e)(c)
function policy_E(x)
    e=getindex(findall(kgrid .== x),1)
    return kprime_E[e]
end

function findkpath_E(T, k_0)
    w = findindex(k_0)
    kpath = zeros(T)
    k = kgrid[w]
        for i in 1:T
            k = policy_E(k)
            kpath[i] = k
        end
    return kpath
end

function plotkpath_E(T, k_0)
    w = findindex(k_0)
    kpath = zeros(T)
    k = kgrid[w]
        for i in 1:T
            k = policy_E(k)
            kpath[i] = k
        end
    time = collect(range(1, T, length = T))
    plot(time, kpath, label = "kpath")
end

findkpath_E(100, 0.5)
plotkpath_E(100, 0.5)

#the steady state level of capital is around 4

#1(d)
findkpath_E(100, 4.5)
plotkpath_E(100, 4.5)

#the steady state level of capital is aslo around 4

#the steady state k is same as before but the capital converges faster
#1/σ is the eslasticity of intertemperal substitution, the bigger eslasticity means it is harder
# to subsitute consumption intertemperally, as a result, the consumer saves more resulting more faster capital accumulation

max_stars_repo_name: paradoxuu/Dynamic-Programming-Problems

max_issues_repo_name: paradoxuu/Dynamic-Programming-Problems

max_forks_repo_name: paradoxuu/Dynamic-Programming-Problems


text:
module DummyModule end

dummymime = MIME("dummy/dummy")
DummyMime = typeof(dummymime)
Genser.genser_type_for_mime(::Type{DummyMime}) = DummyModule

dummyencoding = Encoding("dummy")
DummyEncoding = typeof(dummyencoding)
Genser.genser_converter_for_encoding(::Type{DummyEncoding}) = DummyEncoding

@testset "Registry" begin
    @testset "types" begin
        @testset "With mime type" begin
            @test genser_type_for_mime(DummyMime) == DummyModule
        end

        @testset "With mime value" begin
            @test genser_type_for_mime(MIME("dummy/dummy")) == DummyModule
        end

        @testset "With mime string" begin
            @test genser_type_for_mime("dummy/dummy") == DummyModule
        end

        @testset "Unregistered" begin
            @test_throws ArgumentError genser_type_for_mime("dummy/unregistered")
        end
    end

    @testset "encoding converters" begin
        @testset "With encoding type" begin
            genser_converter_for_encoding(DummyEncoding) == DummyModule
        end

        @testset "With encoding value" begin
            genser_converter_for_encoding(Encoding("dummy")) == DummyModule
        end

        @testset "With encoding string" begin
            genser_converter_for_encoding("dummy") == DummyModule
        end

        @testset "Unregistered" begin
            @test_throws ArgumentError genser_converter_for_encoding("unregistered")
        end
    end
end


text:
# Recombinations
# ==============
function average(population::Vector{T}) where {T <: Vector}
    obj = zeros(eltype(T), length(population[1]))
    l = length(population)
    for i in 1:l
        obj += population[i]
    end
    return obj./l
end

function marriage(population::Vector{T}) where {T <: Vector}
    s = length(population)
    l = length(population[1])
    obj = zeros(eltype(T), l)
    for i in 1:l
        obj[i] = population[rand(1:s)][i]
    end
    return obj
end

# Strategy recombinations
# =======================
function averageSigma1(ss::Vector{S}) where {S <: Strategy}
    s = copy(ss[1])
    σ = 0.0
    l = length(ss)
    for i in 1:l
        σ += ss[i][:σ]
    end
    s[:σ] = σ/l
    return s
end

function averageSigmaN(ss::Vector{S}) where {S <: Strategy}
    s = copy(ss[1])
    σ = zeros(length(ss[1][:σ]))
    l = length(ss)
    for i in 1:l
        σ += ss[i][:σ]
    end
    s[:σ] = σ./l
    return s
end

# ==================
# Genetic algorithms
# ==================

# Binary crossovers
# -----------------

# Single point crossover
function singlepoint(v1::T, v2::T) where {T <: Vector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    pos = rand(1:l)
    for i in pos:l
        vswap!(c1, c2, i)
    end
    return c1, c2
end

# Two point crossover
function twopoint(v1::T, v2::T) where {T <: Vector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    for i in from:to
        vswap!(c1, c2, i)
    end
    return c1, c2
end

# Uniform crossover
function uniform(v1::T, v2::T) where {T <: Vector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    xch = rand(Bool, l)
    for i in 1:l
        if xch[i]
            vswap!(c1, c2, i)
        end
    end
    return c1, c2
end

# Real valued crossovers
# ----------------------

function discrete(v1::T, v2::T) where {T <: Vector}
    l = length(v1)
    c1 = similar(v1)
    c2 = similar(v2)
    sltc = rand(Bool, 2, l)
    for i in 1:l
        c1[i] = sltc[1,i] ? v1[i] : v2[i]
        c2[i] = sltc[2,i] ? v2[i] : v1[i]
    end
    return c1, c2
end

# Weighted arithmetic mean
function waverage(w::Vector{Float64})
    function wavexvr(v1::T, v2::T) where {T <: Vector}
        c1 = (v1+v2)./w
        return c1, copy(c1)
    end
    return wavexvr
end

# Intermediate recombination
function intermediate(d::Float64 = 0.0)
    function intermxvr(v1::T, v2::T) where {T <: Vector}
        l = length(v1)
        α = (1.0+2d) * rand(l) .- d
        c1 = v2 .+ α .* (v1 - v2)
        α = (1.0+2d) * rand(l) .- d
        c2 = v1 .+ α .* (v2 - v1)
        return c1, c2
    end
    return intermxvr
end

# Line recombination
function line(d::Float64 = 0.0)
    function linexvr(v1::T, v2::T) where {T <: Vector}
        α1, α2 = (1.0+2d) * rand(2) .- d
        c1 = v2 .+ α2 * (v1 - v2)
        c2 = v1 .+ α1 * (v2 - v1)
        return c1, c2
    end
    return linexvr
end

# Permutation crossovers
# ----------------------

# Partially mapped crossover
function pmx(v1::T, v2::T) where {T <: Vector}
    s = length(v1)
    from, to = rand(1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = similar(v1)
    c2 = similar(v2)

    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]

    # Fill in from parents
    for i in vcat(1:from-1, to+1:s)
        # Check conflicting offspring
        in1 = inmap(v1[i], c1, from, to)
        if in1 == 0
            c1[i] = v1[i]
        else
            tmpin = in1
            while tmpin > 0
                tmpin = inmap(c2[in1], c1, from, to)
                in1 = tmpin > 0 ? tmpin : in1
            end
            c1[i] = v1[in1]
        end

        in2 = inmap(v2[i], c2, from, to)
        if in2 == 0
            c2[i] = v2[i]
        else
            tmpin = in2
            while tmpin > 0
                tmpin = inmap(c1[in2], c2, from, to)
                in2 = tmpin > 0 ? tmpin : in2
            end
            c2[i] = v2[in2]
        end
    end
    return c1, c2
end

# Order crossover
function ox1(v1::T, v2::T) where {T <: Vector}
    s = length(v1)
    from, to = rand(1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = zeros(v1)
    c2 = zeros(v2)
    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]
    # Fill in from parents
    k = to+1 > s ? 1 : to+1 #child1 index
    j = to+1 > s ? 1 : to+1 #child2 index
    for i in vcat(to+1:s,1:from-1)
        while in(v1[k],c1)
            k = k+1 > s ? 1 : k+1
        end
        c1[i] = v1[k]
        while in(v2[j],c2)
            j = j+1 > s ? 1 : j+1
        end
        c2[i] = v2[j]
    end
    return c1, c2
end

# Cycle crossover
function cx(v1::T, v2::T) where {T <: Vector}
    s = length(v1)
    c1 = zeros(v1)
    c2 = zeros(v2)

    f1 = true #switch
    k = 1
    while k > 0
        idx = k
        if f1
            #cycle from v1
            while c1[idx] == zero(T)
            c1[idx] = v1[idx]
            c2[idx] = v2[idx]
            idx = inmap(v2[idx],v1,1,s)
            end
        else
            #cycle from v2
            while c2[idx] == zero(T)
            c1[idx] = v2[idx]
            c2[idx] = v1[idx]
            idx = inmap(v1[idx],v2,1,s)
            end
        end
        f1 $= true
        k = inmap(zero(T),c2,1,s)
    end
    return c1,c2
end

# Order-based crossover
function ox2(v1::T, v2::T) where {T <: Vector}
    s = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)

    for i in 1:s
        if rand(Bool)
            idx1 = inmap(v2[i],v1,1,s)
            idx2 = inmap(v1[i],v2,1,s)
            c1[idx1] = zero(T)
            c2[idx2] = zero(T)
        end
    end

    for i in 1:s
        if !in(v2[i],c1)
            tmpin = inmap(zero(T),c1,1,s)
            c1[tmpin] = v2[i]
        end
        if !in(v1[i],c2)
            tmpin = inmap(zero(T),c2,1,s)
            c2[tmpin] = v1[i]
        end
    end
    return c1,c2
end

# Position-based crossover
function pos(v1::T, v2::T) where {T <: Vector}
    s = length(v1)
    c1 = zeros(v1)
    c2 = zeros(v2)

    for i in 1:s
        if rand(Bool)
            c1[i] = v2[i]
            c2[i] = v1[i]
        end
    end

    for i in 1:s
        if !in(v1[i],c1)
            tmpin = inmap(zero(T),c1,1,s)
            c1[tmpin] = v1[i]
        end
        if !in(v2[i],c2)
            tmpin = inmap(zero(T),c2,1,s)
            c2[tmpin] = v2[i]
        end
    end
    return c1,c2
end

# Utils
# =====
function vswap!(v1::T, v2::T, idx::Int) where {T <: Vector}
    val = v1[idx]
    v1[idx] = v2[idx]
    v2[idx] = val
end

function inmap(v::T, c::Vector{T}, from::Int, to::Int) where{T}
    exists = 0
    for j in from:to
        if exists == 0 && v == c[j]
            exists = j
        end
    end
    return exists
end


text:
function _JRect(
    cornerpoint::Point,
    w::Real,
    h::Real,
    color,
    linewidth::Real,
    action::Symbol,
)
    sethue(color)
    setline(linewidth)
    rect(cornerpoint, w, h, action)
    return cornerpoint
end

"""
    1. JRect(cornerpoint::Point, w::Real, h::Real; kwargs...)
    2. JRect(xmin::Int64, ymin::Int64, w::Real, h::Real; kwargs...)
        - same as 1. with `cornerpoint = Point(xmin, ymin)`

Create a rectangle with one corner at cornerpoint with width w and height h and do an action.
You can specify the `linewidth` and the `color` of the rectangle.

# Keywords for all
- `color` = "black"
- `linewidth` = 1
- `action` Defines whether the rectangle should be outlined (`:stroke`) or filled (`:fill`)
"""
JRect(
    cornerpoint::Point,
    w::Real,
    h::Real;
    color = "black",
    linewidth = 1,
    action = :stroke,
) =
    (
        args...;
        cornerpoint = cornerpoint,
        w = w,
        h = h,
        color = color,
        linewidth = linewidth,
        action = action,
    ) -> _JRect(cornerpoint, w, h, color, linewidth, action)

JRect(xmin::Int64, ymin::Int64, w::Real, h::Real; kwargs...) =
    JRect(Point(xmin, ymin), w, h; kwargs...)


text:
# DeclContext
function getParentASTContext(x::DeclContext)
    @check_ptrs x
    return ASTContext(clang_DeclContext_getParentASTContext(x))
end

function addDecl(x::DeclContext, decl::AbstractDecl)
    @check_ptrs x decl
    return clang_DeclContext_addDecl(x, decl)
end

function addDeclInternal(x::DeclContext, decl::AbstractDecl)
    @check_ptrs x decl
    return clang_DeclContext_addDeclInternal(x, decl)
end

function addHiddenDecl(x::DeclContext, decl::AbstractDecl)
    @check_ptrs x decl
    return clang_DeclContext_addHiddenDecl(x, decl)
end

function removeDecl(x::DeclContext, decl::AbstractDecl)
    @check_ptrs x decl
    return clang_DeclContext_removeDecl(x, decl)
end

function containsDecl(x::DeclContext, decl::AbstractDecl)
    @check_ptrs x decl
    return clang_DeclContext_containsDecl(x, decl)
end

function getDeclKindName(x::DeclContext)
    @check_ptrs x
    return unsafe_string(clang_DeclContext_getDeclKindName(x))
end

function getParent(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_getParent(x)
end

function getLexicalParent(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_getLexicalParent(x)
end

function getLookupParent(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_getLookupParent(x)
end

function isClosure(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isClosure(x)
end

function isFunctionOrMethod(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isFunctionOrMethod(x)
end

function isLookupContext(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isLookupContext(x)
end

function isFileContext(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isFileContext(x)
end

function isTranslationUnit(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isTranslationUnit(x)
end

function isRecord(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isRecord(x)
end

function isNamespace(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isNamespace(x)
end

function isStdNamespace(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isStdNamespace(x)
end

function isInlineNamespace(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isInlineNamespace(x)
end

function is_dependent_context(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isDependentContext(x)
end

function isTransparentContext(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isTransparentContext(x)
end

function isExternCContext(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isExternCContext(x)
end

function isExternCXXContext(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_isExternCXXContext(x)
end

function Equals(x::DeclContext, y::DeclContext)
    @check_ptrs x y
    return clang_DeclContext_Equals(x, y)
end

function getPrimaryContext(x::DeclContext)
    @check_ptrs x
    return DeclContext(clang_DeclContext_getPrimaryContext(x))
end

function dumpDeclContext(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_dumpDeclContext(x)
end

function dumpLookups(x::DeclContext)
    @check_ptrs x
    return clang_DeclContext_dumpLookups(x)
end

# cast
function TagDecl(x::DeclContext)
    @check_ptrs x
    return TagDecl(clang_DeclContext_castToTagDecl(x))
end

function RecordDecl(x::DeclContext)
    @check_ptrs x
    return RecordDecl(clang_DeclContext_castToRecordDecl(x))
end

function CXXRecordDecl(x::DeclContext)
    @check_ptrs x
    return CXXRecordDecl(clang_DeclContext_castToCXXRecordDecl(x))
end


text:
#=
Vector
=#

using SimpleDataStructures

sv = SimpleVector{Int}()
dump(sv)
for i = 1:4
    push!(sv, i)
    @show sv.size, sv.capacity
end
sv
sv[1] = 10
resize!(sv, 3)
dump(sv)

SimpleVector([1, 2, 3])

sv = SimpleVector{Int}()
sizehint!(sv, 100)
for i = 1:100
    push!(sv, i)
    @show sv.size, sv.capacity
end

sv = SimpleVector{Int}()
push!(sv, 2.)
push!(sv, "hello")

sv = SimpleVector{Complex{Float64}}()
push!(sv, complex(7, 3))
push!(sv, 1//2)
push!(sv, π)

SimpleVector(1, 2, 3, 4)

x = SimpleVector(1:10...)
insert!(x, 3, 100)
deleteat!(x, 3)


text:
module CbcMathProgSolverInterface

using Cbc.CoinMPInterface

importall MathProgBase.SolverInterface

export CbcMathProgModel,
    CbcSolver,
    model,
    loadproblem!,
    writeproblem,
    updatemodel!,
    setsense!,
    getsense,
    numvar,
    numconstr,
    setvartype!,
    optimize!,
    status,
    getobjval,
    getobjbound,
    getsolution,
    getrawsolver

type CbcMathProgModel <: AbstractMathProgModel
    inner::CoinProblem
    sense
end

immutable CbcSolver <: AbstractMathProgSolver
    options
end
CbcSolver(;kwargs...) = CbcSolver(kwargs)

function CbcMathProgModel(;options...)
    c = CoinProblem()
    setOption(c, "LogLevel", 0)
    for (optname, optval) in options
        setOption(c, string(optname), optval)
    end
    return CbcMathProgModel(c,:Min)
end

model(s::CbcSolver) = CbcMathProgModel(;s.options...)

function loadproblem!(m::CbcMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
    @assert sense == :Min || sense == :Max
    dir = 1
    if sense == :Max
        dir = -1
    end
    m.sense = sense
    LoadMatrix(m.inner, dir, 0.0, obj, collb, colub, rowlb, rowub, A)
end

function writeproblem(m::CbcMathProgModel, filename::String)
    if endswith(filename,".mps")
        WriteFile(m.inner, 3, filename)
    else
        error("Only MPS output supported")
    end
end

updatemodel(m::CbcMathProgModel) = nothing

function setsense(m::CbcMathProgModel,sense)
    if sense != m.sense
        error("CoinMP interface does not permit modifying the problem sense")
    end
end

getsense(m::CbcMathProgModel) = m.sense

numvar(m::CbcMathProgModel) = GetColCount(m.inner)
numconstr(m::CbcMathProgModel) = GetRowCount(m.inner)

function setvartype!(m::CbcMathProgModel,vartype::Vector{Symbol})
    ncol = numvar(m)
    @assert length(vartype) == ncol
    coltype = Array(Uint8,ncol)
    for i in 1:ncol
        if vartype[i] == :Int
            coltype[i] = 'I'
        elseif vartype[i] == :Cont
            coltype[i] = 'C'
        elseif vartype[i] == :Bin
            coltype[i] = 'B'
        else
            error("Unsupported variable type $(vartype[i]) present")
        end
    end
    LoadInteger(m.inner,coltype)
end

optimize!(m::CbcMathProgModel) = OptimizeProblem(m.inner)

function status(m::CbcMathProgModel)
    stat = GetSolutionText(m.inner)
    # CoinMP status reporting is faulty,
    # add logic from CbcModel.cpp.
    objval = GetObjectValue(m.inner)
    if stat == "Optimal solution found"
        if objval < 1e30
            return :Optimal
        else
            return :Infeasible
        end
    elseif stat == "Problem primal infeasible"
        return :Infeasible
    elseif stat == "Problem dual infeasible" # what does this mean for MIP??
        return :Unbounded
    elseif stat == "Stopped on iterations" || stat == "Stopped by user"
        return :UserLimit
    elseif stat == "Stopped due to errors"
        return :Error
    else
        error("Internal library error")
    end
end

getobjval(m::CbcMathProgModel) = GetObjectValue(m.inner)

getobjbound(m::CbcMathProgModel) = GetMipBestBound(m.inner)

getsolution(m::CbcMathProgModel) = GetSolutionValues(m.inner)

getrawsolver(m::CbcMathProgModel) = m.inner

end