text:
abstract type AbstractTransition end

"""
    DiscreteTransition{RT, GT, IT⁻, IT⁺, WT} <: AbstractTransition

Type that encodes a discrete transition with an affine assignment of the form:

```math
    post_d(X) = (R(X ∩ G ∩ I⁻) ⊕ W) ∩ I⁺
```
where ``I⁻`` and ``I⁺``  are invariants at the source and the target locations
respectively, ``G ⊆ \\mathbb{R}^n`` is the guard set of the transition, the assignment
is of the form ``x^+ := Rx + w``, ``w ∈ W``, ``x^+ ∈ \\mathbb{R}^m`` are the values after
the transition, ``R ∈ \\mathbb{R}^{m\\times n}`` is the assignment map (or reset map)
and ``W ⊆ \\mathbb{R}^m`` is a closed and bounded convex set of non-deterministic inputs.

### Fields

- `R`  -- assignment map
- `W`  -- non-deterministic inputs
- `G`  -- guard set of the transition from the source location to the target location
- `I⁻` -- invariant at the source location
- `I⁺` -- invariant at the target location
"""
struct DiscreteTransition{RT, WT, GT, IT⁻, IT⁺} <: AbstractTransition
    R::RT
    W::WT
    G::GT
    I⁻::IT⁻
    I⁺::IT⁺
end

reset_map(tr::DiscreteTransition{<:AbstractMatrix})  = ConstrainedAffineMap(tr.R, tr.W)
reset_map(tr::DiscreteTransition{<:IdentityMap}) = ConstrainedIdentityMap(dim(tr.W), tr.W)
HybridSystems.guard(tr::DiscreteTransition) = tr.G
source_invariant(tr::DiscreteTransition) = tr.I⁻
target_invariant(tr::DiscreteTransition) = tr.I⁺

# constructor when the reset map is the identity
function DiscreteTransition(; guard::GT, source_invariant::IT⁻, target_invariant::IT⁺) where {GT, IT⁻, IT⁺}
    n = dim(source_invariant)
    return DiscreteTransition(IdentityMap(n), ZeroSet(n), guard, source_invariant, target_invariant)
end

# constructor from hybrid system given source and target modes id's and the transition t
function DiscreteTransition(H::HybridSystem, transition)
    assignment = resetmap(H, transition)
    R = _state_matrix(assignment)
    W = _affine_term(assignment)
    G = stateset(assignment) # HybridSystems.guard(H, transition) doesn't work
    I⁻ = stateset(H, source(H, transition))
    I⁺ = stateset(H, target(H, transition))
    return DiscreteTransition(R, W, G, I⁻, I⁺)
end

# -----------------------------------------------
# Getter functions for MathematicalSystems maps
# (see also MathematicalSystems#209)
# -----------------------------------------------

_state_matrix(m::IdentityMap) = m       # no-op, to dispatch on identity assignments
_state_matrix(m::ConstrainedIdentityMap) = IdentityMap(statedim(m))
_state_matrix(m::MathematicalSystems.LinearMap) = m.A
_state_matrix(m::ConstrainedLinearMap) = m.A
_state_matrix(m::MathematicalSystems.AffineMap) = m.A
_state_matrix(m::ConstrainedAffineMap) = m.A

_affine_term(m::IdentityMap) = ZeroSet(statedim(m)) # there is no numeric type in these maps; use ZeroSet default
_affine_term(m::ConstrainedIdentityMap) = ZeroSet(statedim(m))
_affine_term(m::MathematicalSystems.LinearMap{N}) where {N} = ZeroSet{N}(statedim(m))
_affine_term(m::ConstrainedLinearMap{N}) where {N} = ZeroSet{N}(statedim(m))
_affine_term(m::MathematicalSystems.AffineMap{N, <:AbstractVector}) where {N} = Singleton(m.c)
_affine_term(m::MathematicalSystems.AffineMap{N, <:LazySet}) where {N} = m.c
_affine_term(m::ConstrainedAffineMap{N, MT, VT, <:AbstractVector}) where {N, MT, VT} = Singleton(m.c)
_affine_term(m::ConstrainedAffineMap{N, MT, VT, <:LazySet}) where {N, MT, VT} = Singleton(m.c)

LazySets.dim(c::ConstrainedResetMap) = statedim(c)

function _state_matrix(m::Union{<:MathematicalSystems.ResetMap{N},
                                <:MathematicalSystems.ConstrainedResetMap{N}}) where {N}
    n = dim(m)
    v = ones(N, n)
    for i in keys(m.dict)
        v[i] = zero(N)
    end
    return Diagonal(v)
end

function _affine_term(m::Union{<:MathematicalSystems.ResetMap{N},
                               <:MathematicalSystems.ConstrainedResetMap{N}}) where {N}
    n = dim(m)
    b = sparsevec(Int[], N[], n)
    for (i, val) in m.dict
        b[i] = val
    end
    return Singleton(b)
end

# -------------------------------------------
# Applying reset maps
# -------------------------------------------

@inline function _apply_reset(R::AbstractMatrix, P::HPolytope, W::ZeroSet, ::HRepIntersection)
    return linear_map(R, P) # TODO pass options here
end

@inline function _apply_reset(R::AbstractMatrix, P::HPolytope, W::AbstractVector, ::HRepIntersection)
    return affine_map(R, P, W) # TODO pass options here
end

@inline function _apply_reset(R::AbstractMatrix, P::HPolytope, W::LazySet, ::HRepIntersection)
    Ylm = linear_map(R, P) # TODO pass options here
    return minkowski_sum(Ylm, W)
end

# we assume that the reach-set is bounded: it's either a polytopic set
# or a union of polytopic sets
# => the return type is either the empty set or an HPolytope unless otherwise stated
# for examples in cases involving unions this method may return a union of polytopes
apply(tr::DiscreteTransition,
      R::AbstractLazyReachSet,
      method::AbstractIntersectionMethod) = apply(tr, set(R), method)

# ==========================================
# Intersection method: HRepIntersection
# ==========================================

# -------------------------------------------
# Cases where the reset map is the identity
# -------------------------------------------

# all sets are polyhedral
function  apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
                X::AbstractPolytope{N},
                method::HRepIntersection) where {N,
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    success, out = _intersection(X, tr.G, tr.I⁻, tr.I⁺, method)
    return success ? HPolytope(out) : EmptySet(dim(X))
end

# the source invariant is the set union of half-spaces, the other sets are polyhedral
function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
               X::AbstractPolytope{N},
               method::HRepIntersection) where {N,
                                GT<:AbstractPolyhedron{N},
                                VN,
                                IT⁻<:UnionSetArray{N, LinearConstraint{N, VN}},
                                IT⁺<:AbstractPolyhedron{N}}

    success, aux = _intersection(X, tr.G, tr.I⁺, method)
    !success && return EmptySet(dim(X))

    clist_I⁻ = tr.I⁻.array
    out = Vector{HPolytope{N, Vector{N}}}()

    # intersect Y := X ∩ G ∩ I⁺ with each half-space ci in I⁻
    # only store those polytpes for which Y ∩ ci is non-empty
    for ci in clist_I⁻
        Y_ci = vcat(aux, ci)
        success = remove_redundant_constraints!(Y_ci)
        !success && continue
        push!(out, HPolytope(copy(Y_ci)))
    end
    return isempty(out) ? EmptySet(dim(X)) : UnionSetArray(out)
end

# the source invariant is the set union of half-spaces, the other sets are polyhedral
# and the set X is the convex hull of polytopes
function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
               X::ConvexHullArray{N, PT},
               method::HRepIntersection) where {N,
                                PT<:AbstractPolytope{N},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    vlist = vertices_list(X) # TODO pass the backend as an option
    Xhrep = tohrep(vlist)
    return apply(tr, Xhrep, method)
end

# the guard is the set union of half-spaces, the source invariant is polyhedral,
# and the target invariant is universal
function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
                  X::AbstractPolytope{N},
                  method::HRepIntersection) where {N, VN,
                                                   GT<:UnionSetArray{N, HalfSpace{N, VN}},
                                                   IT⁻<:AbstractPolyhedron{N},
                                                   IT⁺<:Universe{N}}

    success, aux = _intersection(X, tr.I⁻, method)
    !success && return EmptySet(dim(X))

    clist_G = tr.G.array
    out = Vector{HPolytope{N, Vector{N}}}()

    # intersect Y := X ∩ I⁻ with each half-space ci in G
    # only store those polytpes for which Y ∩ ci is non-empty
    for ci in clist_G
        Y_ci = vcat(aux, ci)
        success = remove_redundant_constraints!(Y_ci)
        !success && continue
        push!(out, HPolytope(copy(Y_ci)))
    end

    return !isempty(out) ? UnionSetArray(out) : EmptySet(dim(X))
end

# -------------------------------------------
# Cases where the reset map is RX ⊕ W
# -------------------------------------------

function _apply_hrep(R, W, G, I⁻, I⁺, X, method::HRepIntersection)
    success, aux = _intersection(X, G, I⁻, method)
    !success && return EmptySet(dim(X))

    P = HPolytope(aux)
    Y = _apply_reset(R, P, W, method)

    success, out = _intersection(Y, I⁺, method)
    return success ? HPolytope(out) : EmptySet(dim(X))
end

# cases with W = 0
function apply(tr::DiscreteTransition{<:AbstractMatrix, <:ZeroSet, GT, IT⁻, IT⁺},
               X::AbstractPolytope{N},
               method::HRepIntersection) where {N, GT<:AbstractPolyhedron{N},
                                                   IT⁻<:AbstractPolyhedron{N},
                                                   IT⁺<:AbstractPolyhedron{N}}
    return _apply_hrep(tr.R, tr.W, tr.G, tr.I⁻, tr.I⁺, X, method)
end

# cases where the inputs are deterministic (W is a vector)
function apply(tr::DiscreteTransition{<:AbstractMatrix, <:AbstractVector, GT, IT⁻, IT⁺},
               X::AbstractPolytope{N},
               method::HRepIntersection) where {N, GT<:AbstractPolyhedron{N},
                                                   IT⁻<:AbstractPolyhedron{N},
                                                   IT⁺<:AbstractPolyhedron{N}}
    return _apply_hrep(tr.R, tr.W, tr.G, tr.I⁻, tr.I⁺, X, method)
end

# cases where the inputs are non-deterministic (W is a set)
function apply(tr::DiscreteTransition{<:AbstractMatrix, <:LazySet, GT, IT⁻, IT⁺},
               X::AbstractPolytope{N},
               method::HRepIntersection) where {N, GT<:AbstractPolyhedron{N},
                                                   IT⁻<:AbstractPolyhedron{N},
                                                   IT⁺<:AbstractPolyhedron{N}}
    return _apply_hrep(tr.R, tr.W, tr.G, tr.I⁻, tr.I⁺, X, method)
end

# ==========================================
# Intersection method: BoxIntersection
# ==========================================

function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
                X::AbstractPolytope{N},
                ::BoxIntersection) where {N,
                                       GT<:AbstractPolyhedron{N},
                                       VN,
                                       IT⁻<:UnionSetArray{N, LinearConstraint{N, VN}},
                                       IT⁺<:AbstractPolyhedron{N}}

    out = apply(tr, X, HRepIntersection())
    return overapproximate(out, Hyperrectangle)
end

# the guard is the set union of half-spaces, the source invariant is polyhedral,
# and the target invariant is the set union of half-spaces
function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
                  X::AbstractPolytope{N},
                  method::BoxIntersection) where {N, VN,
                                                  GT<:UnionSetArray{N, HalfSpace{N, VN}},
                                                  IT⁻<:AbstractPolyhedron{N},
                                                  IT⁺<:Universe{N}}

    out = apply(tr, X, HRepIntersection())
    return overapproximate(out, Hyperrectangle)
end

#=
function _apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
                X::ConvexHullArray{N, PT},
                ::BoxEnclosure) where {N,
                                       PT<:AbstractPolyhedron{N},
                                       GT<:AbstractPolyhedron{N},
                                       VN,
                                       IT⁻<:UnionSetArray{N, LinearConstraint{N, VN}},
                                       IT⁺<:AbstractPolyhedron{N}}
    Xbox = overapproximate(X, Hyperrectangle)
    out = _apply(tr, Xbox, HRepIntersection())
    return overapproximate(out, Hyperrectangle)
end
=#

#=

function _apply(tr::DiscreteTransition{<:Universe, <:ZeroSet, GT, IT⁻, IT⁺},
                X::AbstractPolyhedron{N},
                ::BoxEnclosure) where {N,
                                       GT<:AbstractPolyhedron{N},
                                       VN,
                                       IT⁻<:UnionSetArray{N, LinearConstraint{N, VN}},
                                       IT⁺<:AbstractPolyhedron{N}}
    out = _apply(tr, X, HRepIntersection())
    return overapproximate(out, Hyperrectangle)
end

function _apply(tr::DiscreteTransition{<:Universe, <:ZeroSet, GT, IT⁻, IT⁺},
                X::ConvexHullArray{N, PT},
                ::BoxEnclosure) where {N,
                                       PT<:AbstractPolyhedron{N},
                                       GT<:AbstractPolyhedron{N},
                                       IT⁻<:AbstractPolyhedron{N},
                                       IT⁺<:AbstractPolyhedron{N}}

    Xbox = overapproximate(X, Hyperrectangle)
    out = _apply(tr, Xbox, HRepIntersection())
    return overapproximate(out, Hyperrectangle)
end
=#

# ===============================================
# Intersection method: TemplateHullIntersection
# ===============================================

# compute (R(X ∩ G ∩ I⁻) ⊕ W) ∩ I⁺ first computing the exact intersection
# of Y := X ∩ G ∩ I⁻ all of which are polyhedral, using their list of constraints
# then we compute (RY ⊕ W) ∩ I⁺ using again the given template, Z := [(RY ⊕ W) ∩ I⁺]_dirs
function apply(tr::DiscreteTransition{<:AbstractMatrix, <:LazySet, GT, IT⁻, IT⁺},
               X::PT,
               method::TemplateHullIntersection) where {N,
                                PT<:AbstractPolyhedron{N},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    # intersect X with the guard and the source invariant (NOTE can be cached)
    success, Y = _intersection(X, tr.G, tr.I⁻, HRepIntersection())
    !success && return EmptySet(dim(X))

    # compute the intersection [(RY ⊕ W) ∩ I⁺]_dirs using the template
    K = (tr.R * HPolyhedron(Y)) ⊕ tr.W # lazy affine map

    Z = _intersection(K, tr.I⁺, method)
    return Z
end

# compute (R(X ∩ G ∩ I⁻) ⊕ W) ∩ I⁺ first computing the exact intersection
# of Y := G ∩ I⁻ all of which are polyhedral, then compute K := [X ∩ Y]_dirs,
# finally we compute (RK ⊕ W) ∩ I⁺ using again the given template, Z := [(RK ⊕ W) ∩ I⁺]_dirs
function apply(tr::DiscreteTransition{<:AbstractMatrix, <:LazySet, GT, IT⁻, IT⁺},
               X::PT,
               method::TemplateHullIntersection) where {N,
                                PT<:LazySet{N},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    # intersect the guard and the source invariant
    success, Y = _intersection(tr.G, tr.I⁻, HRepIntersection())
    !success && return EmptySet(dim(X))

    K = _intersection(X, HPolyhedron(Y), method)
    isempty(K) && return EmptySet(dim(X))

    # compute Z := [(RK ⊕ W) ∩ I⁺]_dirs using the template
    Km = (tr.R * K) ⊕ tr.W # lazy affine map
    Z = _intersection(Km, tr.I⁺, method)
    return Z
end

# compute X ∩ G ∩ I⁻ ∩ I⁺ first computing the exact intersection
# of Y := G ∩ I⁻ ∩ I⁺ all of which are polyhedral, then compute K := [X ∩ Y]_dirs
function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
               X::PT,
               method::TemplateHullIntersection) where {N,
                                PT<:LazySet{N},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    # intersect the guard and the source invariant
    success, Y = _intersection(tr.G, tr.I⁻, tr.I⁺, HRepIntersection())
    !success && return EmptySet(dim(X))

    # compute the intersection K := [X ∩ Y]_dirs using the template
    K = _intersection(X, HPolyhedron(Y), method)
    return K
end

# ======================================
# Template hulls with set unions
#
# TODO: now X is a vector of reach-sets,
# we may dispatch on UnionSetArray later
# ==============================

# general case, requires constraints list of each Xi
function apply(tr::DiscreteTransition{<:AbstractMatrix, <:LazySet, GT, IT⁻, IT⁺},
               X::AbstractVector{RT},
               method::TemplateHullIntersection) where {N,
                                RT<:AbstractReachSet{N},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    @unpack R, W, G, I⁻, I⁺ = tr
    got_inv = LazySets.Arrays.isinvertible(R)

    m = length(X)
    Y = Vector{HPolytope{N, Vector{N}}}(undef, m)
    for i in 1:m
        Xi = set(X[i])

        # concrete Xi ∩ G ∩ I⁻
        success, Qi = _intersection(Xi, G, I⁻, HRepIntersection())
        !success && return EmptySet(dim(Xi))

        if got_inv
            # concretely compute R*Qi then lazily add W and lazily intersect with I⁺
            E = linear_map(R, HPolytope(Qi)) ⊕ W
            Y[i] = _intersection(E, I⁺, method)
        else
            Y[i] = _intersection(R * HPolytope(Qi) ⊕ W, I⁺, method)
        end
    end
    return overapproximate(ConvexHullArray(Y), method.dirs)
end

# case with R generic and W = 0
# requires constraints list of X
function apply(tr::DiscreteTransition{<:AbstractMatrix, <:ZeroSet, GT, IT⁻, IT⁺},
               X::AbstractVector{RT},
               method::TemplateHullIntersection) where {N,
                                RT<:AbstractReachSet{N},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    @unpack R, W, G, I⁻, I⁺ = tr
    got_inv = LazySets.Arrays.isinvertible(R)

    m = length(X)
    Y = Vector{HPolytope{N, Vector{N}}}(undef, m)
    for i in 1:m
        Xi = set(X[i])

        # concrete Xi ∩ G ∩ I⁻
        success, Qi = _intersection(Xi, G, I⁻, HRepIntersection())
        !success && return EmptySet(dim(Xi))

        # lazily or concretely compute R*Qi ∩ I⁺, depending on the invertibility of R
        if got_inv
            E = linear_map(R, HPolytope(Qi))
            success, Yi = _intersection(E, I⁺, HRepIntersection())
            !success && return EmptySet(dim(Xi))
            Y[i] = HPolytope(Yi)
        else
            Y[i] = _intersection(R * HPolytope(Qi), I⁺, method)
        end
    end
    return overapproximate(ConvexHullArray(Y), method.dirs)
end

# the R = Id and W = 0
function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
               X::AbstractVector{RT},
               method::TemplateHullIntersection) where {N,
                                RT<:AbstractReachSet{N},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    @unpack R, W, G, I⁻, I⁺ = tr

    # precompute L = G ∩ I⁻ ∩ I⁺
    success, L = _intersection(G, I⁻, I⁺, HRepIntersection())
    !success && return EmptySet(dim(X))

    m = length(X)
    Y = Vector{HPolytope{N, Vector{N}}}(undef, m)

    for i in 1:m
        Xi = set(X[i])

        # concrete intersection of Xi ∩ L
        success, Qi = _intersection(Xi, L, HRepIntersection())
        !success && return EmptySet(dim(Xi))
        Y[i] = HPolytope(Qi)
    end
    return overapproximate(ConvexHullArray(Y), method.dirs)
end

# the R = Id and W = 0
# and AbstractZonotope => use template directions for the intersection
function apply(tr::DiscreteTransition{<:IdentityMap, <:ZeroSet, GT, IT⁻, IT⁺},
               X::AbstractVector{RT},
               method::TemplateHullIntersection) where {N,
                                ZT<:AbstractZonotope{N},
                                RT<:ReachSet{N, ZT},
                                GT<:AbstractPolyhedron{N},
                                IT⁻<:AbstractPolyhedron{N},
                                IT⁺<:AbstractPolyhedron{N}}

    @unpack R, W, G, I⁻, I⁺ = tr

    # precompute L = G ∩ I⁻ ∩ I⁺
    success, L = _intersection(G, I⁻, I⁺, HRepIntersection())
    !success && return EmptySet(dim(X))
    L = HPolyhedron(L)

    m = length(X)
    Y = Vector{HPolytope{N, Vector{N}}}(undef, m)

    for i in 1:m
        Xi = set(X[i])

        # concrete intersection of Xi ∩ L
        Y[i] = _intersection(Xi, L, method)
    end

    return overapproximate(ConvexHullArray(Y), method.dirs)
end


text:
module Bullet3

# package code goes here
const file = joinpath(@__DIR__, "..", "deps", "deps.jl")
isfile(file) || error("Bullet3 not build correctly. Please run Pkg.build(\"Bullet3\")")

include(file)

function __init__()
    check_deps()
end

end # module

ext: jl


text:
function set_volume(nodes,cellxmax,cellymax)
    volume = zeros(cellxmax,cellymax)
    for i in 1:cellxmax
        for j in 1:cellymax
            vec_r1x = nodes[i+1,j+1,1] - nodes[i,j,1]
            vec_r1y = nodes[i+1,j+1,2] - nodes[i,j,2]
            vec_r2x = nodes[i,j+1,1] - nodes[i+1,j,1]
            vec_r2y = nodes[i,j+1,2] - nodes[i+1,j,2]

            volume[i,j] = abs(vec_r1x*vec_r2y - vec_r1y*vec_r2x) /2
        end
    end
    return volume
end

function setup_cell_flux_hat(Qcon,cellxmax,cellymax,cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus)
    cell_E_hat_plas  = zeros(cellxmax,cellymax,4)
    cell_E_hat_minus = zeros(cellxmax,cellymax,4)
    cell_F_hat_plas  = zeros(cellxmax,cellymax,4)
    cell_F_hat_minus = zeros(cellxmax,cellymax,4)

    for i in 1:cellxmax
        for j in 1:cellymax
            for l in 1:4
                for m in 1:4
                    # Ehat = Ahat * Qhat = Ahat * Q/J = 1/J*Ahat * Q
                    # 1/J*Ahatの計算をvecAを使ったため，ここではQconを掛けてる
                    cell_E_hat_plas[i,j,l]  += cell_Ahat_plas[i,j,l,m]*Qcon[i,j,m]
                    cell_E_hat_minus[i,j,l] += cell_Ahat_minus[i,j,l,m]*Qcon[i,j,m]
                    cell_F_hat_plas[i,j,l]  += cell_Bhat_plas[i,j,l,m]*Qcon[i,j,m]
                    cell_F_hat_minus[i,j,l] += cell_Bhat_minus[i,j,l,m]*Qcon[i,j,m]
                end
            end
        end
    end

    return cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus
end

function cal_jacobi(Qbase,Qcon,cellxmax,cellymax,specific_heat_ratio,vecAx,vecAy)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]

    q=u^2+v^2
    c=(rp/rho)^0.5
    theta=Z=U= A^xi_x * u   +   A^xi_y * v
    H=(e+p)/rho
    b1=q^2/2 * (r-1)/(c^2)
    b2=(r-1)/(c^2)

    xi_x_bar= xi_x / (xi_x^2+xi_y^2)^0.5
    xi_y_bar= xi_y / (xi_x^2+xi_y^2)^0.5
    U_bar=Z_bar=theta_bar=

    A = Rleft * Lambda * Rright
    """

    cell_Ahat_plas  = zeros(cellxmax,cellymax,4,4)
    cell_Ahat_minus = zeros(cellxmax,cellymax,4,4)
    cell_Bhat_plas  = zeros(cellxmax,cellymax,4,4)
    cell_Bhat_minus = zeros(cellxmax,cellymax,4,4)
    r=specific_heat_ratio

    for i in 1:cellxmax
        for j in 1:cellymax
            q=(Qbase[i,j,2]^2+Qbase[i,j,3]^2)^0.5
            c=0
            try c=(r*Qbase[i,j,4]/Qbase[i,j,1])^0.5
            catch
                println("\n"*string(i)*","*string(j)*" Qbase error")
                println(Qbase[i,j,:])
                println("\n")
                throw(UndefVarError(:x))
            end

            H  = (Qcon[i,j,4]+Qbase[i,j,4])/Qbase[i,j,1]
            b1 = q^2/2 * (r-1)/(c^2)
            b2 = (r-1)/(c^2)
            u  = Qbase[i,j,2]
            v  = Qbase[i,j,3]

            k_x = 0.5*(vecAx[i,j,1]+vecAx[i+1,j,1])
            k_y = 0.5*(vecAx[i,j,2]+vecAx[i+1,j,2])

            Rleft_E,Lambda_E,Rright_E = Eigenvalue_vector(q,c,H,b1,b2,u,v,k_x,k_y)
            A_plas, A_minus = RLmbdaR(Rleft_E,Lambda_E,Rright_E)

            for l in 1:4
                for m in 1:4
                    cell_Ahat_plas[i,j,l,m]  = A_plas[l,m]
                    cell_Ahat_minus[i,j,l,m] = A_minus[l,m]
                end
            end

            k_x = 0.5*(vecAy[i,j,1]+vecAy[i,j+1,1])
            k_y = 0.5*(vecAy[i,j,2]+vecAy[i,j+1,2])

            Rleft_E,Lambda_E,Rright_E = Eigenvalue_vector(q,c,H,b1,b2,u,v,k_x,k_y)
            A_plas, A_minus = RLmbdaR(Rleft_E,Lambda_E,Rright_E)

            for l in 1:4
                for m in 1:4
                    cell_Bhat_plas[i,j,l,m]  = A_plas[l,m]
                    cell_Bhat_minus[i,j,l,m] = A_minus[l,m]
                end
            end
        end
    end

    return cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus
end

function Eigenvalue_vector(q,c,H,b1,b2,u,v,k_x,k_y)
    Z= k_x*u + k_y*v
    k_x_bar=k_x / (k_x^2+k_y^2)^0.5
    k_y_bar=k_y / (k_x^2+k_y^2)^0.5
    Z_bar=Z / (k_x^2+k_y^2)^0.5

    Lambda = [Z-c*(k_x^2+k_y^2)^0.5
              Z
              Z+c*(k_x^2+k_y^2)^0.5
              Z]

    Rleft = [[1 1 1 0]
             [u-k_x_bar*c u u+k_x_bar*c -k_y_bar]
             [v-k_y_bar*c v v+k_y_bar*c k_x_bar]
             [H-c*Z_bar q^2/2 H+c*Z_bar -(k_y_bar*u-k_x_bar*v)]]

    Rright = [[(b1+Z_bar/c)/2 -(k_x_bar/c+b2*u)/2 -(k_y_bar/c+b2*v)/2 b2/2]
              [1-b1 b2*u b2*v -b2]
              [(b1-Z_bar/c)/2 (k_x_bar/c-b2*u)/2 (k_y_bar/c-b2*v)/2 b2/2]
              [k_y_bar*u-k_x_bar*v -k_y_bar k_x_bar 0]]

    return Rleft,Lambda,Rright
end

function RLmbdaR(R,Lam,Rm)
    Lamp=zeros(4,4)
    Lamm=zeros(4,4)
    for i in 1:length(Lam)
        Lamp[i,i]=(Lam[i]+abs(Lam[i]))/2
        Lamm[i,i]=(Lam[i]-abs(Lam[i]))/2
    end
    Ap=nn_inner_product(nn_inner_product(R,Lamp),Rm)
    Am=nn_inner_product(nn_inner_product(R,Lamm),Rm)
    return Ap,Am
end

function nn_inner_product(a,b)
    temp=zeros(size(a)[1],size(a)[1])
    for i in 1:size(a)[1]
        for j in 1:size(a)[1]
            for k in 1:size(a)[1]
                temp[i,j] += a[i,k]*b[k,j]
            end
        end
    end
    return temp
end

function FVS(cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus,cellxmax,cellymax)
    E_hat = zeros(cellxmax+1,cellymax,4)
    F_hat = zeros(cellxmax,cellymax+1,4)

    Threads.@threads for i in 2:cellxmax+1-1
        for j in 2:cellymax-1
            for k in 1:4
                E_hat[i,j,k] = cell_E_hat_plas[i-1,j,k] + cell_E_hat_minus[i,j,k]
            end
        end
    end
    Threads.@threads for i in 2:cellxmax-1
        for j in 2:cellymax+1-1
            for k in 1:4
                F_hat[i,j,k] = cell_F_hat_plas[i,j-1,k] + cell_F_hat_minus[i,j,k]
            end
        end
    end

    return E_hat,F_hat
end

function setup_RHS(cellxmax,cellymax,E_hat,F_hat,Qcon_hat)
    RHS = zeros(cellxmax,cellymax,4)

    Threads.@threads for i in 2:cellxmax-1
        for j in 2:cellymax-1
            for k in 1:4
                RHS[i,j,k] = -(E_hat[i+1,j,k]-E_hat[i,j,k]+F_hat[i,j+1,k]-F_hat[i,j,k])
            end
        end
    end
    return RHS
end

function set_RHS_inner(cellxmax,cellymax,dt,Qcon_hat_n,Qcon_hat_m1,E_hat,F_hat)
    RHS = zeros(cellxmax,cellymax,4)

    Threads.@threads for i in 2:cellxmax-1
        for j in 2:cellymax-1
            for k in 1:4
                RHS[i,j,k] = -(Qcon_hat_m1[i,j,k]-Qcon_hat_n[i,j,k]) -dt*(E_hat[i+1,j,k]-E_hat[i,j,k]+F_hat[i,j+1,k]-F_hat[i,j,k])
            end
        end
    end
    return RHS
end

function setup_Qcon_hat(Qcon,cellxmax,cellymax,volume)
    Qcon_hat=ones(cellxmax,cellymax,4)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:4
                Qcon_hat[i,j,k] = Qcon[i,j,k]*volume[i,j]
            end
        end
    end
    return Qcon_hat
end

function Qhat_to_Q(Qcon_hat,cellxmax,cellymax,volume)
    Qcon = ones(cellxmax,cellymax,4)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:4
                Qcon[i,j,k] = Qcon_hat[i,j,k]/volume[i,j]
            end
        end
    end
    return Qcon
end


text:
@testset "Time Function Macro Tests" begin

  @test isdefined(SetAndForget, Symbol("@time_function")) == true

end


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

using Test

A = Int64[1, 2, 3, 4]
B = Complex{Int64}[5+6im, 7+8im, 9+10im]
# getindex
@test reinterpret(Complex{Int64}, A) == [1 + 2im, 3 + 4im]
@test reinterpret(Float64, A) == reinterpret.(Float64, A)

@test reinterpret(NTuple{3, Int64}, B) == [(5,6,7),(8,9,10)]

# setindex
let Ac = copy(A), Bc = copy(B)
    reinterpret(Complex{Int64}, Ac)[2] = -1 - 2im
    @test Ac == [1, 2, -1, -2]
    reinterpret(NTuple{3, Int64}, Bc)[2] = (4,5,6)
    @test Bc == Complex{Int64}[5+6im, 7+4im, 5+6im]
    reinterpret(NTuple{3, Int64}, Bc)[1] = (1,2,3)
    @test Bc == Complex{Int64}[1+2im, 3+4im, 5+6im]

    A1 = reinterpret(Float64, A)
    A2 = reinterpret(Complex{Float64}, A)
    A1[1] = 1.0
    @test real(A2[1]) == 1.0
end

# same-size reinterpret where one of the types is non-primitive
let a = NTuple{4,UInt8}[(0x01,0x02,0x03,0x04)]
    @test reinterpret(Float32, a)[1] == reinterpret(Float32, 0x04030201)
    reinterpret(Float32, a)[1] = 2.0
    @test reinterpret(Float32, a)[1] == 2.0
end

# ensure that reinterpret arrays aren't erroneously classified as strided
let A = reshape(1:20, 5, 4)
    V = view(A, :, :)
    R = reinterpret(Int32, V)
    R2 = reinterpret(Int32, A)
    @test !(R isa StridedArray)
    @test !(R2 isa StridedArray)
    @test R * ones(4, 5) == R2 * ones(4,5) == copy(R) * ones(4,5) == copy(R2) * ones(4,5)
end

# but ensure that strided views of strided reinterpret arrays are still strided
let A = collect(reshape(1:20, 5, 4))
    R = reinterpret(Int32, A)
    @test R isa StridedArray
    @test view(R, :, :) isa StridedArray
    @test reshape(R, :) isa StridedArray
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function act_2_solution()
  cur_T = 17.8

  cur_reactor = act_2_deck(cur_T)

  cur_reactor.n_bar = 0.86
  cur_reactor.I_P = 13.98

  cur_reactor.R_0 = 9.75
  cur_reactor.B_0 = 8.75

  cur_reactor.P_F = 2637

  cur_reactor.f_BS = 0.77
  cur_reactor.f_CD = 1.00 - cur_reactor.f_BS
  cur_reactor.f_IN = 0.00

  cur_reactor.beta_N = cur_reactor.max_beta_N
  cur_reactor.q_95 = 8.0
  cur_reactor.P_W = 1.46

  cur_reactor.volume = 2209

  cur_reactor
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_AuthorizationListResult.jl")
include("model_AuthorizationPropertiesFormat.jl")
include("model_CircuitConnectionStatus.jl")
include("model_CloudError.jl")
include("model_CloudErrorBody.jl")
include("model_ExpressRouteCircuit.jl")
include("model_ExpressRouteCircuitArpTable.jl")
include("model_ExpressRouteCircuitAuthorization.jl")
include("model_ExpressRouteCircuitConnection.jl")
include("model_ExpressRouteCircuitConnectionListResult.jl")
include("model_ExpressRouteCircuitConnectionPropertiesFormat.jl")
include("model_ExpressRouteCircuitListResult.jl")
include("model_ExpressRouteCircuitPeering.jl")
include("model_ExpressRouteCircuitPeeringConfig.jl")
include("model_ExpressRouteCircuitPeeringListResult.jl")
include("model_ExpressRouteCircuitPeeringPropertiesFormat.jl")
include("model_ExpressRouteCircuitPropertiesFormat.jl")
include("model_ExpressRouteCircuitRoutesTable.jl")
include("model_ExpressRouteCircuitRoutesTableSummary.jl")
include("model_ExpressRouteCircuitServiceProviderProperties.jl")
include("model_ExpressRouteCircuitSku.jl")
include("model_ExpressRouteCircuitStats.jl")
include("model_ExpressRouteCircuitsArpTableListResult.jl")
include("model_ExpressRouteCircuitsRoutesTableListResult.jl")
include("model_ExpressRouteCircuitsRoutesTableSummaryListResult.jl")
include("model_ExpressRouteConnectionId.jl")
include("model_ExpressRoutePeeringState.jl")
include("model_ExpressRoutePeeringType.jl")
include("model_ExpressRouteServiceProvider.jl")
include("model_ExpressRouteServiceProviderBandwidthsOffered.jl")
include("model_ExpressRouteServiceProviderListResult.jl")
include("model_ExpressRouteServiceProviderPropertiesFormat.jl")
include("model_Ipv6CircuitConnectionConfig.jl")
include("model_Ipv6ExpressRouteCircuitPeeringConfig.jl")
include("model_PeerExpressRouteCircuitConnection.jl")
include("model_PeerExpressRouteCircuitConnectionListResult.jl")
include("model_PeerExpressRouteCircuitConnectionPropertiesFormat.jl")
include("model_ProvisioningState.jl")
include("model_Resource.jl")
include("model_ServiceProviderProvisioningState.jl")
include("model_SubResource.jl")
include("model_TagsObject.jl")


text:
"""
Quadratic Polynomial line search
"""
struct QuadraticState{OBJ,OPT,T} <: LinesearchState where {OBJ <: AbstractObjective, OPT <: Options, T <: Number}
    objective::OBJ
    config::OPT

    α₀::T
    σ₀::T
    σ₁::T
    ϵ::T

    function QuadraticState(objective; config = Options(),
                    α₀::T = DEFAULT_ARMIJO_α₀,
                    σ₀::T = DEFAULT_ARMIJO_σ₀,
                    σ₁::T = DEFAULT_ARMIJO_σ₁,
                    ϵ::T = DEFAULT_WOLFE_ϵ) where {T}
        new{typeof(objective), typeof(config), T}(objective, config, α₀, σ₀, σ₁, ϵ)
    end
end

function QuadraticState(F::Callable, x::AbstractVector; D = nothing, kwargs...)
    objective = MultivariateObjective(F, D, x)
    QuadraticState(objective; kwargs...)
end

function QuadraticState(F::Callable, x::Number; D = nothing, kwargs...)
    objective = UnivariateObjective(F, D, x)
    QuadraticState(objective; kwargs...)
end

Base.show(io::IO, ls::QuadraticState) = print(io, "Polynomial quadratic")

LinesearchState(algorithm::Quadratic, objective; kwargs...) = QuadraticState(objective; kwargs...)

function (ls::QuadraticState{<:UnivariateObjective})()
    T = typeof(ls.α₀)
    local α::T = ls.α₀
    local αₜ::T
    local y₁::T
    local p₀::T
    local p₁::T
    local p₂::T

    # determine constant coefficients of polynomial p(α) = p₀ + p₁α + p₂α²
    # p₀ = f(x₀)    # value of initial solution
    # p₁ = f'(x₀)   # derivative of initial solution
    y₀ = value!(ls.objective, zero(α))
    d₀ = derivative!(ls.objective, zero(α))
    p₀ = y₀
    p₁ = d₀

    for _ in 1:ls.config.max_iterations
        # compute value of new solution
        y₁ = value!(ls.objective, α)

        if y₁ ≥ y₀ + ls.ϵ * α * d₀
            # determine nonconstant coefficient of polynomial p(α) = p₀ + p₁α + p₂α²
            p₂ = (y₁^2 - p₀ - p₁*α) / α^2

            # compute minimum αₜ of p(α)
            αₜ = - p₁ / (2p₂)

            if αₜ < ls.σ₀ * α
                α = ls.σ₀ * α
            elseif αₜ > ls.σ₁ * α
                α = ls.σ₁ * α
            else
                α = αₜ
            end
        else
            break
        end
    end

    return α
end

quadratic(f, x; kwargs...) = QuadraticState(f, x; kwargs...)()


text:
module Plot3d

import Winston.output_surface
if output_surface == :tk
    eval(Expr(:toplevel, Expr(:using, :Tk)))
elseif output_surface == :gtk
    eval(Expr(:toplevel, Expr(:using, :Gtk)))
end
using Color
using Base.Graphics

export plot3d, surf
#export demo_sombrero, demo_sphere

type Canvas3D
    win::Canvas
    ctm::Matrix{Float64}
    sctm::Matrix{Float64}
    lastx::Int
    lasty::Int
    GW::Float64
    GH::Float64
    scalem::Vector{Float64}
    boxv::Matrix{Float64}
    center::Vector{Float64}
    wincenter::Vector{Float64}
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
    zmin::Float64
    zmax::Float64
    colorbg::RGB
    colorcube::RGB
    models_motion::Vector{Any}
    models_release::Vector{Any}

    function Canvas3D(win; xmin=0, xmax=width(win)-1, ymin=0, ymax=height(win)-1,
                      zmin=-10, zmax=10, colorbg=RGB(1,1,1), colorcube=RGB(0,0,0))
        this = new(win)
        this.xmin = xmin; this.xmax = xmax
        this.ymin = ymin; this.ymax = ymax
        this.zmin = zmin; this.zmax = zmax
        this.ctm = eye(3)
        this.colorbg = colorbg
        this.colorcube = colorcube
        this.models_motion = Any[]
        this.models_release = Any[]

        win.mouse.button1press = (c,x,y)->canvas3d_mousedown(this,x,y)
        win.mouse.button1motion = (c,x,y)->canvas3d_button1motion(this,x,y)
        win.mouse.button1release = (c,x,y)->canvas3d_button1release(this,x,y)
        win.resize = w->configure(this)
        win.draw = c->draw(getgc(this.win), this, false)
        this
    end
end

Base.show(io::IO, c3::Canvas3D) = print(io, "Canvas3D")

cube_verts(x0, x1, y0, y1, z0, z1) = [x0 x1 x1 x0 x0 x1 x1 x0
                                      y0 y0 y1 y1 y0 y0 y1 y1
                                      z0 z0 z0 z0 z1 z1 z1 z1]

function configure(this::Canvas3D)
    WW, WH = width(this.win), height(this.win)
    scalem = [ 0.5774 * WW/(this.xmax - this.xmin),
              -0.5774 * WH/(this.ymax - this.ymin),
               0.5774 * WW/(this.xmax - this.xmin)]
    this.scalem = scalem
    this.sctm = scale(this.ctm, scalem)

    this.center = [(this.xmax+this.xmin)/2,
                   (this.ymax+this.ymin)/2,
                   (this.zmax+this.zmin)/2]
    this.boxv = cube_verts(this.xmin, this.xmax, this.ymin, this.ymax,
                           this.zmin, this.zmax)
    this.GW = WW/2
    this.GH = WH/2
    this.wincenter = [this.GW, this.GH, 0]
    this
end

function project(this::Canvas3D, v::AbstractVector)
    this.sctm * (v-this.center) + this.wincenter
end

function project(this::Canvas3D, v::AbstractMatrix)
    broadcast(+, this.sctm * broadcast(-,v,this.center), this.wincenter)
end

const cube4sides = {[1,5,6,2], [2,6,7,3], [3,7,8,4], [4,8,5,1]}

function draw(gc, this::Canvas3D, motion::Bool)
    bv = similar(this.boxv)
    for i=1:size(bv,2)
        bv[:,i] = project(this, this.boxv[:,i])
    end

    edges = sort(cube4sides, by=p->mean(bv[3,:][p]))

    set_source(gc, this.colorbg)
    paint(gc)

    set_source(gc, this.colorcube)
    set_line_width(gc, 0.6)

    polygon(gc, bv, edges[1])
    polygon(gc, bv, edges[2])
    stroke(gc)

    # draw contents here
    if motion
        for m in this.models_motion
            draw(gc, this, m)
        end
    else
        for m in this.models_release
            draw(gc, this, m)
        end
    end

    polygon(gc, bv, edges[3])
    polygon(gc, bv, edges[4])
    stroke(gc)
end

function canvas3d_mousedown(this::Canvas3D, x, y)
    this.lastx = x
    this.lasty = y
end

# construct a rotation by φ around the given vector
function rotation(x,y,z,φ)
    a = 0.5φ
    s = sin(a)/sqrt(x*x+y*y+z*z)
    cos(a), [x*s, y*s, z*s]
end

# apply quaternion to vector
qrotate(w, qv, v) = v - 2*cross(cross(qv, v) - w.*v, qv)

# project window point x,y onto a sphere of radius r and center cx,cy
function sphereproject(r, cx, cy, x, y)
    qx = x-cx; qy = y-cy
    sph = r^2 - qx^2 - qy^2
    if sph < 0
        [r*qx, r*qy, 0]
    else
        [qx, qy, sqrt(sph)]
    end
end

# arcball - translate mouse motion into a rotation, as a 3x3 matrix
# x, y         - new mouse position
# lastx, lasty - old mouse position
# W, H         - window width / 2, height / 2
# ctm          - the current transform matrix (must be orthogonal)
function arcball(x, y, lastx, lasty, W, H, ctm)
    r = max(W,H)
    q0 = sphereproject(r, W, H, lastx, lasty)
    q1 = sphereproject(r, W, H, x, y)

    if q0 == q1
        return ctm
    end

    ictm = ctm'
    rx, ry, rz = ictm * cross(q0,q1)
    w, v = rotation(rx, ry, rz,
                    2acos(dot(q0,q1) / (norm(q0)*norm(q1))))

    (any(isnan,v)||isnan(w)) && return ctm

    xv = qrotate(w, v, ictm[:,1]); xv /= norm(xv)
    yv = qrotate(w, v, ictm[:,2]); yv /= norm(yv)
    zv = qrotate(w, v, ictm[:,3]); zv /= norm(zv)
    [xv[1] xv[2] xv[3]
     yv[1] yv[2] yv[3]
     zv[1] zv[2] zv[3]]
end

function canvas3d_mouseupdate(this::Canvas3D, x, y)
    this.ctm = arcball(x, y, this.lastx, this.lasty, this.GW, this.GH,
                       this.ctm)
    this.sctm = scale(this.ctm, this.scalem)
    this.lastx = x
    this.lasty = y
end

function canvas3d_button1motion(this::Canvas3D, x, y)
    canvas3d_mouseupdate(this, x, y)
    draw(getgc(this.win), this, true)
    reveal(this.win)
end

function canvas3d_button1release(this::Canvas3D, x, y)
    canvas3d_mouseupdate(this, x, y)
    draw(getgc(this.win), this, false)
    reveal(this.win)
end

# connectivity of m x n grid
function grid_polygons(m,n)
    E = Vector{Int}[]
    for k in 0:n-2, j in 0:m-2
        i = k*m+j+1
        push!(E, [i, i+1, i+m+1, i+m])  # quads
        #push!(E, [i, i+1, i+m+1])      # triangles
        #push!(E, [i, i+m+1, i+m])
    end
    E
end

# evaluate a surface over ranges of u,v parameters, giving a 3xN vertex matrix
function evalsurface(xf, yf, zf, ur, vr)
    X = Float64[ xf(u,v) for u in ur, v in vr ]
    Y = Float64[ yf(u,v) for u in ur, v in vr ]
    Z = Float64[ zf(u,v) for u in ur, v in vr ]
    N = length(ur)*length(vr)
    [reshape(X, 1, N)
     reshape(Y, 1, N)
     reshape(Z, 1, N)]
end

type Polygons3D
    V::Matrix{Float64}
    P::Vector{Vector{Int}}
    colors

    Polygons3D(V, P, colors=[ RGB(1,1,1) for i=1:length(P) ]) =
        new(V, P, colors)
    Polygons3D(V, P, coloring::Function) =
        Polygons3D(V, P,
                   [ coloring(V[1,p[1]], V[2,p[1]], V[3,p[1]]) for p in P ])
end

function draw(gc, c::Canvas3D, this::Polygons3D)
    v = project(c, this.V)
    z_ord = sortperm(this.P, by=p->v[3,p[1]])  # z sort
    set_line_width(gc, 0.5)
    for n in z_ord
        p = this.P[n]
        valid = true
        for i in p
            if isnan(v[1,i]) || isnan(v[2,i])
                valid = false; break
            end
        end
        if valid
            polygon(gc, v, p)
            set_source(gc, this.colors[n])
            fill_preserve(gc)
            set_source_rgb(gc, 0, 0, 0)
            stroke(gc)
        end
    end
end

function surf(xf::Function, yf::Function, zf::Function, ur, vr; coloring=false)
    V = evalsurface(xf, yf, zf, ur, vr)
    P = grid_polygons(length(ur), length(vr))
    if coloring === false
        Polygons3D(V, P)
    else
        Polygons3D(V, P, coloring)
    end
end

function surf(X::AbstractMatrix, Y::AbstractMatrix, Z::AbstractMatrix)
    n = length(X)
    Polygons3D([reshape(X,1,n)
                reshape(Y,1,n)
                reshape(Z,1,n)],
               grid_polygons(size(X,1),size(X,2)))
end

function plot3d(o::Polygons3D)
    w = Window("3d plot", 320, 320)
    if output_surface == :tk
        c = Canvas(w)
        pack(c, expand = true, fill = "both")
    else
        c = Canvas()
        push!(w, c)
    end

    xmin = minimum(o.V[1,:]); xmax = maximum(o.V[1,:])
    ymin = minimum(o.V[2,:]); ymax = maximum(o.V[2,:])
    zmin = minimum(o.V[3,:]); zmax = maximum(o.V[3,:])

    c3d = Canvas3D(c, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
                   zmin=zmin, zmax=zmax)
    push!(c3d.models_motion, o)
    push!(c3d.models_release, o)
    c3d
end

demo_sphere() =
    plot3d(surf((u,v)->cos(v*pi)*sin(u),
                (u,v)->-cos(v*pi)*cos(u),
                (u,v)->sin(v*pi),
                0:(2pi/29):2pi, -.5:(1/17):.5,
                coloring = (x,y,z)->RGB((x-y+1)/3+.3,
                                        (z-y+1)/3+.3,
                                        z/1.5+.3)))

demo_sombrero() =
    plot3d(surf((u,v)->u,
                (u,v)->sin(hypot(u,v))/hypot(u,v),
                (u,v)->v,
                -8:.53:8, -8:.53:8))

end


text:
using ImageCore, ImageFiltering, ImageTransformations, CoordinateTransformations, Interpolations, OffsetArrays, StaticArrays, ColorTypes, FixedPointNumbers, TestImages, IdentityRanges, MappedArrays, ComputationalResources, MLDataPattern, ImageInTerminal, ReferenceTests, Base.Test

# check for ambiguities
refambs = detect_ambiguities(ImageTransformations, Base, Core)
using Augmentor
ambs = detect_ambiguities(Augmentor, ImageTransformations, Base, Core)
# The 1 is from plain_indices with a Tuple{} (so its spurious)
@test length(setdiff(ambs, refambs)) == 1

str_show(obj) = @io2str Base.show(::IO, obj)
str_showcompact(obj) = @io2str Base.showcompact(::IO, obj)
str_showconst(obj) = @io2str Augmentor.showconstruction(::IO, obj)

camera = testimage("cameraman")
cameras = similar(camera, size(camera)..., 2)
copy!(view(cameras,:,:,1), camera)
copy!(view(cameras,:,:,2), camera)
square = Gray{N0f8}[0.1 0.2 0.3; 0.4 0.5 0.6; 0.7 0.6 0.9]
square2 = rand(Gray{N0f8}, 4, 4)
rect = Gray{N0f8}[0.1 0.2 0.3; 0.4 0.5 0.6]
checkers = Gray{N0f8}[1 0 1 0 1; 0 1 0 1 0; 1 0 1 0 1]
rgb_rect = rand(RGB{N0f8}, 2, 3)

tests = [
    "tst_utils.jl",
    "operations/tst_channels.jl",
    "operations/tst_dims.jl",
    "operations/tst_convert.jl",
    "operations/tst_mapfun.jl",
    "operations/tst_noop.jl",
    "operations/tst_cache.jl",
    "operations/tst_rotation.jl",
    "operations/tst_shear.jl",
    "operations/tst_flip.jl",
    "operations/tst_crop.jl",
    "operations/tst_resize.jl",
    "operations/tst_scale.jl",
    "operations/tst_zoom.jl",
    "operations/tst_distortions.jl",
    "operations/tst_either.jl",
    "tst_operations.jl",
    "tst_pipeline.jl",
    "tst_augment.jl",
    "tst_augmentbatch.jl",
    "tst_distortedview.jl",
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
function Base.select(f, db::DB, table::String, sel::String="*")
    select_sql = "SELECT $sel FROM $table"
    # extract body of f as Expr object
    ast = Base.uncompressed_ast(f.code)
    e = ast.args[3]
    # extract arguments from call to `where`
    # NOTE: `where` serves only as a parsing landmark; it has no implementation
    clauses = parse_where(e)
    where_sql = " WHERE"
    for (i, clause) in enumerate(clauses)
        i > 1 ? where_sql = where_sql * " AND $clause" :
                where_sql = where_sql * " $clause"
    end
    res_sql = select_sql * where_sql
    return Cursor(Stmt(db, res_sql))
end

parse_where(e) = nothing
function parse_where(e::Expr)
    e = de_global(e) # preprocess by converting GlobalRefs to Symbols
    if (e.head == :call) & (e.args[1] == :where)
        res = e.args[2:end]
    else
        for arg in e.args
            res = parse_where(arg)
        end
    end
    return res
end

de_global(e) = e
de_global(g::GlobalRef) = g.name
function de_global(e::Expr)
    return Expr(e.head, [ de_global(arg) for arg in e.args ]...)
end


text:
# Generate documentation with this command:
# (cd docs && julia --color=yes make.jl)

push!(LOAD_PATH, "..")

using Documenter
using ADIOS2

makedocs(; sitename="ADIOS2", format=Documenter.HTML(), modules=[ADIOS2])

deploydocs(; repo="github.com/eschnett/ADIOS2.jl.git", devbranch="main",
           push_preview=true)


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0201 - Data Wrangling

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will hower be lost)
## Pkg.resolve()
Pkg.instantiate()
using Random
Random.seed!(123)

# ## Introduction

# This segment will be mostly based on [DataFrames](https://github.com/JuliaData/DataFrames.jl) and related packages

# !!! info DataFrames vs Matrix
#     DataFrames are popular format for **in-memory tabular data**. Their main advantages over Arrays are that they can efficiently store different types of data on each column (indeed each column is a wrapper over an `Array{T,1}` where `T` is specific to each column) and, thanks also to their named columns, provide convenient API for data operations, like indexing, querying , joining, split-apply-combine, etc.

# !!! info
#     In most circunstances we can refer to dataframe columns either by using their name as a string, e.g. `"Region"`, or as symbol, e.g. `:Region`. In the rest of the segment we will use the strings approach.

# ## Data import

# Our example: Forest volumes and area by country and year
# Source: Eurostat; units: forarea: Milion hectars, forvol: Milion cubic metres

# ### Built-in solution: CSV --> Matrix
using DelimitedFiles  # in the stdlib
data = convert(Array{Float64,2},readdlm("data.csv",';')[2:end,3:end]) # Skip the first 1 row and the first 2 columns - out is a Matrix

# ### CSV.jl: {CSV file, hardwritten data} --> DataFrame
using CSV, DataFrames
data = CSV.read("data.csv",DataFrame) # source, sink, kword options
data = CSV.read("data.csv",NamedTuple)

data = CSV.read(IOBuffer("""
Country	Year	forarea	forvol
Germany	2000	11.354	3381
France	2000	15.288	2254.28
Italy	2000	8.36925	1058.71
Sweden	2000	28.163	3184.67
Germany	2020	11.419	3663
France	2020	17.253	3055.83
Italy	2020	9.56613	1424.4
Sweden	2020	27.98	3653.91
"""), DataFrame, copycols=true)

# Some common CSV.jl options: `delim` (use `'\t'` for tab delimited files), `quotechar`, `openquotechar`, `closequotechar`, `escapechar`, `missingstring`, `dateformat`, `append`, `writeheader`, `header`, `newline`, `quotestrings`, `decimal`, `header`, `normalizenames`, `datarow`, `skipto`, `footerskip`, `limit`, `transpose`, `comment`, `use_mmap`, `type`, `types` (e.g. `types=Dict("fieldFoo" => Union{Missing,Int64})`), `typemap`, `pool`, `categorical`, `strict`, `silencewarnings`, `ignorerepeated`

# ### XLSX.jl: xlsx -> {Matrix, DataFrame}
using XLSX
sheetNames = XLSX.sheetnames(XLSX.readxlsx("data.xlsx"))
data = XLSX.readxlsx("data.xlsx")["Sheet1"]["A1:D9"]
data = XLSX.readxlsx("data.xlsx")["Sheet1"][:]
data = XLSX.readdata("data.xlsx", "Sheet1", "A1:D9")
XLSX.readtable("data.xlsx", "Sheet1") # tuple vector of (data) vectors, vector of symbols, usable as DF constructor
data = DataFrame(XLSX.readtable("data.xlsx", "Sheet1")...)

# ### OdsIO.jl: ods -> {Matrix, DataFrame}
using OdsIO
ods_read("data.ods";sheetName="Sheet1",retType="DataFrame")
ods_read("data.ods";sheetName="Sheet1",retType="Matrix",range=[(2,3),(9,4)]) # [(tlr,tlc),(brr,brc)]

# ### HTTP.jl: from internet
import HTTP
using Pipe, ZipFile, Tar
urlData = "https://github.com/sylvaticus/IntroSPMLJuliaCourse/raw/main/lessonsSources/02_-_JULIA2_-_Scientific_programming_with_Julia/data.csv"
urlDataZ = "https://github.com/sylvaticus/IntroSPMLJuliaCourse/raw/main/lessonsSources/02_-_JULIA2_-_Scientific_programming_with_Julia/data.zip"

data = @pipe HTTP.get(urlData).body        |>
             replace!(_, UInt8(';') => UInt8(' ')) |>  # if we need to do modifications to the file before importing
             CSV.File(_, delim=' ')                |>
             DataFrame;

# ### ZipFile.jl : from a single zipped csv file..
# ...on disk...
data = @pipe ZipFile.Reader("data.zip").files[1] |>
             CSV.File(read(_), delim=';') |>
             DataFrame

# ...or on internet:
data = @pipe HTTP.get(urlDataZ).body       |>
             IOBuffer(_)                      |>
             ZipFile.Reader(_).files[1]    |>
             CSV.File(read(_), delim=';')  |>
             DataFrame

# ### DataFrame constructor
## named tuple of (colName,colData)..
data = DataFrame(
    country = ["Germany", "France", "Italy", "Sweden", "Germany", "France", "Italy", "Sweden"],
    year    = [2000,2000,2000,2000,2020,2020,2020,2020],
    forarea = [11.354, 15.288, 8.36925, 28.163, 11.419, 17.253, 9.56613, 27.98],
    forvol  = [3381, 2254.28, 1058.71, 3184.67, 3663, 3055.83, 1424.4, 3653.91]
)

# ### Matrix -> DataFrame

# Headers and data separated:

M = ["Germany"	2000	11.354	3381
     "France"	2000	15.288	2254.28
     "Italy"	2000	8.36925	1058.71
     "Sweden"	2000	28.163	3184.67
     "Germany"	2020	11.419	3663
     "France"	2020	17.253	3055.83
     "Italy"	2020	9.56613	1424.4
     "Sweden"	2020	27.98	3653.91]
headers = ["Country", "Year", "forarea", "forvol"]
## array of colData arrays, array of headers
data = DataFrame([[M[:,i]...] for i in 1:size(M,2)], Symbol.(headers))

# Headers on the first row of the matrix:
M = ["Country"	"Year"	"forarea" "forvol"
     "Germany"	2000	11.354	  3381
     "France"	2000	15.288	  2254.28
     "Italy"	2000	8.36925	  1058.71
     "Sweden"	2000	28.163	  3184.67
     "Germany"	2020	11.419	  3663
     "France"	2020	17.253	  3055.83
     "Italy"	2020	9.56613	  1424.4
     "Sweden"	2020	27.98	  3653.91]
data = DataFrame([[M[2:end,i]...] for i in 1:size(M,2)], Symbol.(M[1,:])) # note the autorecognision of col types

# ## Getting insights on the data

# In VSCode, we can also use the workpanel for a nice sortable tabular view of a df

show(data,allrows=true,allcols=true)
first(data,6)
last(data, 6)
describe(data)
nR,nC = size(data)
names(data)
for r in eachrow(data)
    println(r) # note is it a "DataFrameRow"
end
for c in eachcol(data)
    println(c) # an array
    println(nonmissingtype(eltype(c)))
end

# ## Selection and querying data

# ### Column(s) selection
# In general we can select: (a) by name (strings or symbols) or by position; (b) a single or a vector of columns, (c) copying or making a view (a reference without copying) of the underlying data
data[:,["Country","Year"]] # copy, using strings
data[!,[:Country,:Year]]   # view, using symbols
data.Year                  # equiv. to `data[!,Year]`
data[:,1]                  # also `data[!,1]`
data[:,Not(["Year"])]

# ### Row(s) selection
data[1,:]    # DataFrameRow
data[1:3,:]  # DataFrame
# Note rows have no title names as colums do.

# ### Cell(s) selection
data[2,[2,4]]
data[2,"forarea"]

# Note that the returned selection is:
# * an `Array{T,1}` if it is a single column;
# * a `DataFrameRow` (similar in behaviour to a `DataFrame`) if a single row;
# * `T` if a single cell;
# * an other `DataFrame` otherwise.

# ### Boolean selection
# Both rows and column of a DataFrame (but also of an Matrix) can be selected by passing an array of booleans as column or row mask (and, only for Matrices, also a matrix of booleans)

mask = [false, false, false, false, true, true, true, true]
data[mask,:]
mask = fill(false,nR,nC)
mask[2:end,2:3] .= true
mask
## data[mask] # error !
Matrix(data)[mask]

# Boolean selection can be used to filter on conditions, e.g.:
data[data.Year .>= 2020,:]
data[[i in ["France", "Italy"] for i in data.Country] .&& (data.Year .== 2000),:] # note the parhenthesis

# ### Filtering using the @subset macro from the `DataFramesMacro` package
using DataFramesMeta
@subset(data, :Year .> 2010 )
colToFilter = :Country
@subset(data, :Year .> 2010, cols(colToFilter) .== "France" ) # Conditions are "end" by default. If the column name is embedded in a varaible we eed to use `cols(varname)`

# ### Filtering using the `Query` package
using Query
dfOut = @from i in data begin # `i` is a single row
    @where i.Country == "France" .&& i.Year >= 2000
    ## Select a group of columns, eventually changing their name:
    @select {i.Year, FranceForArea=i.forarea} # or just `i` for the whole row
    @collect DataFrame
end
# long but flexible

# ### Managing missing values

# !!! tip
#     See also the section [`Missingness implementations`](@ref) for a general discussion on missing values

df = copy(data)
## df[3,"forarea"]  = missing # Error, type is Flat64, not Union{Float64,Missing}
df.forarea = allowmissing(df.forarea) # also disallowmissing
allowmissing!(df)
df[3,"forarea"]  = missing
df[6,"forarea"]  = missing
df[6,"Country"]  = missing
nMissings        = length(findall(x -> ismissing(x), df.forarea)) # Count `missing` values in a column.
dropmissing(df)
dropmissing(df[:,["forarea","forvol"]])
collect(skipmissing(df.forarea))
completecases(df)
completecases(df[!,["forarea","forvol"]])
[df[ismissing.(df[!,col]), col] .= 0 for col in names(df) if nonmissingtype(eltype(df[!,col])) <: Number] # Replace `missing` with `0` values in all numeric columns, like `Float64` and `Int64`;
[df[ismissing.(df[!,col]), col] .= "" for col in names(df) if nonmissingtype(eltype(df[!,col])) <: AbstractString] # Replace `missing` with `""` values in all string columns;
df

# ## Editing data
df = copy(data)
df[1,"forarea"] = 11.3
df[[2,4],"forarea"] .= 10
df
push!(df,["UK",2020,5.0,800.0]) # add row
sort!(df,["Country","Year"], rev=false)
df2 = similar(df) # rubish inside
df = similar(df,0) # empty a dataframe. The second parameter is the number of rows desired

# ## Work on dataframe structure

df       = copy(data)
df.foo   = [1,2,3,4,5,6,7,8] # existing or new column
df.volHa = data.forvol ./ data.forarea
df.goo   = Array{Union{Missing,Float64},1}(missing,size(data,1))
select!(df,Not(["foo","volHa","goo"])) # remove  cols by name
rename!(df, ["Country2", "Year2", "forarea2", "forvol2"])
rename!(df, Dict("Country2" => "Country"))
df       = df[:,["Year2", "Country","forarea2","forvol2"] ] #  change column order
insertcols!(df, 2, :foo => [1,2,3,4,5,6,7,8] ) # insert a column at position 2

df.namedYear = map(string,df.Year2)
stringToInt(str) = try parse(Int64, str) catch; return(missing) end; df.Year3 = map(stringToInt, df.namedYear)

df2 = hcat(df,data,makeunique=true)
df3 = copy(data)
df4 = vcat(data,df3)

# ### Categorical data
using CategoricalArrays # You may want to consider also PooledArrays
df.Year2 = categorical(df.Year2)
transform!(df, names(df, AbstractString) .=> categorical, renamecols=false) # transform to categorical all string columns

# !!! warning
#     Attention that while the memory to store the data decreases, and grouping is way more efficient, filtering with categorical values is not necessarily quicker (indeed it can be a bit slower)

levels(df.Year2)
levels!(df.Country,["Sweden","Germany","France","Italy"]) # Let you define a personalised order, useful for ordered data
sort(df.Country)
sort!(df,"Country")
df.Years2 = unwrap.(df.Year2) # convert a categorical array into a normal one.

# ### Joining dataframes
df1,df2 = copy(data),copy(data)
push!(df1,["US",2020,5.0,1000.0])
push!(df2,["China",2020,50.0,1000.0])
rename!(df2,"Year"=>"year")
innerjoin(df1,df2,on=["Country","Year"=>"year"],makeunique=true) # common records only
# Also available: `leftjoin` (all records on left df), `rightjoin` (all on right df), `outerjoin`(all records returned), `semijoin` (like inner join by only with columns from the right df), `antijoin` (left not on right df) and `crossjoin` (like the cartesian product, each on the right by each on the left)

# ## Pivoting data

# _long_ and _wide_ are two kind of layout for equivalent representation of tabled data.
# In a _long_ layout we have each row being an observation of a single variable, and each record is represented as dim1, dim2, ..., value. As the name implies, "long" layouts tend to be relativelly long and hard to analyse by an human, but are very easity to handle.
# At the opposite, A _wide_ layout represents multiple observations on the same row, eventually using multiple horizontal axis as in the next figure (but Julia dataframes handle only a single horizzontal axis):

# |          |              |       |                  |      |
# | -------- | ------------ | ----- | ---------------- | ---- |
# |          | Forest Area  |       |   Forest Volumes |      |
# |          | 2000         | 2020  | 2000             | 2020 |
# | Germany  | 11.35        | 11.42 | 3381             | 3663 |
# | France   | 15.29        | 17.25 | 2254             | 3056 |

# _wide layout is easier to visually analise for a human mind but much more hard to analyse in a sistermatic way. We will learn now how to move from one type of layout to the other.

# ### Stacking columns: from _wide_ to _long_
longDf  = stack(data,["forarea","forvol"])    # we specify the variable to stack (the "measured" variables)
longDf2 = stack(data,Not(["Country","Year"])) # we specify the variables _not_ to stack (the id variables)
longDf3 = stack(data)                         # automatically stack all numerical variables
longDf == longDf2 == longDf3
# Note how the columns `variable` and `value` have been added automatically to host the stachked data

# ### Unstacking columns: from _wide_ to _long_
wideDf = unstack(longDf,["Country","Year"],"variable","value") # args: df, [cols to remains cols also in the wide layout], column with the ids to expand horizontally and column with the relative values
wideDf2 = unstack(longDf,"variable","value") # cols to remains cols also in the wide layout omitted: all cols not to expand and relative value col remains as col
wideDf == wideDf2 == data

# While the DataFrames package doesn't support multiple axis we can still arrive to the table below with a little bit of work by unstacking different columns in separate wide dataframes and then joining or horizontally concatenating them:
wideArea = unstack(data,"Country","Year","forarea")
wideVols = unstack(data,"Country","Year","forvol")
rename!(wideArea,["Country","area_2000","area_2020"])
rename!(wideVols,["Country","vol_2000","vol_2020"])
wideWideDf = outerjoin(wideArea,wideVols,on="Country")

# ## The Split-Apply-Combine strategy

# Aka "divide and conquer". Rather than try to modify the dataset direclty, we first split it in subparts, we work on each subpart and then we recombine them in a target dataset

using Statistics # for `mean`
groupby(data,["Country","Year"]) # The "split" part

# Aggregation:
combine(groupby(data,["Year"]) ,  "forarea" => sum => "sum_area", "forvol" => sum => "sum_vol", nrow)
# ...or...
combine(groupby(data,["Year"])) do subdf # slower
    (sumarea = sum(subdf.forarea), sumvol = sum(subdf.forvol), nCountries = size(subdf,1))
end
# Cumulative computation:
a = combine(groupby(data,["Year"])) do subdf # slower
    (country = subdf.Country, area = subdf.forarea, cumArea = cumsum(subdf.forarea))
end

# Note in these examples that while in the aggregation we was returning a _single record_ for each subgroup (hence we did some dimensionality reduction) in the cumulative compuation we still output the whole subgroup, so the combined dataframe in output has the same number of rows as the original dataframe.

# An alternative approach is to use the `@linq` macro from the `DatAFrameMEta` package that provide a R's `dplyr`-like query language using piped data:
using DataFramesMeta
dfCum = @linq data |>
            groupby([:Year]) |>
            transform(:cumArea = cumsum(:forarea))

# ## Export and saving

# ### DataFrame to Matrix

M = Matrix(data)
# !!! warning
#     Attention that if the dataframe contains different types across columns, the inner type of the matrix will be `Any`

M = Matrix{Union{Float64,Int64,String}}(data)

# ### DataFrame to Dictionary

function toDict(df, dimCols, valueCol)
    toReturn = Dict()
    for r in eachrow(df)
        keyValues = []
        [push!(keyValues,r[d]) for d in dimCols]
        toReturn[(keyValues...,)] = r[valueCol]
    end
    return toReturn
end
dict = toDict(data,["Country","Year"],["forarea","forvol"])
dict["Germany",2000][1]
dict["Germany",2000]["forvol"]
toDict(data,["Country","Year"],"forarea")

# ### DataFrame to NamedTuple

nT = NamedTuple(Dict([Symbol(c) => data[:,c]  for c in names(data)]))        # order not necessarily preserved
using DataStructures
nT = NamedTuple(OrderedDict([Symbol(c) => data[:,c]  for c in names(data)])) # order preserved

# ### Saving as CSV file

CSV.write("outdata.csv",data) # see options at the beginning of segment in the import section and `? CSV.write` for specific export options

# ### Saving as OpenDocument spreadsheet

ods_write("outdata.ods",Dict(("myData",3,2) => data)) # exported starting on cell B3 of sheet "myData"

# ### Saving as Excel spreadsheet
XLSX.writetable("outdata.xlsx",myData = (collect(eachcol(data)),names(data)))

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
using GitCommand

cpp_itensor_version = get(ENV, "CPP_ITENSOR_VERSION", nothing)
if isnothing(cpp_itensor_version)
  error("""Environment variable CPP_ITENSOR_VERSION is not set so C++ ITensor will not be built automatically. Set it from within Julia with `ENV["CPP_ITENSOR_VERSION"] = "3.1.6"`, for example. Or start Julia with `CPP_ITENSOR_VERSION="3.1.6" julia` or `export CPP_ITENSOR_VERSION="3.1.6"; julia`. Once you have set the variable CPP_ITENSOR_VERSION, you can try rebuilding with `] build ITensorsBenchmarks` at the Julia command prompt.""")
end

const deps_dir = @__DIR__
const cpp_itensor_dir = joinpath(deps_dir, "itensor_v$cpp_itensor_version")
if !isfile(joinpath(cpp_itensor_dir, "itensor", "core.h"))
  println("Cloning C++ ITensor and checking out the desired version v$cpp_itensor_version.")
  git() do git
    run(`$git clone https://github.com/ITensor/ITensor.git $cpp_itensor_dir`)
  end
  git() do git
    run(`$git -C $cpp_itensor_dir checkout v$cpp_itensor_version`)
  end
  println("Finished cloning C++ ITensor, next we will try to build it.")
end
if !isfile(joinpath(cpp_itensor_dir, "lib", "libitensor.a"))
  options = joinpath(cpp_itensor_dir, "options.mk")
  if !isfile(options)
    options_sample = joinpath(deps_dir, "options.mk.sample")
    println("Need options.mk to build C++ ITensor.")
    println("First we are copying the file $options_sample to $options.")
    cp(options_sample, options)
  end
  println("Edit $options (for example with `edit(\"$options\")` at the Julia command line) to point to the proper location of your BLAS installation. You can use the instructions at http://www.itensor.org/docs.cgi?page=install&vers=cppv3.")
  println("Building C++ ITensor, this may take some time. Running command `make --directory=$cpp_itensor_dir`")
  cd(cpp_itensor_dir)
  run(`make --directory=$cpp_itensor_dir`)
end


text:
# Dirichlet distribution

import StatsFuns
export Dirichlet

@parameterized Dirichlet(α)

TV.as(d::Dirichlet{(:α,)}) = TV.UnitSimplex(length(d.α))

function basemeasure(μ::Dirichlet{(:α,)})
    t = as(μ)
    d = dimension(t)
    return Pushforward(t, Lebesgue(ℝ)^d, false)
end

@kwstruct Dirichlet(α)

function logdensity(d::Dirichlet{(:α,)}, x)
    α = d.α
    s = 0.0
    for j in eachindex(x)
        s += α[j] * log(x[j])
    end
    return s
end

Base.rand(rng::AbstractRNG, T::Type, μ::Dirichlet) = rand(rng, Dists.Dirichlet(μ.α))

distproxy(d::Dirichlet{(:α,)}) = Dists.Dirichlet(d.α)


text:
module {{{PKG}}}

export question1,
       question2,
       question3,
       question4,
       question5

# using BISC195

include("assignment.jl")
# Write your package code here.

end


text:
@testset "MochaClassification" begin
    using Mocha
    using SGMCMC
    using MLUtilities
    using MochaClassificationWrapper
    using MochaClassification
    using LogisticRegression
    using MochaModelFactories
    # generate artificial log reg model.
    dm = LogisticRegression.artificial_logreg_dm(200,25)
    xtrain = dm.x[1:100,:]'
    ytrain = (dm.y[1:100]+1)/2
    xtest = dm.x[101:end,:]'
    ytest = (dm.y[101:end]+1)/2

    backend = initMochaBackend(false)
    nunits = 1
    model,name = make_dense_nn([nunits],2)
    dm = MochaClassificationModel(xtrain,ytrain,model,backend)
    dmtest = MochaClassificationModel(xtest,ytest,model,backend,do_accuracy=true)# f,or test set accuracy
    dmtraintest = MochaClassificationModel(xtrain,ytrain,model,backend,do_accuracy=true)# for training set accuracy

    # initialised from models
    @test 25nunits+nunits+nunits*2+2 == MochaClassification.fetchnparams(dm)

    # test setparams! and getparams
    d = MochaClassification.fetchnparams(dm)
    x = randn(d)
    MochaClassification.setparams!(dm,x)
    @test MochaClassification.getparams(dm) == x
    g = DataModel.getgrad(dm)
    l = DataModel.getllik(dm)
    for i in 1:10
        x = 0.01*randn(d)
        @test abs(checkgrad(x,l,g)) < 1e-4
    end

    dmtest1 = MochaClassificationModel(xtest,ytest,model,backend,do_accuracy=true)
    dmtest2 = MochaClassificationModel(xtest,ytest,model,backend,do_accuracy=false, do_predprob = true)
    # test evaluate functions
    predprobs = MochaClassification.getpredprobs(dmtest2,x)
    manual_results = MochaClassification.evaluate(dmtest2,predprobs[:predictiveprobs])
    mocha_results = MochaClassification.evaluate(dmtest,x)
    @test manual_results[:accuracy] == mocha_results[:accuracy]
    @test abs(manual_results[:loglikelihood] -mocha_results[:loglikelihood]) < 1
end


text:
module TestSharedResources

using ContextManagers.++
using Test

function test_channel()
    @sync begin
        ch = Channel{Int}(0)
        @with(handle = SharedResource(ch)) do
            for i in 1:3
                context = open(handle)
                @async try
                    @with(ch = context) do
                        put!(ch, i)
                    end
                catch
                    close(ch)
                    rethrow()
                end
            end
        end
        @test sort!(collect(ch)) == 1:3
    end
end

end  # module


text:
int_rules_9_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 9.3*Piecewise*linear*functions =#
    @apply_utils Antiderivative((~u) ^ ~(m'), ~x) => With([c = Simplify(D(~u, ~x))], (1 / c) * Subst(Antiderivative((~x) ^ ~m, ~x), ~x, ~u)) <-- FreeQ(~m, ~x) && PiecewiseLinearQ(~u, ~x)
    @apply_utils Antiderivative(~v / ~u, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (b * ~x) / a - ((b * ~u - a * ~v) / a) * Antiderivative(1 / ~u, ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x)
    @apply_utils Antiderivative((~v) ^ ~n / ~u, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (~v) ^ ~n / (a * ~n) - ((b * ~u - a * ~v) / a) * Antiderivative((~v) ^ (~n - 1) / ~u, ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && (GtQ(~n, 0) && NeQ(~n, 1))
    @apply_utils Antiderivative(1 / (~u * ~v), ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (b / (b * ~u - a * ~v)) * Antiderivative(1 / ~v, ~x) - (a / (b * ~u - a * ~v)) * Antiderivative(1 / ~u, ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x)
    @apply_utils Antiderivative(1 / (~u * sqrt(~v)), ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (2 * atan(sqrt(~v) / Rt((b * ~u - a * ~v) / a, 2))) / (a * Rt((b * ~u - a * ~v) / a, 2)) <-- NeQ(b * ~u - a * ~v, 0) && PosQ((b * ~u - a * ~v) / a)) <-- PiecewiseLinearQ(~u, ~v, ~x)
    @apply_utils Antiderivative(1 / (~u * sqrt(~v)), ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (-2 * atanh(sqrt(~v) / Rt(-((b * ~u - a * ~v)) / a, 2))) / (a * Rt(-((b * ~u - a * ~v)) / a, 2)) <-- NeQ(b * ~u - a * ~v, 0) && NegQ((b * ~u - a * ~v) / a)) <-- PiecewiseLinearQ(~u, ~v, ~x)
    @apply_utils Antiderivative((~v) ^ ~n / ~u, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (~v) ^ (~n + 1) / ((~n + 1) * (b * ~u - a * ~v)) - ((a * (~n + 1)) / ((~n + 1) * (b * ~u - a * ~v))) * Antiderivative((~v) ^ (~n + 1) / ~u, ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && LtQ(~n, -1)
    @apply_utils Antiderivative((~v) ^ ~n / ~u, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], ((~v) ^ (~n + 1) / ((~n + 1) * (b * ~u - a * ~v))) * HypergeometricFunctions._₂F₁(1, ~n + 1, ~n + 2, (-a * ~v) / (b * ~u - a * ~v)) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative(1 / (sqrt(~u) * sqrt(~v)), ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (2 / Rt(a * b, 2)) * atanh((Rt(a * b, 2) * sqrt(~u)) / (a * sqrt(~v))) <-- NeQ(b * ~u - a * ~v, 0) && PosQ(a * b)) <-- PiecewiseLinearQ(~u, ~v, ~x)
    @apply_utils Antiderivative(1 / (sqrt(~u) * sqrt(~v)), ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (2 / Rt(-a * b, 2)) * atan((Rt(-a * b, 2) * sqrt(~u)) / (a * sqrt(~v))) <-- NeQ(b * ~u - a * ~v, 0) && NegQ(a * b)) <-- PiecewiseLinearQ(~u, ~v, ~x)
    @apply_utils Antiderivative((~u) ^ ~m * (~v) ^ ~n, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (-((~u) ^ (~m + 1)) * (~v) ^ (~n + 1)) / ((~m + 1) * (b * ~u - a * ~v)) <-- NeQ(b * ~u - a * ~v, 0)) <-- FreeQ([~m, ~n], ~x) && (PiecewiseLinearQ(~u, ~v, ~x) && (EqQ(~m + ~n + 2, 0) && NeQ(~m, -1)))
    @apply_utils Antiderivative((~u) ^ ~m * (~v) ^ ~(n'), ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], ((~u) ^ (~m + 1) * (~v) ^ ~n) / (a * (~m + 1)) - ((b * ~n) / (a * (~m + 1))) * Antiderivative((~u) ^ (~m + 1) * (~v) ^ (~n - 1), ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- FreeQ([~m, ~n], ~x) && (PiecewiseLinearQ(~u, ~v, ~x) && (NeQ(~m, -1) && (LtQ(~m, -1) && (GtQ(~n, 0) && Not(ILtQ(~m + ~n, -2) && (FractionQ(~m) || GeQ(2 * ~n + ~m + 1, 0)))) || (IGtQ(~n, 0) && (IGtQ(~m, 0) && LeQ(~n, ~m)) || (IGtQ(~n, 0) && Not(IntegerQ(~m)) || ILtQ(~m, 0) && Not(IntegerQ(~n)))))))
    @apply_utils Antiderivative((~u) ^ ~m * (~v) ^ ~(n'), ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], ((~u) ^ (~m + 1) * (~v) ^ ~n) / (a * (~m + ~n + 1)) - ((~n * (b * ~u - a * ~v)) / (a * (~m + ~n + 1))) * Antiderivative((~u) ^ ~m * (~v) ^ (~n - 1), ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && (NeQ(~m + ~n + 2, 0) && (GtQ(~n, 0) && (NeQ(~m + ~n + 1, 0) && (Not(IGtQ(~m, 0) && (Not(IntegerQ(~n)) || LtQ(0, ~m, ~n))) && Not(ILtQ(~m + ~n, -2))))))
    @apply_utils Antiderivative((~u) ^ ~m * (~v) ^ ~n, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], ((~u) ^ (~m + 1) * (~v) ^ ~n) / (a * (~m + ~n + 1)) - ((~n * (b * ~u - a * ~v)) / (a * (~m + ~n + 1))) * Antiderivative((~u) ^ ~m * (~v) ^ Simplify(~n - 1), ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && (NeQ(~m + ~n + 1, 0) && (Not(RationalQ(~n)) && SumSimplerQ(~n, -1)))
    @apply_utils Antiderivative((~u) ^ ~m * (~v) ^ ~n, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (-((~u) ^ (~m + 1)) * (~v) ^ (~n + 1)) / ((~m + 1) * (b * ~u - a * ~v)) + ((b * (~m + ~n + 2)) / ((~m + 1) * (b * ~u - a * ~v))) * Antiderivative((~u) ^ (~m + 1) * (~v) ^ ~n, ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && (NeQ(~m + ~n + 2, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~u) ^ ~m * (~v) ^ ~n, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (-((~u) ^ (~m + 1)) * (~v) ^ (~n + 1)) / ((~m + 1) * (b * ~u - a * ~v)) + ((b * (~m + ~n + 2)) / ((~m + 1) * (b * ~u - a * ~v))) * Antiderivative((~u) ^ Simplify(~m + 1) * (~v) ^ ~n, ~x) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && (Not(RationalQ(~m)) && SumSimplerQ(~m, 1))
    @apply_utils Antiderivative((~u) ^ ~m * (~v) ^ ~n, ~x) => With([a = Simplify(D(~u, ~x)), b = Simplify(D(~v, ~x))], (((~u) ^ ~m * (~v) ^ (~n + 1)) / (b * (~n + 1) * ((b * ~u) / (b * ~u - a * ~v)) ^ ~m)) * HypergeometricFunctions._₂F₁(-(~m), ~n + 1, ~n + 2, (-a * ~v) / (b * ~u - a * ~v)) <-- NeQ(b * ~u - a * ~v, 0)) <-- PiecewiseLinearQ(~u, ~v, ~x) && (Not(IntegerQ(~m)) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~u) ^ ~(n') * log(~(a') + ~(b') * ~x), ~x) => With([c = Simplify(D(~u, ~x))], (((~u) ^ ~n * (~a + ~b * ~x) * log(~a + ~b * ~x)) / ~b - Antiderivative((~u) ^ ~n, ~x)) - ((c * ~n) / ~b) * Antiderivative((~u) ^ (~n - 1) * (~a + ~b * ~x) * log(~a + ~b * ~x), ~x)) <-- FreeQ([~a, ~b], ~x) && (PiecewiseLinearQ(~u, ~x) && (Not(LinearQ(~u, ~x)) && GtQ(~n, 0)))
    #= Antiderivative((~u)^(~n')*log((~a')+(~b')*(~x))/((~a')+(~b')*(~x)),~x) :=  With([c=Simplify(D(u,x))], u^n*log(a+b*x)^2/(2*b) - c*n/(2*b)*Antiderivative(u^(n-1)*log(a+b*x)^2,x)) <-- FreeQ([a,b],x) &&  PiecewiseLinearQ(u,x) && GtQ(n,0) =#
    @apply_utils Antiderivative((~u) ^ ~(n') * (~(a') + ~(b') * ~x) ^ ~(m') * log(~(a') + ~(b') * ~x), ~x) => With([c = Simplify(D(~u, ~x))], (((~u) ^ ~n * (~a + ~b * ~x) ^ (~m + 1) * log(~a + ~b * ~x)) / (~b * (~m + 1)) - (1 / (~m + 1)) * Antiderivative((~u) ^ ~n * (~a + ~b * ~x) ^ ~m, ~x)) - ((c * ~n) / (~b * (~m + 1))) * Antiderivative((~u) ^ (~n - 1) * (~a + ~b * ~x) ^ (~m + 1) * log(~a + ~b * ~x), ~x)) <-- FreeQ([~a, ~b, ~m], ~x) && (PiecewiseLinearQ(~u, ~x) && (Not(LinearQ(~u, ~x)) && (GtQ(~n, 0) && NeQ(~m, -1))))
end

max_stars_repo_path: src/rules/9 Miscellaneous/9.3 Piecewise linear functions.jl


text:
using FFTW

# projections and energies
@testset "projections" begin
    a = ones(2,2)
    a[1,1] = -1
    a[2,2] = 2
    sqrtI = abs.(fft(a))
    b = zeros(2,2)
    b[1,1] = 1

    # projections
    @test FPRS.pM(b,sqrtI) == [2.5 0.5; 0.5 -0.5]
    @test FPRS.pM0(b,sqrtI) == [2.0 0.0; 0.0 -1.0]
    @test FPRS.pP(a) == [0.0 1.0; 1.0 2.0]
    @test FPRS.pTa(a, 1.6) == [0.0 0.0; 0.0 2.0]
    @test FPRS.pTn(a, 1) == [0.0 0.0; 0.0 2.0]
    @test FPRS.pTn(a, 3) == [0.0 1.0; 1.0 2.0]

    # energies
    @test FPRS.eM(b,sqrtI) == 1.0
    @test FPRS.eM0(b,sqrtI) == 1.0
    @test FPRS.eP(a) == 0.5
    @test FPRS.eTa(a, 1.6) == 1.5
    @test FPRS.eTn(a, 1) == 1.5
    @test FPRS.eTn(a, 3) == 0.5
end

max_stars_repo_name: tsipenyuk/FPRS

max_issues_repo_name: tsipenyuk/FPRS

max_forks_repo_name: tsipenyuk/FPRS


text:
#!/usr/bin/env julia

const DATA_DIR = "../10_data"
const INPUT_FILE    = "input.txt"
# const INPUT_FILE    = "input_test1.txt"  # 24
# const INPUT_FILE    = "input_test2.txt"

macro assert(boolean)
    message = string("Assertion: ", boolean, " failed")
    :($(esc(boolean)) || error($message))
end

function get_groups(jds)
    groups = Array{Tuple{Int, Int}, 1}()
    current_val = nothing
    current_count = 0
    for k in 1:length(jds)
        if jds[k] == current_val
            current_count += 1
        else
            if current_count > 0
                push!(groups, (current_val, current_count))
            end
            current_val = jds[k]
            current_count = 1
        end
    end
    # Add the last group
    push!(groups, (current_val, current_count))
    return groups
end

function part1(diffs)
    println("========== PART 1 ==========")
    diff1_count = length(filter(x -> x == 1, diffs))
    diff3_count = length(filter(x -> x == 3, diffs))
    result = diff1_count * diff3_count
    println("part1 result=$(result)")
end

function part2(diffs)
    println("========== PART 2 ==========")
    groups = get_groups(diffs)
    group1_sizes = [count for (num, count) in groups if num == 1]
    comb_map = Dict(1=>1, 2=>2, 3=>4, 4=>7)  # 3:111,12,21,3. 4:1111,22,13,31,121,211,112
    result = prod(map(n->comb_map[n], group1_sizes))
    println("part2 result=$(result)")
end

function main(filename::String)
    lines = readlines(filename)
    jolts = sort(map(s->parse(Int64, s), lines))
    # diffs = [jolts[k+1] - jolts[k] for k in 1:length(jolts) - 1]
    diffs = [jolts; jolts[end] + 3] - [0; jolts]  # All 1s & 3s (up to 4)

    part1(diffs)
    println()
    part2(diffs)
end

main("$DATA_DIR/$INPUT_FILE")

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
function quadrotor_dynamics!(ẋ,X,u)
      #TODO change concatentations to make faster!
      # Quaternion representation
      # Modified from D. Mellinger, N. Michael, and V. Kumar,
      # "Trajectory generation and control for precise aggressive maneuvers with quadrotors",
      # In Proceedings of the 12th International Symposium on Experimental Robotics (ISER 2010), 2010.

      ## States: X ∈ R^13; q = [s;v]
      # x
      # y
      # z
      # q0
      # q1
      # q2
      # q3
      # xdot
      # ydot
      # zdot
      # omega1
      # omega2
      # omega3

      x = X[1:3]
      q = X[4:7]./norm(X[4:7]) #normalize quaternion
      v = X[8:10]
      omega = X[11:13]

      # Parameters
      m = .5 # mass
      J = Matrix(Diagonal([0.0023; 0.0023; 0.004])) # inertia matrix
      Jinv = Matrix(Diagonal(1.0./[0.0023; 0.0023; 0.004])) # inverted inertia matrix
      g = 9.81 # gravity
      L = 0.1750 # distance between motors

      w1 = u[1]
      w2 = u[2]
      w3 = u[3]
      w4 = u[4]

      kf = 1; # 6.11*10^-8;
      F1 = kf*w1;
      F2 = kf*w2;
      F3 = kf*w3;
      F4 = kf*w4;
      F = [0;0;F1+F2+F3+F4] #total rotor force in body frame

      km = 0.0245;
      M1 = km*w1;
      M2 = km*w2;
      M3 = km*w3;
      M4 = km*w4;
      tau = [L*(F2-F4);L*(F3-F1);(M1-M2+M3-M4)] #total rotor torque in body frame

      ẋ[1:3] = v # velocity in world frame
      ẋ[4:7] = 0.5*qmult(q,[0;omega]) #quaternion derivative
      ẋ[8:10] = [0;0;-g] + (1/m)*qrot(q,F) #acceleration in world frame
      ẋ[11:13] = Jinv*(tau - cross(omega,J*omega)) #Euler's equation: I*ω + ω x I*ω = constraint_decrease_ratio
end

function quadrotor_dynamics_mrp!(ẋ, X, u, params)
      #TODO change concatentations to make faster!
      # Modified Rodrigues Parameter (mrp) representation
      # Modified from D. Mellinger, N. Michael, and V. Kumar,
      # "Trajectory generation and control for precise aggressive maneuvers with quadrotors",
      # In Proceedings of the 12th International Symposium on Experimental Robotics (ISER 2010), 2010.
      # States: X ∈ R^12; r = modified Rodrigues parameters
      # x
      # y
      # z
      # r1
      # r2
      # r3
      # xdot
      # ydot
      # zdot
      # omega1
      # omega2
      # omega3

      x = X[1:3]
      r = X[4:6] # mrp
      v = X[7:9]
      omega = X[10:12]

      # Rotation
      R = mrptodcm(r)

      # Parameters
      m = .5 # mass
      J = Matrix(Diagonal([0.0023; 0.0023; 0.004])) # inertia matrix
      Jinv = Matrix(Diagonal(1.0./[0.0023; 0.0023; 0.004])) # inverted inertia matrix
      g = 9.81 # gravity
      L = 0.1750 # distance between motors

      w1 = u[1]
      w2 = u[2]
      w3 = u[3]
      w4 = u[4]

      kf = 1; # 6.11*10^-8;
      F1 = kf*w1;
      F2 = kf*w2;
      F3 = kf*w3;
      F4 = kf*w4;
      F = [0;0;F1+F2+F3+F4] #total rotor force in body frame

      km = 0.0245;
      M1 = km*w1;
      M2 = km*w2;
      M3 = km*w3;
      M4 = km*w4;
      tau = [L*(F2-F4);L*(F3-F1);(M1-M2+M3-M4)] #total rotor torque in body frame

      ẋ[1:3] = v # velocity in world frame
      ẋ[4:6] = 0.25*((1-r'*r)*ω - 2*cross(ω,r) + 2*(ω'*r)*r) #mrp derivative
      ẋ[7:9] = [0;0;-g]  + (1/m)*(R*F) #acceleration in world frame
      ẋ[11:13] = Jinv * (tau - cross(omega, J * omega)) #Euler's equation: I*ω + ω x I*ω = constraint_decrease_ratio
end

function quadrotor_dynamics(X,u)
      ẋ = zeros(13,1)
      quadrotor_dynamics!(ẋ,X,u)
      ẋ
end

function quadrotor_dynamics_mrp(X,u)
      ẋ = zeros(13,1)
      quadrotor_dynamics_mrp!(ẋ,X,u)
      ẋ
end

## Utilities
"""
@(SIGNATURES)
    Rotate a vector by a quaternion
"""
function qrot(q,r)
      r + 2*cross(q[2:4],cross(q[2:4],r) + q[1]*r)
end

"""
@(SIGNATURES)
    Multiplication of two quaternions (q = [s;v])
"""
function qmult(q1,q2)
      [q1[1]*q2[1] - q1[2:4]'*q2[2:4]; q1[1]*q2[2:4] + q2[1]*q1[2:4] + cross(q1[2:4],q2[2:4])]
end

"""
@(SIGNATURES)
  Transformation from modified Rodrigues parameters (mrp) to rotation matrix.
"""
function mrptodcm(mrp)
    # Converts a vector of Modified Rodrigues Parameters to a Rotation Matrix
    mrp2 = mrp' * mrp;
    S = [0       -mrp[3]   mrp[2];
         mrp[3]   0       -mrp[1];
        -mrp[2]   mrp[1]   0]
    I = diagm(0=>fill(1., 3))
    R = I + (8 * S * S + 4 * (1 - mrp2) * S) / ((1 + mrp2)^2)
    return R
end

# Model
n = 13
m = 4

model = Model(quadrotor_dynamics!,n,m)

# Unconstrained
Q = (1e-1)*Matrix(I,n,n)
Q[4,4] = 1.0; Q[5,5] = 1.0; Q[6,6] = 1.0; Q[7,7] = 1.0
R = (1.0)*Matrix(I,m,m)
# R = (1e-1)*Matrix(I,m,m)
Qf = (1000.0)*Matrix(I,n,n)
tf = 5.0
dt = 0.05

# -initial state
x0 = zeros(n)
x0[1:3] = [0.; 0.; 0.]
q0 = [1.;0.;0.;0.]
x0[4:7] = q0

# -final state
xf = copy(x0)
xf[1:3] = [0.;40.;0.] # xyz position
xf[4:7] = q0

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)

# Model + objective
quadrotor = [model, obj_uncon]

## Constrained

r_quad = 3.0
r_sphere = 3.0
spheres = ((0.,10.,0.,r_sphere),(0.,20.,0.,r_sphere),(0.,30.,0.,r_sphere))
n_spheres = 3

# -control limits
u_min = 0.0
u_max = 10.0

# 3 sphere obstacles
function cI_3obs_quad(c,x,u)
    for i = 1:n_spheres
        c[i] = TrajectoryOptimization.sphere_constraint(x,spheres[i][1],spheres[i][2],spheres[i][3],spheres[i][4]+r_quad)
    end
    c
end

# unit quaternion constraint
function unit_quaternion(c,x,u)
    c[1] = sqrt(x[4]^2 + x[5]^2 + x[6]^2 + x[7]^2) - 1.0
end

obj_uq = TrajectoryOptimization.ConstrainedObjective(obj_uncon,u_min=u_min,u_max=u_max,cE=unit_quaternion)
obj_3obs = TrajectoryOptimization.ConstrainedObjective(obj_uncon,u_min=u_min,u_max=u_max,cI=cI_3obs_quad)#,cE=unit_quaternion)

quadrotor_unit_quaternion = [model, obj_uq]
quadrotor_3obs = [model, obj_3obs, spheres]


text:
module DateLib
using Iterators
import Base.==,Base.-,Base.+,Base.println,Base.show,Base.display; # we must import a method to add methods (as opposed to replacing it)
function MonthMatch(i::Integer)

outString=0;
if i==1
	outString="January"
elseif(i==2)
	outString="February"
elseif(i==3)
	outString="March"
elseif(i==4)
	outString="April"
elseif(i==5)
	outString="May"
elseif(i==6)
	outString="June"
elseif(i==7)
	outString="July"
elseif(i==8)
	outString="August"
elseif(i==9)
	outString="September"
elseif(i==10)
	outString="October"
elseif(i==11)
	outString="November"
elseif(i==12)
	outString="December"
end
return outString;
end
#export DataCheck
function DataCheckInternal(d::Integer,m::Integer,y::Integer)
flag=isLeapYear(y);
resFlag=false;
if (d<1||d>31)
	error("The following Day is Invalid : $d")
elseif(m>12||m<1)
	error("The following Month is Invalid : $m")
elseif(m==2)
	if(flag)
		if (d < 30 && d > 0)
			resFlag = true;
		end
	else
		if (d < 29 && d > 0)
			resFlag = true;
		end
	end
else
	if((m==11)||(m==9)||(m==6)||(m==4))
		if (d < 31 && d > 0)
			resFlag = true;
		end
	else
		if (d < 32 && d > 0)
			resFlag = true;
		end
	end
end
if (!resFlag)
	error("The following Date is Invalid: $d / $m / $y")
end
return resFlag;
end

export DateClass
# package code goes here
type DateClass
   day::Integer
   month::Integer
   year::Integer
   DateClass(day,month,year)= DataCheck(day,month,year) ? new(day,month,year) : error("Invalid Costructor");
end

export println
function println(inDate::DateClass)
println("$(inDate.day) $(MonthMatch(inDate.month)) $(inDate.year)")
end

export display
function display(inDate::DateClass)
display("$(inDate.day) $(MonthMatch(inDate.month)) $(inDate.year)")
end

export show
function show(inDate::DateClass)
show("$(inDate.day) $(MonthMatch(inDate.month)) $(inDate.year)")
end

export dayNumber
function dayNumber(inDate::DateClass)

	y=Integer(inDate.year);
	m=Integer(inDate.month);
	d=Integer(inDate.day);
	m = (m + 9) % 12;
	y = y - div(m,10);
	return (365*y + div(y,4) - div(y,100) + div(y,400) + div((m*306 + 5),10) + ( d - 1 ));

end
export isLeapYear
function isLeapYear(year::Integer)

	return ((year & 3) == 0 && ((year % 25) != 0 || (year & 15) == 0));

end

export is30Month
function is30Month(m::Integer)
	return (((m==11)||(m==9)||(m==6)||(m==4)));
end

export lengthMonth
function lengthMonth(m::Integer,y::Integer)
	return is30Month(m) ? 30 : (m==2 ? isLeapYear(y)+28 : 31);
end

export isEndMonth
function isEndMonth(inDate::DateClass)
		return (lengthMonth(inDate.month,inDate.year)==inDate.day);
end

export dayact
function dayact(Data1::DateClass,Data2::DateClass)
	D1= dayNumber(Data1);
	D2= dayNumber(Data2);
	dayCount=0;
	if (D1>D2)
		dayCount = D1 - D2;
	else
		dayCount=D2-D1;
	end
	return dayCount;
end
#########OPERATORS
==(x::DateClass, y::DateClass) = ((x.day==y.day)&&(x.month==y.month)&&(x.year==y.year))
-(x::DateClass, y::DateClass) = dayact(y,x);
+(x::DateClass, y::Integer) = nth(iterate(addoneday,x),y+1);
-(x::DateClass, y::Integer) = nth(iterate(suboneday,x),y+1);

export isLastOfFebruary
function isLastOfFebruary(inDate::DateClass)
	return (isLeapYear(inDate.year)&&(inDate.day==29)&&(inDate.month==2))||((!isLeapYear(inDate.year))&&(inDate.day==28)&&(inDate.month==2));
end

export yearfrac
function yearfrac(startDate::DateClass,endDate::DateClass,convention::Integer)
yearFrac=0.0;
tmpDate=1;
if (convention==0) 		#(ACT/ACT)
	Nday=dayact(startDate,endDate);
	yearFrac=(Nday)/ dayact(startDate,DateClass(startDate.day,startDate.month,startDate.year+1));
elseif(convention==1)  	#(30/360 SIA)
	y1=startDate.year;
	m1=startDate.month;
	y2=endDate.year;
	m2=endDate.month;
	d1=startDate.day;
	d2=endDate.day;
	if(isLastOfFebruary(startDate)&&isLastOfFebruary(endDate))
		d2=30;
	end
	if(isLastOfFebruary(startDate)||startDate.day==31)
		d1=30;
	end
	if(d1==30&&d2==31)
		d2=30;
	end
	dy=y2-y1;
	dm=m2-m1;
	dd=d2-d1;
	yearFrac=(360.0*dy+30.0*dm+dd)/360.0;
elseif(convention==2)		#(ACT/360)
	Nday=dayact(startDate,endDate);
	yearFrac=(Nday)/360.0;
elseif(convention==3)		# (ACT/365)
	Nday=dayact(startDate,endDate);
	yearFrac=( Nday)/365.0;
elseif(convention==4)		# (30/360 PSA)
	y1=startDate.year;
	m1=startDate.month;
	y2=endDate.year;
	m2=endDate.month;
	d1=startDate.day;
	d2=endDate.day;
	if((startDate.day==31)||isLastOfFebruary(startDate))
		d1=30;
	end
	if((startDate.day==30||isLastOfFebruary(startDate))&&endDate.day==31)
		d2=30;
	end
	dy=y2-y1;
	dm=m2-m1;
	dd=d2-d1;
	yearFrac=(360.0*dy+30.0*dm+dd)/360.0;

	#########
elseif(convention==5)		#(30/360 ISDA)
	y1=startDate.year;
	m1=startDate.month;
	y2=endDate.year;
	m2=endDate.month;
	if(startDate.day<31)
		d1=startDate.day;
	else
		d1=30;
	end
	if ((endDate.day==31)&(d1>29))
		d2=30;
	else
		d2=endDate.day;
	end
	yearFrac=(360.0*((y2-y1))+30.0*((m2-m1))+(d2-d1))/360.0;

elseif(convention==6)		#(30E/360)
	y1=startDate.year;
	m1=startDate.month;
	y2=endDate.year;
	m2=endDate.month;
	d1=startDate.day;
	d2=endDate.day;
	if(d1==31)
		d1=30;
	end
	if(d2==31)
		d2=30;
	end
	dy=y2-y1;
	dm=m2-m1;
	dd=d2-d1;
	yearFrac=(360.0*dy+30.0*dm+dd)/360.0;

	############
end
return yearFrac;

end

currMaxImplemented=6;

export addXMonth
function addXMonth(inDate::DateClass,shiftM::Integer=1)
if(shiftM<0)
	error("It is possible only to add month")
end
if(shiftM==0)
	return inDate;
end
	d1=inDate.day;
	m1=inDate.month;
	y1=inDate.year;
	shiftY=div(shiftM,12);
	shiftM=shiftM%12;
	if(shiftM>12-m1)
		shiftM-=12;
		shiftY=shiftY+1;
	end
	m1+=shiftM;
	y1+=shiftY;
	d1=d1<=lengthMonth(m1,y1)?d1:lengthMonth(m1,y1);
	return DateClass(d1,m1,y1);
end

export addoneday;
function addoneday(inDate::DateClass)

	d=0;#inDate.day;
	m=inDate.month;
	y=inDate.year;
	if(!isEndMonth(inDate))
		d=inDate.day+1;
	elseif(inDate.month!=12)
		d=1;
		m=inDate.month+1;
	else
		m=1;
		d=1;
		y=inDate.year+1;
	end
	return DateClass(d,m,y);
end

export suboneday;
function suboneday(inDate::DateClass)
	d=0;#inDate.day;
	m=inDate.month;
	y=inDate.year;
	if(inDate.day!=1)
		d=inDate.day-1;
	elseif(inDate.month!=1)
		m=inDate.month-1;
		if(is30Month(m))
			d=30;
		elseif(m==2)

			if(isLeap(DateClass(1,m,y)))
				d=29;
			else
				d=28;
			end
		else
			d=31;
		end
	else
		m=12;
		d=31;
		y=inDate.year-1;
	end
	return DateClass(d,m,y);
end

function DataCheck{T1 <: Number,T2 <: Number,T3 <: Number}(d::T1,m::T2,y::T3)
	if !isa(d,Integer)
		error("Day must be an Integer")
	elseif !isa(m,Integer)
		error("Month must be an Integer")
	elseif !isa(y,Integer)
		error("Year must be an Integer")
	end
	resFlag=DataCheckInternal(d,m,y);
	return resFlag;
end

end # module