text:
@testset "Core" begin
    T = Float64
    @testset "DSProblem" begin
        N = 3
        p = DSProblem{T}(N)
        @test p.N == N
        @test typeof(p.config.search) == NullSearch
        @test typeof(p.config.poll) == LTMADS{T}
        @test p.status.optimization_status == DS.Unoptimized
        @test p.status.optimization_status_string == "Unoptimized"
        @test p.sense == DS.Min
        @test p.config.max_simultaneous_evaluations == 1
        @test typeof(p.config.mesh) == DS.Mesh{T}

        @test typeof(p.cache) == DS.PointCache{T}
        @test isempty(p.cache.costs)
        @test isempty(p.cache.order)

        @test typeof(p.constraints) == DS.Constraints{T}
    end

    @testset "Mesh" begin
        N = 3
        p = DSProblem{T}(N)
        m = DS.Mesh{T}(N)
        @test p.config.mesh.G == m.G == [1 0 0; 0 1 0; 0 0 1]
        @test p.config.mesh.D == m.D == [1 0 0 -1 0 0; 0 1 0 0 -1 0; 0 0 1 0 0 -1]
        @test p.config.mesh.Δᵐ == m.Δᵐ == 1.0
        # More rigourous tests for MeshUpdate! are in test_LTMADS etc.
        DS.MeshUpdate!(p, DS.Unsuccessful)
        @test p.config.mesh.Δᵐ == 0.25
    end

    @testset "min_mesh_size" begin
        N = 3
        p = DSProblem{Float64}(N)
        @test DS.min_mesh_size(p) ≈ 1.1102230246251565e-16
    end

    @testset "max_evals" begin
        #This will likely be changing, so don't bother test yet
    end

    @testset "Setters" begin
        test_point = [1, 0.3]
        test_out = 49
        p = DSProblem{T}(2)
        f = DS.rosenbrock

        # SetObjective
        @test !isdefined(p, :objective)
        SetObjective(p, f)
        @test isdefined(p, :objective)
        @test p.objective(test_point) ≈ test_out

        p = DSProblem{T}(2, sense=DS.Max)
        SetObjective(p, f)
        @test p.objective(test_point) ≈ -test_out

        # SetInitialPoint
        @test_throws ErrorException SetInitialPoint(p, [1.0, 1.0, 1.0])
        @test_throws ErrorException SetInitialPoint(p, [1.0])
        SetInitialPoint(p, [5.0, 5.0])
        @test p.user_initial_point == [5.0, 5.0]

        # Iteration limits
        @test p.stoppingconditions[1].limit == 1000
        @test p.status.iteration == 0

        Optimize!(p)

        @test p.status.iteration == 1000
        @test_throws ErrorException SetIterationLimit(p, 900)
        SetIterationLimit(p, 1100)
        @test p.stoppingconditions[1].limit == 1100
        BumpIterationLimit(p, 200)
        @test p.stoppingconditions[1].limit == 1300

        # Variable Ranges
        p = DSProblem{T}(3)
        @test p.config.meshscale == [1.0, 1.0, 1.0]
        @test_throws ErrorException SetVariableRange(p, 4, -5.0, 5.0)
        SetVariableRange(p, 3, -5.0, 5.0)
        @test p.config.meshscale == [1.0, 1.0, 1.0]
        SetVariableRange(p, 3, -10.0, 10.0)
        @test p.config.meshscale == [1.0, 1.0, 2.0]

        @test_throws ErrorException SetVariableRanges(p, [-5.0, -5.0, -5.0, -5.0], [5.0, 5.0, 5.0, 5.0])
        SetVariableRanges(p, [-5.0, -5.0, -5.0], [5.0, 5.0, 5.0])
        @test p.config.meshscale == [1.0, 1.0, 1.0]
        SetVariableRanges(p, [-10.0, -10.0, -10.0], [10.0, 10.0, 10.0])
        @test p.config.meshscale == [2.0, 2.0, 2.0]

        # Sense
        p = DSProblem{T}(3)
        @test p.sense == DS.Min
        p = DSProblem{T}(3, sense=DS.Min)
        @test p.sense == DS.Min
        p = DSProblem{T}(3, sense=DS.Max)
        @test p.sense == DS.Max

        p = DSProblem{T}(3)
        SetSense(p, DS.Max)
        @test p.sense == DS.Max
        SetSense(p, DS.Min)
        @test p.sense == DS.Min
    end

    @testset "Optimize!" begin
    # TODO need to think about how to test this properly.
    # It is run several times during testing, so if it is broken then it should come up,
    # but it should still be tested.
    end

    @testset "EvaluatePoint!" begin
        function setup(con; initial=[4.0, 4.0, 4.0])
            p = DSProblem{T}(3)
            f = DS.rosenbrock
            SetObjective(p, f)
            SetInitialPoint(p, initial)
            con(p)
            DS.EvaluateInitialPoint(p)
            return p
        end
        p = setup(x->Nothing)
        # empty
        @test DS.EvaluatePoint!(p, []) == DS.Unsuccessful
        # Extreme constraints
        AddExtremeConstraint(p, x->x[1] > 0)
        @test DS.EvaluatePoint!(p, [[3.0, 3.0, 3.0]]) == DS.Dominating
        @test DS.EvaluatePoint!(p, [[-1.0, 3.0, 3.0]]) == DS.Unsuccessful

        p = setup(p->AddExtremeConstraint(p, x->x[1] > 0))
        @test DS.EvaluatePoint!(p, [[3.0, 3.0, 3.0],[-1.0, 3.0, 3.0]]) == DS.Dominating

        # Progressive Constraints

        #Feasible point with cost reduction is dominating
        p = setup(p->AddProgressiveConstraint(p, x->x[1] > 0 ? 0 : -x[1]))
        @test DS.EvaluatePoint!(p, [[3.0, 3.0, 3.0]]) == DS.Dominating

        #Feasible point with cost increase is unsuccessful
        p = setup(p->AddProgressiveConstraint(p, x->x[1] > 0 ? 0 : -x[1]))
        @test DS.EvaluatePoint!(p, [[5.0, 4.0, 4.0]]) == DS.Unsuccessful

        #Infeasible point with worse cost but less violation is improving
        p = setup(p->AddProgressiveConstraint(p, x->-x[1]), initial=[-6.0, 4.0, 4.0])
        @test DS.EvaluatePoint!(p, [[-5.0, 11.0, 4.0]]) == DS.Improving

        #Infeasible point with improving cost and less violation is dominating
        p = setup(p->AddProgressiveConstraint(p, x->x[1] > 0 ? 0 : -x[1]))
        @test DS.EvaluatePoint!(p, [[-1.0, 4.0, 4.0]]) == DS.Dominating

        #Infeasible point with improving cost and greater violation is unsuccessful
        p = setup(p->AddProgressiveConstraint(p, x->x[1] > 0 ? 0 : -x[1]))
        p.constraints.collections[2].h_max = -4.0
        @test DS.EvaluatePoint!(p, [[-1.0, 4.0, 4.0]]) == DS.Unsuccessful
    end

    @testset "Opportunistic Evaluation" begin
        #Initialisation
        p = DSProblem(3)
        @test p.config.opportunistic == false

        p = DSProblem(3, opportunistic=false)
        @test p.config.opportunistic == false

        p = DSProblem(3, opportunistic=true)
        @test p.config.opportunistic == true

        #Setter
        p = DSProblem(3)
        @test p.config.opportunistic == false
        SetOpportunisticEvaluation(p)
        @test p.config.opportunistic == true
        SetOpportunisticEvaluation(p, opportunistic=true)
        @test p.config.opportunistic == true
        SetOpportunisticEvaluation(p, opportunistic=false)
        @test p.config.opportunistic == false

        #Efficacy

        #Confirm that all points are checked and best is taken
        p = DSProblem(3; objective=x->sum(x.^2), initial_point=[10, 10, 10])
        SetIterationLimit(p, 1000) #Needed for setup to run
        DS.Setup(p)
        SetOpportunisticEvaluation(p, opportunistic=false)
        @test p.x == [10.0, 10.0, 10.0]
        DS.EvaluatePoint!(p, [[9.0,9.0,9.0], [8.0,8.0,8.0], [5.0,5.0,5.0], [11.0,11.0,11.0]])
        @test p.x == [5.0,5.0,5.0]

        p = DSProblem(3; objective=x->sum(x.^2), initial_point=[10, 10, 10])
        SetIterationLimit(p, 1000) #Needed for setup to run
        DS.Setup(p)
        SetOpportunisticEvaluation(p, opportunistic=true)
        @test p.x == [10.0, 10.0, 10.0]
        DS.EvaluatePoint!(p, [[9.0,9.0,9.0], [8.0,8.0,8.0], [5.0,5.0,5.0], [11.0,11.0,11.0]])
        @test p.x == [9.0,9.0,9.0]
    end
end


text:
function FieldRotation(arg0::Field, arg1::Rotation)
    return FieldRotation((Field, Rotation), arg0, arg1)
end

function FieldRotation(arg0::FieldVector3D, arg1::FieldVector3D)
    return FieldRotation((FieldVector3D, FieldVector3D), arg0, arg1)
end

function FieldRotation(arg0::FieldVector3D, arg1::FieldVector3D, arg2::FieldVector3D, arg3::FieldVector3D)
    return FieldRotation((FieldVector3D, FieldVector3D, FieldVector3D, FieldVector3D), arg0, arg1, arg2, arg3)
end

function FieldRotation(arg0::FieldVector3D, arg1::RealFieldElement, arg2::RotationConvention)
    return FieldRotation((FieldVector3D, RealFieldElement, RotationConvention), arg0, arg1, arg2)
end

function FieldRotation(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::jboolean)
    return FieldRotation((RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, jboolean), arg0, arg1, arg2, arg3, arg4)
end

function FieldRotation(arg0::RotationOrder, arg1::RotationConvention, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement)
    return FieldRotation((RotationOrder, RotationConvention, RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2, arg3, arg4)
end

function FieldRotation(arg0::Vector{Vector{RealFieldElement}}, arg1::jdouble)
    return FieldRotation((Vector{Vector{RealFieldElement}}, jdouble), arg0, arg1)
end

function apply_inverse_to(::Type{FieldRotation}, arg0::Rotation, arg1::FieldRotation)
    return jcall(FieldRotation, "applyInverseTo", FieldRotation, (Rotation, FieldRotation), arg0, arg1)
end

function apply_inverse_to(::Type{FieldRotation}, arg0::Rotation, arg1::FieldVector3D)
    return jcall(FieldRotation, "applyInverseTo", FieldVector3D, (Rotation, FieldVector3D), arg0, arg1)
end

function apply_inverse_to(obj::FieldRotation, arg0::FieldRotation)
    return jcall(obj, "applyInverseTo", FieldRotation, (FieldRotation,), arg0)
end

function apply_inverse_to(obj::FieldRotation, arg0::FieldVector3D)
    return jcall(obj, "applyInverseTo", FieldVector3D, (FieldVector3D,), arg0)
end

function apply_inverse_to(obj::FieldRotation, arg0::Rotation)
    return jcall(obj, "applyInverseTo", FieldRotation, (Rotation,), arg0)
end

function apply_inverse_to(obj::FieldRotation, arg0::Vector3D)
    return jcall(obj, "applyInverseTo", FieldVector3D, (Vector3D,), arg0)
end

function apply_inverse_to(obj::FieldRotation, arg0::Vector{RealFieldElement}, arg1::Vector{RealFieldElement})
    return jcall(obj, "applyInverseTo", void, (Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1)
end

function apply_inverse_to(obj::FieldRotation, arg0::Vector{jdouble}, arg1::Vector{RealFieldElement})
    return jcall(obj, "applyInverseTo", void, (Vector{jdouble}, Vector{RealFieldElement}), arg0, arg1)
end

function apply_to(::Type{FieldRotation}, arg0::Rotation, arg1::FieldRotation)
    return jcall(FieldRotation, "applyTo", FieldRotation, (Rotation, FieldRotation), arg0, arg1)
end

function apply_to(::Type{FieldRotation}, arg0::Rotation, arg1::FieldVector3D)
    return jcall(FieldRotation, "applyTo", FieldVector3D, (Rotation, FieldVector3D), arg0, arg1)
end

function apply_to(obj::FieldRotation, arg0::FieldRotation)
    return jcall(obj, "applyTo", FieldRotation, (FieldRotation,), arg0)
end

function apply_to(obj::FieldRotation, arg0::FieldVector3D)
    return jcall(obj, "applyTo", FieldVector3D, (FieldVector3D,), arg0)
end

function apply_to(obj::FieldRotation, arg0::Rotation)
    return jcall(obj, "applyTo", FieldRotation, (Rotation,), arg0)
end

function apply_to(obj::FieldRotation, arg0::Vector3D)
    return jcall(obj, "applyTo", FieldVector3D, (Vector3D,), arg0)
end

function apply_to(obj::FieldRotation, arg0::Vector{RealFieldElement}, arg1::Vector{RealFieldElement})
    return jcall(obj, "applyTo", void, (Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1)
end

function apply_to(obj::FieldRotation, arg0::Vector{jdouble}, arg1::Vector{RealFieldElement})
    return jcall(obj, "applyTo", void, (Vector{jdouble}, Vector{RealFieldElement}), arg0, arg1)
end

function compose(obj::FieldRotation, arg0::FieldRotation, arg1::RotationConvention)
    return jcall(obj, "compose", FieldRotation, (FieldRotation, RotationConvention), arg0, arg1)
end

function compose(obj::FieldRotation, arg0::Rotation, arg1::RotationConvention)
    return jcall(obj, "compose", FieldRotation, (Rotation, RotationConvention), arg0, arg1)
end

function compose_inverse(obj::FieldRotation, arg0::FieldRotation, arg1::RotationConvention)
    return jcall(obj, "composeInverse", FieldRotation, (FieldRotation, RotationConvention), arg0, arg1)
end

function compose_inverse(obj::FieldRotation, arg0::Rotation, arg1::RotationConvention)
    return jcall(obj, "composeInverse", FieldRotation, (Rotation, RotationConvention), arg0, arg1)
end

function distance(::Type{FieldRotation}, arg0::FieldRotation, arg1::FieldRotation)
    return jcall(FieldRotation, "distance", RealFieldElement, (FieldRotation, FieldRotation), arg0, arg1)
end

function get_angle(obj::FieldRotation)
    return jcall(obj, "getAngle", RealFieldElement, ())
end

function get_angles(obj::FieldRotation, arg0::RotationOrder, arg1::RotationConvention)
    return jcall(obj, "getAngles", Vector{RealFieldElement}, (RotationOrder, RotationConvention), arg0, arg1)
end

function get_axis(obj::FieldRotation, arg0::RotationConvention)
    return jcall(obj, "getAxis", FieldVector3D, (RotationConvention,), arg0)
end

function get_identity(::Type{FieldRotation}, arg0::Field)
    return jcall(FieldRotation, "getIdentity", FieldRotation, (Field,), arg0)
end

function get_matrix(obj::FieldRotation)
    return jcall(obj, "getMatrix", Vector{Vector{RealFieldElement}}, ())
end

function get_q0(obj::FieldRotation)
    return jcall(obj, "getQ0", RealFieldElement, ())
end

function get_q1(obj::FieldRotation)
    return jcall(obj, "getQ1", RealFieldElement, ())
end

function get_q2(obj::FieldRotation)
    return jcall(obj, "getQ2", RealFieldElement, ())
end

function get_q3(obj::FieldRotation)
    return jcall(obj, "getQ3", RealFieldElement, ())
end

function revert(obj::FieldRotation)
    return jcall(obj, "revert", FieldRotation, ())
end

function to_rotation(obj::FieldRotation)
    return jcall(obj, "toRotation", Rotation, ())
end

ext: jl


text:
module Addition

require("Corpora.jl")
require("CommonTools.jl")
require("CVMAbstract.jl")
require("DataTypes.jl")

## Imports
using DataTypes
using Base.LinAlg.BLAS.axpy!
using Corpora.Token
using Corpora.text
using CommonTools: tanh, square!, l2sq, l2sq_distance, mult_by_0!, add!
using CommonTools.CVMIntermediates

using CVMAbstractModel
import CVMAbstractModel: compose!, store_gradients!

## Exports
export AdditionModel
export compose, step!, compose_w_l2
export compose!
export store_gradients!
export AdditionParameters

immutable AdditionModel <: AbstractModel
    wordvecdim::Int64;
    margin::Float;
    regularizer::Float;
    lr_init::Float;
    lr_min::Float;
    niter::Int64;
    gamma::Float;
    nnoise::Int64;
    batchsize::Int64;
end

type AdditionParameters <: AbstractParameters
    wordvecs::Matrix{Float};
    gradients::Matrix{Float};
end

AdditionParameters() = AdditionParameters(Array(Float,0,0), Array(Float,0,0))

## Composes the given sentence (s) and stores the result in o
function compose_w_l2!(model::AdditionModel, s::Vector{Token},
                            o::Vector{Float})
    """
        Compose a vector representation for a sentence while calculating
        the l2 norm of the sentence's word vectors at the same time.
        The sentence is given as a list of tokens.
    """

    # Make sure that o is initialized to 0.0 (assumes that o does
    # not contains NaN or Inf values
    Base.LinAlg.BLAS.scal!(model.wordvecdim,0.0,o,1);
    l2reg = 0.0
    for iword=1:size(s,1)
        v = s[iword].vector
        add!(v, o)
        l2reg += l2sq(v)*0.5
    end
    return o, l2reg
end

function compose_w_l2(model::AdditionModel, s::Vector{Token})
    return compose_w_l2(model, s, zeros(Float, model.wordvecdim));
end

# Composes the given sentence s and store the result in o.
# Composing meas simply adding up all the word vectors
function compose!(model::AdditionModel, s::Vector{Token}, o::Vector{Float})
    """
        Compose a sentence vector from a list of tokens while
        storing intermediates results.
        In-place
    """
    # Make sure that the result is initialized to all zeros
    # Undefined behavior if o contains NaN or Inf
    mult_by_0!(o);
    for iword=1:size(s,1)
        add!(s[iword].vector, o)
    end
    return o
end

function compose!(model::AdditionModel, s::Vector{Token},
                                intermediate::CVMIntermediates)
    return compose!(model, s, intermediate.result)
end

# Compose the given sentence and return the result
function compose(model::AdditionModel, s::Vector{Token})
    """
        Compose a sentence vector from a list of tokens while
        storing intermediates results.
    """
    o = zeros(Float, model.wordvecdim)
    for iword=1:size(s,1)
        o += s[iword].vector
    end
    return o
end

# Propagate the given gradient down to the word representations
function store_gradients!(model::AdditionModel,
                    pre_grad::Vector{Float},
                    s::Vector{Token}, intermediate::CVMIntermediates;
                    intermediates_precomputed::Bool=false)
    """
        Given a gradient from a higher level in the neural network,
        back-propagate the gradient through to the word representations.
    """
    # Since the composition function is addition, we simply
    # need to add the given gradient to the word repr gradients
    for token in s
        add!(pre_grad, token.gradient)
    end
end

end

max_stars_repo_name: ogh/binclusion

max_issues_repo_name: ogh/binclusion

max_forks_repo_name: ogh/binclusion


text:
function autograd{T <: Number}(
    f::Function,
    x::T,
)
    return epsilon(f(Dual(x, one(T))))
end

function autograd!{T <: Real}(
    gr::Vector{T},
    x_d::Vector{Dual{T}},
    f::Function,
    x::Vector{T},
)
    n = length(x)
    for i in 1:n
        x_d[i] = Dual(x[i], zero(T))
    end
    for i in 1:n
        # Why isn't this x[i]
        x_d[i] = Dual(real(x_d[i]), one(T))
        result = f(x_d)
        gr[i] = epsilon(result)
        # Why isn't this x[i]
        x_d[i] = Dual(real(x_d[i]), zero(T))
    end
    return gr
end


text:
using Millboard
using Test

Millboard.set_table_mode(:grid)
Base.eval(:(have_color = false))

types = [AbstractString, String]
board = [T.abstract for T=types]

@test """
+----------------+----------+
|                | abstract |
+================+==========+
| AbstractString |     true |
+----------------+----------+
|         String |    false |
+----------------+----------+""" == table(board, colnames=["abstract"], rownames=types) |> string

Base.eval(:(have_color = true))


text:
# probSNFE has to receive Input1D/Input2D as inputs
# probSNFE throws the ProbOutput1D/ProbOutput2D as outputs (which is of the inputs of solveSNFE)

# Structures for the probSNFE input
struct Input1D{T<:Real,P<:Signed,O,fI,fK,fS,fV0} <: OneDim
          α    :: T
          v    :: T
          V0   :: fV0
          L    :: O
          N    :: P
          T    :: T
          n    :: P
    extInput   :: fI
    kernel     :: fK
    firingRate :: fS
end

struct Input2D{T<:Real,P<:Signed,O,fI,fK,fS,fV0} <: TwoDim
    α          :: T
    v          :: T
    V0         :: fV0
    L          :: O
    N          :: P
    T          :: T
    n          :: P
    extInput   :: fI
    kernel     :: fK
    firingRate :: fS
end

# Structures for the probSNFE output
struct ProbOutput1D{T<:AbstractFloat,P<:Vector{<:Complex},A<:Domain,fI,fS} <: OneDim
    Plan    :: FFTW.rFFTWPlan{Float64,-1,false,1,UnitRange{Int64}}
    PlanInv :: AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.rFFTWPlan{Complex{Float64},1,false,1,UnitRange{Int64}},Float64}
    krings  :: P
    V0      :: Vector{T}
    v0      :: P
    sv      :: P
    Ω       :: A
    α       :: T
    I       :: fI
    S       :: fS
end

struct ProbOutput2D{T<:AbstractFloat,P<:Matrix{<:Complex},A<:Domain,fI,fS} <: TwoDim
    Plan    :: FFTW.rFFTWPlan{Float64,-1,false,2,UnitRange{Int64}}
    PlanInv :: AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.rFFTWPlan{Complex{Float64},1,false,2,UnitRange{Int64}},Float64}
    krings  :: P
    V0      :: Matrix{T}
    v0      :: P
    sv      :: P
    Ω       :: A
    α       :: T
    I       :: fI
    S       :: fS
end


text:
# ----------- Auxiliary functions to generate sparse grid type stuff

"""
`function init1pspec!` : initialize the specification of the 1-particle basis,
generates all possible 1-p basis functions, sorted by degree.
"""
function init1pspec!(B1p::OneParticleBasis,
                     Bsel::DownsetBasisSelector = MaxBasis(1))
   syms = tuple(symbols(B1p)...)
   rgs = indexrange(B1p)
   lens = [ length(rgs[sym]) for sym in syms ]
   spec = []
   maxlev = maxlevel1(Bsel, B1p)
   for I in CartesianIndices(ntuple(i -> 1:lens[i], length(syms)))
      J = ntuple(i -> rgs[syms[i]][I.I[i]], length(syms))
      b = NamedTuple{syms}(J)
      # check whether valid
      if isadmissible(b, B1p)
         if level1(b, Bsel, B1p) <= maxlev
            push!(spec, b)
         end
      end
   end
   sort!(spec, by = b -> level(b, Bsel, B1p))
   return set_spec!(B1p, spec)
end

gensparse(N::Integer, deg::Real; degfun = ν -> sum(ν), kwargs...) =
      gensparse(; NU=N, admissible = ν -> (degfun(ν) <= deg), kwargs...)

"""
`gensparse(...)` : utility function to generate high-dimensional sparse grids
which are downsets.

All arguments are keyword arguments (with defaults):
* `NU` : maximum correlation order
* `minvv = 0` : `minvv[i] gives the minimum value for `vv[i]`
* `maxvv = Inf` : `maxvv[i] gives the minimum value for `vv[i]`
* `tup2b = vv -> vv` :
* `admissible = _ -> false` : determines whether a tuple belongs to the downset
* `filter = _ -> true` : a callable object that returns true of tuple is to be kept and
false otherwise (whether or not it is part of the downset!) This is used, e.g.
to enfore conditions such as ∑ lₐ = even or |∑ mₐ| ≦ M
* `INT = Int` : integer type to be used
* `ordered = false` : whether only ordered tuples are produced; ordered tuples
correspond to  permutation-invariant basis functions
"""
gensparse(; NU::Integer = nothing,
            minvv = [0 for _=1:NU],
            maxvv = [Inf for _=1:NU],
            tup2b = vv -> vv,
            admissible = _-> false,
            filter = _-> true,
            INT = Int,
            ordered = false) =
      _gensparse(Val(NU), tup2b, admissible, filter, INT, ordered,
                 SVector(minvv...), SVector(maxvv...))

"""
`_gensparse` : function barrier for `gensparse`
"""
function _gensparse(::Val{NU}, tup2b, admissible, filter, INT, ordered,
                    minvv, maxvv) where {NU}
   @assert INT <: Integer

   lastidx = 0
   vv = @MVector zeros(INT, NU)
   for i = 1:NU; vv[i] = minvv[i]; end

   spec = SVector{NU, INT}[]
   orig_spec = SVector{NU, INT}[]

   # special trivial case - this should actually never occur :/
   # here, we just push an empty vector provided that the constant term
   # is even allowed.
   if NU == 0
      if all(minvv .== 0) && admissible(vv) && filter(vv)
         push!(spec, SVector(vv))
      end
      return spec
   end

   while true
      # check whether the current vv tuple is admissible
      # the first condition is that its max index is small enough
      # we want to increment `curindex`, but if we've reach the maximum degree
      # then we need to move to the next index down
      isadmissible = true
      if any(vv .> maxvv)
         isadmissible = false
      else
         bb = tup2b(vv)
         isadmissible = admissible(bb)
      end

      if isadmissible
         # ... then we add it to the stack  ...
         # (unless some filtering mechanism prevents it)
         if filter(bb)
            push!(spec, SVector(vv))
            push!(orig_spec, copy(SVector(vv)))
         end
         # ... and increment it
         lastidx = NU
         vv[lastidx] += 1
      else
         if lastidx == 0
            error("""lastidx == 0 should never occur; this means that the
                     smallest basis function is already inadmissible and therefore
                     the basis is empty.""")
         end

         # we have overshot, e.g. level(vv) > maxlevel or something like this
         # we must go back down, by decreasing the index at which we increment
         if lastidx == 1
            # if we have gone all the way down to lastindex==1 and are still
            # inadmissible then this means we are done
            break
         end
         # reset
         vv[lastidx-1] += 1
         if ordered   # ordered tuples (permutation symmetry)
            vv[lastidx:end] .= vv[lastidx-1]
         else         # unordered tuples (no permutation symmetry)
            vv[lastidx:end] .= 0
         end
         lastidx -= 1
      end
   end

   if ordered
      # sanity check, to make sure all is as intended...
      @assert all(issorted, orig_spec)
      @assert length(unique(orig_spec)) == length(orig_spec)
   end

   # here we used to remove the constant term in the past, but this should now
   # be done via the filtering mechanism.

   return spec
end


text:
export dcs_lp, charikar_lp, dcs_lp_extract, extract_primal_solution, charikar_JuMP_LP
export LP_METHODS

using JuMP
using Gurobi
using Convex
LP_METHODS=["DCS"; "MIN"; "AVG"]

function dcs_convert{T}(G::Array{SparseMatrixCSC{T, Int64}, 1})
  m = length(G)
  H = SparseMatrixCSC{Float64, Int64}[]
  for j=1:m
    push!(H, similar(G[j], Float64))
    ii, jj = findn(G[j])
    for k=1:length(ii)
       H[j][ii[k], jj[k] ]  = 1.0
    end
  end
  return H
end

function old_dcs_lp{T}(H::Array{SparseMatrixCSC{T, Int64}, 1})
  G = dcs_convert(H)
  m = length(G)
  n = size(G[1], 1)
  X = Convex.Variable[]
  for j=1:m
    push!(X, Variable(n, n,Positive()))
  end
  Y = Variable(n, Positive())
  t = Variable()
  z = sum(Y)
  constraints = Convex.Constraint[]
  push!(constraints, z <= 1.0)

  for j=1:m
    push!(constraints, trace(G[j] * X[j]) >= t)
    ii, jj = findn(G[j])
    for k=1:length(ii)
      crow = ii[k]
      ccol = jj[k]
      push!(constraints, X[j][crow, ccol] <= Y[crow])
      push!(constraints, X[j][crow, ccol] <= Y[ccol])
    end
  end

  problem  = maximize(t, constraints);
  solve!(problem);
  # println("Optval: ", problem.optval)
  return problem, X, Y, t
end

function dcs_lp{T}(G::Array{SparseMatrixCSC{T, Int64}, 1})

  m = length(G)
  n = size(G[1], 1)
  nz = Int64[]
  @inbounds begin
    iv = {}
    jv = {}
    gv = {}
    for j=1:m
      ii, jj= findn(G[j])
      push!(iv, ii)
      push!(jv, jj)
      push!(nz, length(ii))
      ss = similar(ii, Float64)
      for k=1:length(ii)
        ss[k] = G[j][ii[k], jj[k] ]
      end
      push!(gv, ss)
    end
    ni = deepcopy(nz)
    insert!(ni, 1, 1)
    ni = cumsum(ni)

    nv = sum(nz)
    X = Convex.Variable(nv, 1 , Positive())

    Y = Variable(n, Positive())
    t = Variable()
    w = sum(Y)
    constraints = Convex.Constraint[]
    push!(constraints, w <= 1.0)

    for j=1:m
      cb = ni[j]
      ce = ni[j+1]-1
      # println("gv: $(length(gv[j])) ni: [$cb:$ce]")
      # @assert length(gv[j]) == length([cb:ce;])
      # @assert cb >= 1
      # @assert ce <= length(X)
      push!(constraints, sum(X[cb:ce] .* gv[j] ) >= t)
      # @assert length(gv[j]) == nz[j]
      for k=1:nz[j]
        crow = iv[j][k]
        ccol = jv[j][k]
        ci = ni[j] - 1 + k
        # @assert (ci >= cb) && (ci <= ce)
        push!(constraints, X[ci] <= Y[crow])
        push!(constraints, X[ci] <= Y[ccol])
        # println("j: $j X[$ci] <= Y[$crow], Y[$ccol]")
      end
    end
  end
  problem  = maximize(t, constraints);
  solve!(problem);

  Z = SparseMatrixCSC{Float64, Int64}[]
  typeof(Z)
  for j=1:m
    push!(Z, spzeros(Float64, n, n))
  end
  for j=1:m
    for k=1:nz[j]
      crow = iv[j][k]
      ccol = jv[j][k]
      ci = ni[j] - 1 + k
      # @assert crow >= 1 && crow <= n
      # @assert ccol >= 1 && ccol <= n
      # @assert G[j][crow, ccol] != 0
      # println(X.value[ci])
      Z[j][crow, ccol] = X.value[ci]
    end
  end

  # println("Optval: ", problem.optval)
  return problem, Z, Y.value, t.value
end

function faster_dcs_lp{T}(G::Array{SparseMatrixCSC{T, Int64}, 1})
  m = length(G)
  n = size(G[1], 1)
  nz = Int64[]
  @inbounds begin
    for j=1:m
      cval = nnz(G[j])
      push!(nz, cval)
    end

    ni = deepcopy(nz)
    insert!(ni, 1, 1)
    ni = cumsum(ni)

    X = Convex.Variable(sum(nz), 1 , Positive())
    Y = Variable(n, Positive())
    t = Variable()
    z = sum(Y)
    constraints = Convex.Constraint[]
    push!(constraints, z <= 1.0)

    for j=1:m
      cb = ni[j]
      ce = (ni[j+1]-1)
      gi = find(G[j])
      gval = map((x) -> convert(Float64, x), G[j][gi])
      push!(constraints, sum(X[cb:ce] .* gval) >= t)
      for k=1:nz[j]
        ccol = iceil(gi[k]/n)
        crow = (gi[k] % n)
        if crow == 0
          crow = n
        end
        ci = ni[j] - 1 + k

        push!(constraints, X[ci] <= Y[crow])
        push!(constraints, X[ci] <= Y[ccol])
        # println("j: $j X[$ci] <= Y[$crow], Y[$ccol]")
      end
    end
  end
  problem  = maximize(t, constraints);
  solve!(problem);
  return problem, X, Y, t
end

# function DCSP(G::Graph[])
#   m = length(G)
#   n = size(G[1], 1)
#   R = Model(solver=GurobiSolver(OutputFlag=1))
#   @defVar(R, Y[1:n] >= 0)
#   for i=1:m
#     [gi, gj] = findn(G[i])
#     ne = length(gi)
#     @defVar(R, X[i=i, j=1:ne] >= 0)

# end

function charikar_JuMP_LP(G::Graph)
  n = size(G, 1)
  tic()
  ii, jj = findn(G)
  ne = length(ii)
  println("graph with $n nodes and $ne edges. Setting up Charikar LP ... ")
  R = Model(solver=GurobiSolver(OutputFlag=0))
  @defVar(R, X[1:ne] >= 0)
  @defVar(R, Y[1:n] >= 0)
  @addConstraint(R, gamma, sum{Y[i], i=1:n} <= 1)

  @addConstraint(R, alpha[k=1:ne], X[k] <= Y[ ii[k] ] )
  @addConstraint(R, beta[k=1:ne], X[k] <= Y[ jj[k] ] )
#   for k=1:ne
#     crow = ii[k]
#     ccol = jj[k]
#     @addConstraint(R, X[crow, ccol] <= Y[crow])
#     @addConstraint(R, X[crow, ccol] <= Y[ccol])
#   end
  @setObjective(R, Max, sum{X[k], k=1:ne })
  t_elapsed = toc();
  println("Completed setup, starting to solve LP with Gurobi ")
  status = solve(R);
  if status == :Optimal
    println("solved LP successfully, returning primal variables")
  end

  t = R.objVal
  X1 = spzeros(n, n)
  Y1 = spzeros(n, 1)
  for k=1:ne
    X1[ii[k], jj[k]]= getValue(X[k])
  end
  for j=1:n
    Y1[j] = getValue(Y[j])
  end
  return X1, Y1, t, status
end

function charikar_lp{T}(G::SparseMatrixCSC{T, Int64})

  set_default_solver(GurobiSolver())
  n = size(G, 1)
  X = Variable(n, n, Positive())
  Y = Variable(n, Positive())
  constraints = Convex.Constraint[]

  push!(constraints, sum(Y) <= 1)

  ii, jj = findn(G)
  for k=1:length(ii)
    crow = ii[k]
    ccol = jj[k]
    push!(constraints, X[crow, ccol] <= Y[crow])
    push!(constraints, X[crow, ccol] <= Y[ccol])
  end

  problem = maximize(trace(G * X), constraints);
  solve!(problem);
  t= problem.optval
  return problem, X.value, Y.value, t
end

function extract_primal_solution(Y)
  V = [j > 0.0 for j in Y]
  S = find(V)
  return S
end

# Uses LP primal solution to find dense common subgraph
#
#
function dcs_lp_extract(G::GraphVec; method::String = "DCS" )
  m = length(G)
  @assert m >= 1
  n = size(G{1}, 1)

  if method == "DCS"
    p, X, Y, t = dcs_lp(G)
  elseif method == "MIN"
    H = dcs_intersection_graph(G)
    p, X, Y, t =  charikar_lp(H)
  elseif method == "AVG"
    H = dcs_avg_graph(G)
    p, X, Y, t = charikar_lp(H)
  else
    error("Undefined method: $method") #  ('dcs', 'charikar', 'weighted')"
  end
  S = extract_primal_solution(Y)
  d = common_subgraph_density(S)
  return (S, d, t)
end

# function dcs_charikar_lp{T}(H::Array{SparseMatrixCSC{T, Int64}, 1})
#   m = length(H)
#   n = size(H[1], 1)
#   G = spzeros(T, n, n)

#   ii, jj = findn(H[1])
#   for k=1:length(ii)
#     cr = ii[k]
#     cc = jj[k]
#     G[cr, cc] = H[1][cr, cc]
#     for j=2:m
#       G[cr, cc] = min(G[cr, cc], H[j][cr, cc])
#     end
#   end
#   return charikar_lp(G)

  # end


text:
using Documenter
using LempelZiv

makedocs(
    format = :html,
    sitename = "LempelZiv.jl",
    pages = [
        "index.md",
        "About" => "README.md"
    ]
)


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

import Base: push!
import StatsBase: fit

export Sequential, search, push!, fit, optimize!

"""
    Sequential{T}

A simple exhaustive search index
"""
struct Sequential{T} <: Index
    db::Vector{T}
end

"""
    fit(::Type{Sequential}, db::Vector)

Creates a sequential-exhaustive index
"""
function fit(::Type{Sequential}, db::Vector)
    Sequential(db)
end

"""
    search(index::Sequential, dist::Function, q, res::KnnResult)

Solves the query evaluating ``dist(q,u) \\forall u \\in index`` against
"""
function search(index::Sequential, dist::Function, q, res::KnnResult)
    i::Int32 = 1
    d::Float64 = 0.0
    for obj in index.db
        d = dist(q, obj)
        push!(res, i, convert(Float32, d))
        i += 1
    end

    res
end

"""
   push!(index::Sequential, dist::Function, item)

Interts an item into the index
"""
function push!(index::Sequential, dist::Function, item)
    push!(index.db, item)
    length(index.db)
end

"""
    optimize!(index::Sequential{T}, dist::Function, recall::Float64)

Optimizes the index for the required quality
"""
function optimize!(index::Sequential{T}, dist::Function, recall::Float64) where T
    # do nothing for sequential
end


text:
using Demo
using Test

#2x + y

@testset "Demo.jl" begin
    @test my_f(1, 1) == 3
    @test my_f(1, 2) == 4
    @test my_ff(1, 2) == 5
end

max_stars_repo_name: aivodji/Yg

max_issues_repo_name: aivodji/Yg

max_forks_repo_name: aivodji/Yg


text:
using BilevelBenchmark
using Test

@testset "Bilevel unconstrained" begin
    D_ul = 5
    D_ll = 5
    for i = 1:8
        x, y = SMD_solutions(D_ul, D_ll, i)
        Fx, fy = SMD_leader(x, y, i), SMD_leader(x, y, i)
        @test  abs(Fx) + abs(fy) < 1e-10
    end
end

@testset "Bilevel Ranges" begin
    D_ul = 5
    D_ll = 5
    for i = 1:8
        bounds_ul, bounds_ll = SMD_ranges(D_ul, D_ll, i)
        x = bounds_ul[1,:] + (bounds_ul[2,:] - bounds_ul[1,:]) .* rand(D_ul)
        y = bounds_ll[1,:] + (bounds_ll[2,:] - bounds_ll[1,:]) .* rand(D_ul)

        Fx, fy = SMD_leader(x, y, i), SMD_leader(x, y, i)
        @test  abs(Fx) + abs(fy) > -Inf
    end
end

@testset "SMD" begin
    p, q, r, s = 3, 3, 2, 0
    for i = 1:8
        @test SMD_leader(ones(1000), ones(1000), i, p, q, r, s) ≈ SMD_leader(ones(1000), ones(1000), i, p, q, r, s)
        @test SMD_follower(ones(1000), ones(1000), i, p, q, r, s) ≈ SMD_follower(ones(1000), ones(1000), i, p, q, r, s)
    end
end

@testset "TP" begin
    for i = 1:10
        @test TP_test(i) >= 0

        # bounds_ul, bounds_ll = TP_ranges(i)

    end
end

@testset "PMM" begin
    for i = 1:6
        F, f = BilevelBenchmark.PMM_test(10, 10, i)
        @test abs(f + F) <= 1e-16
    end

    x = zeros(9)
    for fnum = 1:6
        y = PMM_Ψ(x, 9, fnum)
        if fnum < 7
            F = PMM_leader(x, y, fnum)
            f = PMM_follower(x, y, fnum)
        else
            F, _ = PMM_leader(x, y, fnum)
            f, _ = PMM_follower(x, y, fnum)
        end

        @test abs(F) <= 1e-16 &&  abs(f) <= 1e-16
    end
end

@testset "bilevel constrained" begin
    D_ul = 5
    D_ll = 5

    for i = 9:12
        p,q,r,s,lenG, leng = SMD_settings(D_ul, D_ll, i)

        bounds_ul, bounds_ll = SMD_ranges(D_ul, D_ll, i)
        x = bounds_ul[1,:] + (bounds_ul[2,:] - bounds_ul[1,:]) .* rand(D_ul)
        y = bounds_ll[1,:] + (bounds_ll[2,:] - bounds_ll[1,:]) .* rand(D_ul)

        F, G = SMD_leader(x, y, i)
        f, g = SMD_follower(x, y, i)

        @test length(G) == lenG
        @test length(g) == leng
        @test abs(F + sum(G)) >= 0
        @test abs(f + sum(g)) >= 0
    end
end

@testset "bilevel get problem" begin
    for benchmark in [:SMD, :PMM, :TP]
        for fnum = 1:6
            F, f, bounds_ul, bounds_ll, Ψ, z = get_problem(benchmark, fnum)
            x = bounds_ul[1,:]
            y = bounds_ll[1,:]
            Fxy = F(x, y)
            fxy = f(x,y)
            @test abs(z[1]) + abs(z[2]) >= 0
        end
    end

end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 0b19c578-762b-11eb-34b6-01e80cef1406
using Pkg, DrWatson

# ╔═╡ 3f814e9e-762b-11eb-1340-91617ca7b58a
begin
    #@quickactivate "StatisticalRethinkingStan"
    using StanSample
    using StatisticalRethinking
end

# ╔═╡ af1efcf4-779a-11eb-0d62-e3c63876a09b
md" ## Clip-09-22-24s.jl"

# ╔═╡ cb88768a-779b-11eb-1a4e-99f92ded29ca
begin
	y = [-1, 1]
	data = (N = length(y), y=y)
end

# ╔═╡ e37331ea-779b-11eb-1b02-1fdb61aae124
stan9_2 = "
data{
	int N;
    vector[N] y;
}
parameters{
    real alpha;
    real<lower=0> sigma;
}
model{
    real mu;
    sigma ~ exponential( 1e-04 );
    alpha ~ normal( 0 , 1000 );
    mu = alpha;
    y ~ normal( mu , sigma );
}
";

# ╔═╡ 3fafc7fe-762b-11eb-37da-254b6d75f5ca
begin
	m9_2s = SampleModel("m9.2s", stan9_2)
	rc9_2s = stan_sample(m9_2s; data)

	if success(rc9_2s)
		m9_2s_df = read_samples(m9_2s, :dataframe)
		PRECIS(m9_2s_df)
	end
end

# ╔═╡ 316b814a-7ebd-11eb-28a9-27ab6d02d9f9
read_summary(m9_2s)

# ╔═╡ b4d020b8-762d-11eb-167a-9fd79d948647
trankplot(m9_2s, :alpha)[1]

# ╔═╡ 728ad3aa-7eb7-11eb-2fea-634a82832501
trankplot(m9_2s, :sigma)[1]

# ╔═╡ 059b3ada-7eb8-11eb-10aa-459e3ec2bc55
begin
	chns9_2s = read_samples(m9_2s, :mcmcchains)
	Text(sprint(show, "text/plain", summarize(chns9_2s)))
end

# ╔═╡ 6fbfa790-7eb9-11eb-0360-ed6fc7fa11a3
plot(chns9_2s)

# ╔═╡ a93a71e4-7eb9-11eb-2a36-d33e7b88e2d6
stan9_3 = "
data {
	int N;
	vector[N] y;
}
parameters{
    real alpha;
    real<lower=0> sigma;
}
model{
    real mu;
    alpha ~ normal( 1 , 10 );
    sigma ~ exponential( 1 );
    mu = alpha;
    y ~ normal( mu , sigma );
}
";

# ╔═╡ b4ea8d58-7eb9-11eb-172f-3941a1b8fe98
begin
	m9_3s = SampleModel("m9.3s", stan9_3)
	rc9_3s = stan_sample(m9_3s; data)

	if success(rc9_3s)
		m9_3s_df = read_samples(m9_3s, :dataframe)
		PRECIS(m9_3s_df)
	end
end

# ╔═╡ ef7953e8-7eb9-11eb-3344-8d762d7e3584
trankplot(m9_3s, :alpha)[1]

# ╔═╡ 29b1d698-7eba-11eb-00aa-e95796579cc6
trankplot(m9_3s, :sigma)[1]

# ╔═╡ 39bb0f76-7eba-11eb-1bc2-293259d12dfc
begin
	chns9_3s = read_samples(m9_3s, :mcmcchains)
	Text(sprint(show, "text/plain", summarize(chns9_3s)))
end

# ╔═╡ 43b362c6-7eba-11eb-389e-fb53b3783e8d
plot(chns9_3s)

# ╔═╡ ca61c622-779a-11eb-1292-3f272953867f
md" ## End of clip-09-22-24s.jl"

# ╔═╡ Cell order:
# ╟─af1efcf4-779a-11eb-0d62-e3c63876a09b
# ╠═0b19c578-762b-11eb-34b6-01e80cef1406
# ╠═3f814e9e-762b-11eb-1340-91617ca7b58a
# ╠═cb88768a-779b-11eb-1a4e-99f92ded29ca
# ╠═e37331ea-779b-11eb-1b02-1fdb61aae124
# ╠═3fafc7fe-762b-11eb-37da-254b6d75f5ca
# ╠═316b814a-7ebd-11eb-28a9-27ab6d02d9f9
# ╠═b4d020b8-762d-11eb-167a-9fd79d948647
# ╠═728ad3aa-7eb7-11eb-2fea-634a82832501
# ╠═059b3ada-7eb8-11eb-10aa-459e3ec2bc55
# ╠═6fbfa790-7eb9-11eb-0360-ed6fc7fa11a3
# ╠═a93a71e4-7eb9-11eb-2a36-d33e7b88e2d6
# ╠═b4ea8d58-7eb9-11eb-172f-3941a1b8fe98
# ╠═ef7953e8-7eb9-11eb-3344-8d762d7e3584
# ╠═29b1d698-7eba-11eb-00aa-e95796579cc6
# ╠═39bb0f76-7eba-11eb-1bc2-293259d12dfc
# ╠═43b362c6-7eba-11eb-389e-fb53b3783e8d
# ╟─ca61c622-779a-11eb-1292-3f272953867f


text:
# # Instructions for Assignment05
#
# ## Introduction
#
# Note: this file is runnable in its current state,
# but is incomplete.
# You can run the file from the command line,
# or use the VS Code julia extension to run individual lines.

## To view additional instructions for this assignment, #jl
## see https://wellesley-bisc195.github.io/assignments/Assignment05/ #jl

## To view a rendered version of this document, #jl
## see https://wellesley-bisc195.github.io/assignments/Assignment05/#assignment05_code #jl

# ## Question 1 - Dubugging Practice
#
# At the end of [Chapter 8 of _Think Julia_](https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#deb08),
# there is an exercise (8-6) intended to help you debug a string traversal function.
# Below is the initial code with errors.
# Fix it up as described in the chapter.
#
# Note: you should always be completing the exercises from the book,
# by actually writing and executing the code,
# either in a file or in the REPL.

function isreverse(word1, word2)
    if length(word1) != length(word2)
        return false
    end
    i = firstindex(word1)
    j = lastindex(word2)
    while j >= 0
        j = prevind(word2, j)
        if word1[i] != word2[j]
            return false
        end
        i = nextind(word1, i)
    end
    true
end

# ## Question 2 - Reverse complement
#
# In Assignment03, you created a function called `complement()`,
# that took a DNA base as an argument,
# and returned its complement (`A <-> T` or `G <-> C`).
# This function does the same thing,
# using a slightly different method.
# Look it over and see if it makes sense.
# If it doesn't, try running some individual pieces in the REPL,
# But don't modify the code (or if you do, make sure you know how to get it back)

"""
    complement(base)

Get the DNA complement of the provided base:

    A <-> T
    G <-> C

Accepts uppercase or lowercase `String` or `Char`,
but always returns an uppercase `Char`.
If a valid base is not provided, the function throws an error.
"""
function complement(base)
    complements = Dict("A" => 'T',
                       "T" => 'A',
                       "G" => 'C',
                       "C" => 'G')

    base = uppercase(string(base))

    !(base in keys(complements)) && error("Invalid base $base")
    return complements[base]
end

# The "reverse complement" of a DNA sequence is valuable to know,
# since it represents the complementary sequence
# in the expected orientation (5' -> 3').
# If you haven't taken molecular biology and this sounds like gibberish,
# you may want to read up a bit on it [here](https://en.wikipedia.org/wiki/Complementarity_(molecular_biology))
#
# Now, write a function like `isreverse()` above,
# but that checks whether two sequences are the reverse AND complement
# of one another.
# Note one other difference here - rather than returning `false`
# if the sequences are different lengths, this function should error.

"""
    isreversecomplement(seq1, seq2)

Boolean function that checks whether two DNA seqences
are the reverse complement of one another, irrespective of capitalization.
Returns true if yes, false otherwise.

If any invalid bases are encountered,
or if sequences are different length, throws an error.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> isreversecomplement("aaatttcg", "cgaaattt")
    true

    julia> if isreversecomplement("C", "A")
               println("Yes!")
           else
               println("No!")
           end
    No!

    julia> isreversecomplement("TX", "AG")
    Error: Invalid base X

    julia> isreversecomplement("G", "CC")
    Error: Cannot compare sequuences of different length
"""
function isreversecomplement(seq1, seq2)
    ## your code here
end

# ## Question 3
#
# Now let's work on creating a function that _makes_ the reverse complement of a sequence,
# instead of just checking it.
# First, it's probably helpful to know that julia has a built in `reverse()` function
# that can take a `String` or `Vector` and return it with the indecies reversed.
#
# ```julia
# julia> reverse([4,3,2,1])
# 4-element Array{Int64,1}:
#  1
#  2
#  3
#  4
#
# julia> reverse("beep boop")
# "poob peeb"
# ```
#
# Now, create a function that takes a sequence
# and returns the reverse complement.
# The easiest way to do this is probably to perform a "map"
# of the `complement()` function to each character in the sequence,
# then reverse it (or reverse, then map).
# You can write the mapping by hand as in the book,
# or use `map()` as described in the lesson.

"""
    reverse_complement(sequence)

Takes a DNA sequence and returns the reverse complement
of that sequence.

Takes lowercase or uppercase sequences,
but always returns uppercase.

Examples
≡≡≡≡≡≡≡≡≡≡
    julia> reverse_complement("AAATTT")
    "AAATTT"

    julia> reverse_complement("GCAT")
    "ATGC"

    julia> rc = reverse_complement("TTGGG");

    julia> println(rc)
    CCCAA
"""
function reverse_complement(sequence)
    ## your code here
end

# ## Question 4a
#
# [Chapter 9](https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#reading_word_lists)
# begins to describe how to work with files.
# The file `words.txt` from _Think Julia_ is in the `data/` directory of the assignment repository.
# First, look at the file from the shell (not a julia REPL)
# using the `head` command.
#
# ```sh
# $ head data/words.txt
# ```
#
# By default, `head` prints the first 10 lines.
# How can you change this behavior to show the first 15 lines?
# (Hint: you may need a search engine).
# Put the *shell command* to show the first 15 words of `data/words.txt`
# in a String assigned to the variable `first15`

first15 = "#= put shell command here =#"

# ## Question 4b
#
# What about the shell command to view the last 15 lines of `data/words.txt`?
# Put this command in a string assigned to the variable `last15`

last15 = "#= put shell command here =#"

# ## Question 5
#
# Write a function that reads through the words in a file
# with a `for` loop and returns an array with only the words
# that are valid DNA sequences.
# Note: there may be one or more functions from previous Assignments
# that are useful to reuse here.

"""
    find_dna(file)

In a file with one word per line,
returns an array with all words that are valid DNA sequences
(that is, they contain only a, t, g, or c).

Words may be upper- or lowercase,
but returned array will contain only uppercase letters.

Example
≡≡≡≡≡≡≡≡≡
    julia> find_dna("data/words.txt") # assumes you're in Assignment05 directory
    13-element Array{Any,1}:
      "AA"
      "ACT"
      "ACTA"
      "AGA"
      "AT"
      "CAT"
      "GAG"
      "GAGA"
      "GAT"
      "TA"
      "TACT"
      "TAG"
      "TAT"
"""
function find_dna(file)
    ## 1. You probably want to make an empty array

    for line in eachline(file)
        ## if the line is valid dna, put it in the array
    end
end


text:
println("loading packages... ")

using ArgParse
using Common
using FileManager

println("packages OK")

"""
	get_potree_dirs(txtpotreedirs::String)

Return collection of potree directories.
"""
function get_potree_dirs(txtpotreedirs::String)
    if isfile(txtpotreedirs)
        return FileManager.readlines(txtpotreedirs)
    elseif isdir(txtpotreedirs)
        return [txtpotreedirs]
    end
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "potrees"
            help = "Potree projects directory"
            arg_type = String
            required = true
        "--output", "-o"
            help = "output file .las"
            arg_type = String
            required = true
    end

    return parse_args(s)
end

function main()
    args = parse_commandline()

    txtpotreedirs = args["potrees"]
    output = args["output"]

    sources = get_potree_dirs(txtpotreedirs)
    println("")
    println("== PARAMETERS ==")

    for i in eachindex(sources)
        println("source[$i]: $(sources[i])")
    end

    println("Output  =>  $output")

    # creo l'header: mi serve il bb che contiene tutto il modello e il numero di punti totali
    full_aabb = AABB()
    n_points_total = 0
    for potree in sources
        cloudmetadata = FileManager.CloudMetadata(potree)
        aabb = cloudmetadata.tightBoundingBox
        Common.update_boundingbox!(full_aabb, aabb)
        n_points_total += cloudmetadata.points
    end

    println("")
    println("AABB:")
    println("min: [$(full_aabb.x_min),$(full_aabb.y_min),$(full_aabb.z_min)]")
    println("max: [$(full_aabb.x_max),$(full_aabb.y_max),$(full_aabb.z_max)]")
    flush(stdout)
    println("")
    println("Points: $n_points_total")
    flush(stdout)
    mainHeader = FileManager.newHeader(
        full_aabb,
        "Potree2Las",
        FileManager.SIZE_DATARECORD,
        n_points_total
    )

    t = open(output, "w")
    write(t, FileManager.LasIO.magic(FileManager.LasIO.format"LAS"))
    write(t, mainHeader)
    flush(t)

    println("")
    println("== PROCESSING ==")

    # per ogni potree
    for potree in sources
        pointsProcessed = 0
        cloudmetadata = FileManager.CloudMetadata(potree)
        n_points = cloudmetadata.points
        trie = FileManager.potree2trie(potree)
        node_files = FileManager.get_all_values(trie)

        # apro il las
        for file in node_files
            s = open(file, "r")
            h, laspoints = FileManager.read_LAS_LAZ(file) # read file
            for laspoint in laspoints # read each point
                plas = FileManager.newPointRecord(
                    laspoint,
                    h,
                    FileManager.LasIO.LasPoint2,
                    mainHeader,
                )
                write(t, plas) # write this record on temporary file

                pointsProcessed += 1
                # if pointsProcessed % 1_000_000 == 0
                #     println("$pointsProcessed points processed of $n_points")
                # end
            end
            close(s)
            flush(t)

        end

    end

    close(t)
    FileManager.successful(true, splitdir(output)[1]; message="{\n\"points\": $n_points_total\n}")
end

@time main()


text:
using Test
using WiSARD

function tests()
    @testset "Encoding: Thermometer" begin
        @test thermometer([0.2, 0.5, 0.7], 5) == [thermometer(0.2, 5); thermometer(0.5, 5); thermometer(0.7, 5)] == Bool[1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0]
    end

    @testset "Simple Character Learning" begin
        wnn = WNN{Symbol, UInt}(3, 4)

        X = Bool[1 0 1
                 0 1 0
                 0 1 0
                 1 0 1]

        Y = Bool[1 0 1
                 0 1 0
                 0 1 0
                 0 1 0]

        Z = Bool[1 1 1
                 0 1 1
                 1 1 0
                 1 1 1]

        train!(wnn, :x, X)
        train!(wnn, :y, Y)
        train!(wnn, :z, Z)

        @test classify(wnn, X) == :x
        @test classify(wnn, Y) == :y
        @test classify(wnn, Z) == :z
    end

    @testset "Broadcasting" begin
        wnn = WNN{Symbol, UInt}(3, 4)

        x = [
            Bool[1 0 1
                 0 1 0
                 0 1 0
                 1 0 1],
            Bool[1 0 1
                 0 1 0
                 0 1 0
                 0 1 0],
            Bool[1 1 1
                 0 1 1
                 1 1 0
                 1 1 1],
        ]
        y = [:x, :y, :z]

        train!(wnn, y, x)

        @test classify(wnn, x) == y
    end
end

# Here we go!
tests()


text:
string_to_numbers(T::DataType,str::AbstractString) = parse.(T,split(str," ",keep=false))

function readneutral(inputfile::AbstractString)
    info("Reading $inputfile file")
    open(inputfile) do input

        NN = parse(Int,readline(input))
        #nodes = Vector{Vec2D{Float64}}(NN) #Pre-allocation a vector of Vec2D{Float64} with NN elemetns
        nodes = Vec2DArray{Float64}(NN) #Pre-allocation a vector of Vec2D{Float64} with NN elemetns

        for i=1:NN
            nodes[i] = Vec2D(string_to_numbers(Float64,readline(input))...)
        end

        NC = parse(Int,readline(input))

        cells = Vector{TriangleCell}(NC)
        bc_cell = Vector{UInt}(NC)

        for i=1:NC
            a = string_to_numbers(UInt,readline(input))
            bc_cell[i] = a[1]
            cells[i] = TriangleCell((a[2],a[3],a[4]))
        end

        NBF = parse(Int,readline(input))

        bf2n =  Vector{Face2Node{2}}(NBF)
        bc_face = Vector{UInt}(NBF)

        for i=1:NBF
            a = string_to_numbers(UInt,readline(input))
            bc_face[i] = a[1]
            bf2n[i] = Face2Node{2}((a[2],a[3]))
        end
        info("Number of nodes: $NN\nNumber of Cells: $NC\nNumber of Boundary Faces: $NBF\n\n")
        return nodes, cells, bc_cell, bf2n, bc_face
    end
end


text:
abstract type DynamicInverterComponent <: DynamicComponent end
abstract type Converter <: DynamicInverterComponent end
abstract type DCSource <: DynamicInverterComponent end
abstract type Filter <: DynamicInverterComponent end
abstract type FrequencyEstimator <: DynamicInverterComponent end
abstract type VSControl <: DynamicInverterComponent end

abstract type OuterControl <: DynamicInverterComponent end
abstract type ActivePowerControl <: DeviceParameter end
abstract type ReactivePowerControl <: DeviceParameter end

"""
Parameters of a Outer-Loop controller using a virtual inertia with VSM for active power controller
and a reactive power droop controller.
# Conmutable structor
```julia
VirtualInertiaQDroop(A, R)
```
# Arguments
*  `A`::Float64 : Active power controller using virtual inertia with VSM
*  `R`::Float64 : Reactive power controller using reactive power droop
"""
mutable struct VirtualInertiaQdroop{A <: ActivePowerControl, R <: ReactivePowerControl} <:
               OuterControl
    active_power::A
    reactive_power::R
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function VirtualInertiaQdroop(
    active_power::A,
    reactive_power::R,
    ext = Dict{String, Any}(),
) where {A <: ActivePowerControl, R <: ReactivePowerControl}
    VirtualInertiaQdroop(
        active_power,
        reactive_power,
        ext,
        vcat(active_power.states, reactive_power.states),
        active_power.n_states + reactive_power.n_states,
        InfrastructureSystemsInternal(),
    )
end

function VirtualInertiaQdroop(; active_power, reactive_power, ext = Dict{String, Any}())
    VirtualInertiaQdroop(active_power, reactive_power, ext)
end

get_active_power(value::VirtualInertiaQdroop) = value.active_power
get_reactive_power(value::VirtualInertiaQdroop) = value.reactive_power
get_ext(value::VirtualInertiaQdroop) = value.ext
get_states(value::VirtualInertiaQdroop) = value.states
get_n_states(value::VirtualInertiaQdroop) = value.n_states
get_internal(value::VirtualInertiaQdroop) = value.internal


text:
function gradient_descent_trace!(tr::OptimizationTrace,
                                 x::Vector,
                                 f_x::Real,
                                 gradient::Vector,
                                 iteration::Integer,
                                 store_trace::Bool,
                                 show_trace::Bool)
    dt = Dict()
    dt["g(x)"] = copy(gradient)
    dt["|g(x)|"] = norm(gradient, Inf)
    os = OptimizationState(copy(x), f_x, iteration, dt)
    if store_trace
        push!(tr, os)
    end
    if show_trace
        println(os)
    end
end

function gradient_descent(d::DifferentiableFunction,
                          initial_x::Vector;
                          tolerance::Real = 1e-8,
                          iterations::Integer = 1_000,
                          store_trace::Bool = false,
                          show_trace::Bool = false)

    # Keep a copy of the initial state of the system
    x = copy(initial_x)

    # Count the number of gradient descent steps we perform
    iteration = 0

    # Track calls to f and g!
    f_calls, g_calls = 0, 0

    # Allocate vectors for re-use by gradient calls
    n = length(x)
    gradient = Array(Float64, n)

    # Allocate vector for step direction
    p = Array(Float64, n)

    # Allocate vectors for line search
    ls_x = Array(Float64, n)
    ls_gradient = Array(Float64, n)

    # Maintain current value of f and g
    f_calls += 1
    g_calls += 1
    f_x = d.fg!(x, gradient)

    # Show trace
    tr = OptimizationTrace()
    if store_trace || show_trace
        gradient_descent_trace!(tr,
                                x,
                                f_x,
                                gradient,
                                iteration,
                                store_trace,
                                show_trace)
    end

    # Monitor convergence
    converged = false

    # Iterate until the norm of the gradient is within tolerance of zero
    while !converged && iteration < iterations
        # Increment the number of steps we've had to perform
        iteration = iteration + 1

        # Use a back-tracking line search to select a step-size
        #  Direction is always negative gradient
        for i in 1:n
            p[i] = -gradient[i]
        end
        step_size, f_update, g_update =
          backtracking_line_search!(d, x, p, ls_x, ls_gradient)
        f_calls += f_update
        g_calls += g_update

        # Move in the direction of the gradient
        for i in 1:n
            x[i] = x[i] - step_size * gradient[i]
        end

        # Update the function value and gradient
        f_calls += 1
        g_calls += 1
        f_x = d.fg!(x, gradient)

        # Assess convergence
        if norm(gradient, Inf) <= tolerance
            converged = true
        end

        # Show trace
        if store_trace || show_trace
            gradient_descent_trace!(tr,
                                    x,
                                    f_x,
                                    gradient,
                                    iteration,
                                    store_trace,
                                    show_trace)
        end
    end

    OptimizationResults("Gradient Descent w/ Backtracking Line Search",
                        initial_x,
                        x,
                        f_x,
                        iteration,
                        converged,
                        tr,
                        f_calls,
                        g_calls)
end


text:
function natural(x, y)
    k(x) = [occursin(r"\d+", s) ? parse(Int, s) : s
            for s in split(replace(x, r"\d+" => s->" $s "))]
    A = k(x); B= k(y)
    for (a, b) in zip(A, B)
        if !isequal(a, b)
            return typeof(a) <: typeof(b) ? isless(a, b) :
                   isa(a,Int) ? true : false
        end
    end
    return length(A) < length(B)
end

sort(["a1", "a2", "a10"]) |> display

sort(["a1", "a2", "a10"], lt=natural) |> display


text:
seed!(42)
L = 3
D = 10
N = 20
nSamples = 10
b = 5
x = rand(N, D)
y = rand(N)
@testset "HMC Sampling" begin
    i = HMCSampling(nBurnin = 0, thinning = 1)

    @test repr(i) == "Hamilton Monte Carlo Sampler"
    i = AGP.tuple_inference(i, L, D, N, b, [], [])

    @test AGP.getρ(i) == 1
    @test AGP.isStochastic(i) == false

    i = AGP.init_sampler!(i, L, N, nSamples, false)

    @test i.sample_store == zeros(Float64, nSamples, N, L)

    i.nIter = nSamples
    i = AGP.init_sampler!(i, L, N, nSamples, true)

    @test i.sample_store == zeros(Float64, 2*nSamples, N, L)

    i.nIter = 2*nSamples
    i = AGP.init_sampler!(i, L, N, 2, false)

    @test i.sample_store == zeros(Float64, 2, N, L)

    i.nIter = 1
    fs = [rand(N) for _ in 1:L]
    AGP.store_variables!(i, fs)
    @test i.sample_store[1, :, 1] == fs[1]
end


text:
function quandlget(id::String; order="des", rows=100, frequency="daily", transformation="none", from="", to="", format="TimeArray", auth_token="")

    # Create a dictionary with the Query arguments that we pass to get() function
    query_args = {"sort_order" => order, "rows" => rows, "collapse" => frequency, "transformation" => transformation, "trim_start" => from, "trim_end" => to}
    # delete definition above and uncomment below for v0.4
    #query_args = Dict{Any,Any}("sort_order" => order, "rows" => rows, "collapse" => frequency, "transformation" => transformation, "trim_from" => from, "trim_to" => to, auth_token => "")

    # Open the auth_token file and add to query_args if it exists
    auth_token = open(readall, Pkg.dir("Quandl/token/auth_token"))
    auth_token != "" ? query_args["auth_token"] = auth_token : nothing

    # Do not use rows if start or end date range specified
    if from != "" || to != ""
        delete!(query_args, "rows")
    end

    # Get the response from Quandl's API, using Query arguments (see Response.jl README)
    response = get("http://www.quandl.com/api/v1/datasets/$id.csv", query = query_args)

    if response.status != 200
        error("Dataset not found")
    end

    # Convert the response to the right DataType
    if format == "TimeArray"
        timearray(response)
    elseif format == "DataFrame"
        dataframe(response)
    else
        error("Invalid $format format. If you want this format implemented, please report an issue or submit a pull request.")
    end
end

# alias quandl/quandlget
quandl = quandlget

function quandlsearch(query::ASCIIString; page=1, results=20, format="DataFrame", quiet=false)

    query = replace(query, " ", "+")

    # Create a dictionary with the Query arguments that we pass to get() function
    query_args = {"query" => query, "page" => page, "per_page" => results}
    # delete definition above and uncomment below for v0.4
    #query_args = Dict{Any,Any}("query" => query, "page" => page, "per_page" => results)

    # Getting response from Quandl and parsing it
    response = get("http://www.quandl.com/api/v1/datasets.json", query = query_args)
    jsondict = JSON.parse(response.data)

    data = jsondict["docs"]
    totalcount = jsondict["total_count"]

    # Printing summary
    if !quiet
	    println("Returning $results results of $totalcount from page $page")
	end

    # Convert the response to the right DataType
    if format == "DataFrame"
        # Create an NA 1x5 DataFrame
        df = convert(DataFrame, ["" "" "" "" ""])

        # Constructiong the DataFrame
        for elem in data
            newline = convert(DataFrame, ["$(elem["source_code"])/$(elem["code"])" "$(elem["name"])" "$(elem["frequency"])" "$(elem["from_date"])" "$(elem["to_date"])"])
            df = vcat(df,newline)
        end

        # Naming the columns
        names!(df, [:Code, :Name, :Frequency, :From, :To])

        # Not returning the NA line
        return df[2:end,:]
    elseif format == "Dict"
        return data
    else
        error("Invalid $format format. Currently only DataFrame and Dict are supported. If you want this format implemented, please report an issue or submit a pull request.")
    end
end

function interactivequandl(query::ASCIIString; page="1", results="20", order="des", rows=100, frequency="daily", transformation="none", format="TimeArray", auth_token="")

	# Get search results
    searchres = quandlsearch(query, page = page, results = results, format="DataFrame", quiet=true)

    # Print results
    for i in 1:size(searchres, 1)
    	print_with_color(:yellow, string(i) * " ")
    	print_with_color(:blue, string(searchres[i, :Code]))
    	print_with_color(:black, "   From " * string(searchres[i, :From]) * " to " * string(searchres[i, :To]) * "\n")
    	print_with_color(:white, "    " * string(searchres[i, :Name]) * "\n")
    end

    # Print prompt
    print_with_color(:yellow, "\n==> ")
    print_with_color(:white, "Enter n° of the dataset to be downloaded | (N)ext page | (Q)uit\n")
    println("-------------------------------------------------------------------")
    print_with_color(:yellow, "==> ")

    # Receive input from user
    input = readline()[1:end-1] # Remove \n from end

    # Return desired output
    if input == "q" || input == "Q" || input == "" # Quit
    	return nothing
    elseif input == "n" || input == "N" # Next page
    	return interactivequandl(query, page = string(int(page) + 1), results = results, format = format, auth_token = auth_token)
    else  # Get and return result
	    return quandl(searchres[int(input), :Code],
	        order = order, rows = rows, frequency = frequency, transformation = transformation, format = format, auth_token = auth_token)
	end
end

function set_auth_token(token::String)

    # Check the validity of the token
	if length(token) != 20 && length(token) != 0
        error("Invalid Token: must be 20 characters long or be an empty string")
    end

    # Create the token directory if needed
    if !ispath(Pkg.dir("Quandl/token/"))
        run(`mkdir $(Pkg.dir("Quandl/token/"))`)
    end

    # Write to the file
    open(Pkg.dir("Quandl/token/auth_token"), "w") do token_file
        write(token_file, token)
    end

	return nothing
end


text:
using JuMP, Base.Test, AmplNLWriter

## Solve test problem 1 (Synthesis of processing system) in
 #  M. Duran & I.E. Grossmann, "An outer approximation algorithm for
 #  a class of mixed integer nonlinear programs", Mathematical
 #  Programming 36, pp. 307-339, 1986.  The problem also appears as
 #  problem synthes1 in the MacMINLP test set.
 #
 #  min   5 x4 + 6 x5 + 8 x6 + 10 x1 - 7 x3 -18 math.log(x2 + 1)
 #       - 19.2 math.log(x1 - x2 + 1) + 10
 #  s.t.  0.8 math.log(x2 + 1) + 0.96 math.log(x1 - x2 + 1) - 0.8 x3 >= 0
 #        math.log(x2 + 1) + 1.2 math.log(x1 - x2 + 1) - x3 - 2 x6 >= -2
 #        x2 - x1 <= 0
 #        x2 - 2 x4 <= 0
 #        x1 - x2 - 2 x5 <= 0
 #        x4 + x5 <= 1
 #        0 <= x1 <= 2
 #        0 <= x2 <= 2
 #        0 <= x3 <= 1
 #        x1, x2, x3 continuous
 #        x4, x5, x6 binary
 #
 #
 #  The solution is (1.30098, 0, 1, 0, 1, 0).
 ##

 # solver = AmplNLSolver(Ipopt.amplexe, ["print_level=0"])

@testset "example: jump_minlp" begin
    m = Model(solver=solver)
    x_U = [2,2,1]
    @variable(m, x_U[i] >= x[i=1:3] >= 0)
    @variable(m, y[i=4:6], Bin)

    @NLobjective(m, Min, 10 + 10*x[1] - 7*x[3] + 5*y[4] + 6*y[5] + 8*y[6] - 18*log(x[2]+1) - 19.2*log(x[1]-x[2]+1))
    @NLconstraint(m, 0.8*log(x[2] + 1) + 0.96*log(x[1] - x[2] + 1) - 0.8*x[3] >= 0)
    @NLconstraint(m, log(x[2] + 1) + 1.2*log(x[1] - x[2] + 1) - x[3] - 2*y[6] >= -2)
    @NLconstraint(m, x[2] - x[1] <= 0)
    @NLconstraint(m, x[2] - 2*y[4] <= 0)
    @NLconstraint(m, x[1] - x[2] - 2*y[5] <= 0)
    @NLconstraint(m, y[4] + 2*y[5]*0.5 <= 1)

    @test solve(m) == :Optimal

    if contains(getsolvername(solver), "ipopt")
        # Ipopt solves the relaxation
        @test isapprox(getvalue(x)[:], [1.14652, 0.546596, 1.0], atol=1e-5)
        @test isapprox(getvalue(y)[:], [0.27330, 0.299959, 0.0], atol=1e-5)
        @test isapprox(getobjectivevalue(m), 0.75928, atol=1e-5)
    else
        @test isapprox(getvalue(x)[:], [1.30098, 0.0, 1.0], atol=1e-5)
        @test isapprox(getvalue(y)[:], [0.0, 1.0, 0.0], atol=1e-5)
        @test isapprox(getobjectivevalue(m), 6.00975, atol=1e-5)
    end
end


text:
function baz()
  return 999
end


text:
### comparisons-tests.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file contains definition of comparison and test functions that support
# Measurement objects.
#
### Code:

import Base: ==, isless, <, <=, isnan, isfinite, isinf, isinteger

# Two measurements are equal if they have same value and same uncertainty.
# NB: Two measurements are egal (===) if they are exactly the same thing,
# i.e. if they share all the same fields, including the (hopefully) unique tag.
# If you need stricter equality use "===" instead of "==".
==(a::Measurement, b::Measurement) = (a.val==b.val && a.err==b.err)

# Comparison with Real: they are equal if the value of Measurement is equal to
# the number.  If you want to treat the Real like a measurement convert it with
# `Measurement'.
==(a::Measurement, b::Irrational) = a.val==float(b)
==(a::Measurement, b::Rational) = a.val==float(b)
==(a::Measurement, b::Real) = a.val==b
==(a::Irrational, b::Measurement) = float(a)==b.val
==(a::Rational, b::Measurement) = float(a)==b.val
==(a::Real, b::Measurement) = a==b.val

# Order relation is based on the value of measurements, uncertainties are ignored
<(a::Measurement, b::Measurement) = <(a.val, b.val)
<=(a::Measurement, b::Measurement) = <=(a.val, b.val)
isless(a::Measurement, b::Measurement) = isless(a.val, b.val)

isnan(a::Measurement) = isnan(a.val)
isfinite(a::Measurement) = isfinite(a.val)
isinf(a::Measurement) = isinf(a.val)
isinteger(a::Measurement) = isinteger(a.val)


text:
using ArtifactUtils, Artifacts, TOML

function add_all_artifacts()
    artifacts_file = abspath(joinpath(@__DIR__, "..", "Artifacts.toml"))

    metadata_hash = add_artifact!(
        artifacts_file,
        "metadata",
        "https://github.com/umd-otb/OpenTabularDataBenchmark/releases/download/v0.0.7-pre/data.toml.tar.gz",
        force=true,
    )

    metadata = TOML.parsefile(joinpath(artifact_path(metadata_hash), "data.toml"))

    for (name, data) in metadata
        println("Adding artifacts for $(name)")

        # add an artifact for the dataset file
        add_artifact!(
            artifacts_file,
            "$(name)-npz",
            string(data["data_location"], ".tar.gz"),
            lazy=true,
            force=true
        )

        # add an artifact for the extras file if present
        if "extras_location" in keys(data)
            add_artifact!(
                artifacts_file,
                "$(name)-extras",
                string(data["extras_location"], ".tar.gz"),
                lazy=true,
                force=true
            )
        end
    end
end

add_all_artifacts()


text:
# This file includes the Plugins module

abstract type AbstractPlugin end

# Define generic plugin functions.
function process end
function enable end
function disable end
function check end
function add end
function remove end

function search(rootpath::AbstractString, filename::AbstractString)
    paths = String[]
    for (root, dirs, files) in walkdir(rootpath)
        for file in files
            if occursin(filename, file)
                push!(paths, joinpath(root, file))
            end
        end
    end
    paths
end

const remote_repo_url = "https://imel.eee.deu.edu.tr/git/CausalPlugins.jl.git"

function add(name::AbstractString, url::AbstractString=remote_repo_url)
    startswith(name, ".") && error("Name of plugin should not start with `.`")
    startswith(name, "Plugins") && error("Name of plugin cannot be `Plugins`")
    startswith(".jl", name) || (name *= ".jl")

    repopath = joinpath("/tmp", "CausalPlugins", randstring())
    ispath(repopath) || mkpath(repopath)
    @info "Cloning avaliable plugins from $url"
    LibGit2.clone(url, repopath)
    @info "Done..."

    @info "Searching for $name in plugins repo."
    srcpath = search(joinpath(repopath, "src"), name)[1]
    if isempty(srcpath)
        error("$name could not be found in avaliable plugins")
    else
        dstdir = joinpath(@__DIR__, "additionals")
        dstpath = joinpath(dstdir, name)
        cp(srcpath, dstpath, force=true)
        include(dstpath)
        @info "$name is added to Causal.Plugins"
    end
end

macro def_plugin(ex)
    ex.head == :struct ||
        error("Invalid usage. The expression should start with `struct`.\n$ex")
    ex.args[2].head == :(<:) && ex.args[2].args[2] == :AbstractPlugin ||
        error("Invalid usage. The type should be a subtype of AbstractPlugin.\n$ex")
    quote
        Base.@kwdef $ex
    end |> esc
end

# # Includes essential plugins from Causal
# foreach(include, search(joinpath(@__DIR__, "essentials"), ".jl"))

# Includes additional plugins from Causal
foreach(include, search(joinpath(@__DIR__, "additionals"), ".jl"))

# Include plugins from user working directory.
user_plugins_path = joinpath(pwd(), "plugins")
ispath(user_plugins_path) && foreach(include, search(joinpath(user_plugins_path), ".jl"))


text:
using Documenter, LoopVectorization

makedocs(;
    modules=[LoopVectorization],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/chriselrod/LoopVectorization.jl/blob/{commit}{path}#L{line}",
    sitename="LoopVectorization.jl",
    authors="Chris Elrod",
    assets=[],
)

deploydocs(;
    repo="github.com/chriselrod/LoopVectorization.jl",
)


text:
using ExpressionTemplates
using PyPlot

Niters = 1000
NN = collect(1:.0625/2:7.25)

et_times = zeros(length(NN))
native_times = zeros(length(NN))
loop_times = zeros(length(NN))

function loop_kernel(w,a,x,b,y,c,z,d,N)
    res = zeros(N)
    for j = 1:N
        @inbounds res[j] = w*a[j] + x*b[j] + y*c[j] + z*d[j]
    end
end

for Ni = 1:length(NN)
    N = Int64(floor(10^NN[Ni]))

    println(N)

    a = rand(N)
    b = rand(N)
    c = rand(N)
    d = rand(N)

    w = rand()-.5
    x = rand()-.5
    y = rand()-.5
    z = rand()-.5

    tmp_et_times = zeros(Niters)
    tmp_nat_times = zeros(Niters)
    tmp_loop_times = zeros(Niters)
    for iter = 1:Niters

        #Disable garbage collection so it doesn't fire in the middle of a timed test
        gc_enable(false)
        tmp_et_times[iter] = @elapsed res = @et w*a + x*b + y*c + z*d
        tmp_nat_times[iter] = @elapsed res = w*a + x*b + y*c + z*d
        tmp_loop_times[iter] = @elapsed loop_kernel(w,a,x,b,y,c,z,d,N)
        gc_enable(true); #re-enable GC to clean up after the native array expression
    end

    et_times[Ni] = mean(tmp_et_times)
    native_times[Ni] = mean(tmp_nat_times)
    loop_times[Ni] = mean(tmp_loop_times)
end

println(et_times)
println(native_times)
println(loop_times)

table = zeros(length(NN),3)

table[:,1] = NN
table[:,2] = native_times./et_times
table[:,3] = native_times./loop_times

println()
println(table)

semilogx(floor(10.^table[:,1]), table[:,2])
semilogx(floor(10.^table[:,1]), table[:,3])
title("Speedup")
xlabel("N")
ylabel("Speedup = T_native / T")
legend(["Expression Templates", "Hand Loops"], "upper left")
savefig("speedup.svg")


text:
function generate(obj::Scheduler, arg0::List)
    return jcall(obj, "generate", SortedSet, (List,), arg0)
end

function get_builder(obj::AbstractScheduler)
    return jcall(obj, "getBuilder", MeasurementBuilder, ())
end

function get_selector(obj::AbstractScheduler)
    return jcall(obj, "getSelector", DatesSelector, ())
end

function init(obj::AbstractScheduler, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
#= ForwardDiff as a dependency is removed from this version
import MechanicalSketch
import MechanicalSketch: background, sethue, O, WI, HE, EM, FS
import MechanicalSketch: PALETTE, setfont, settext, setline, fontsize, color_with_lumin, color_with_lumin
import MechanicalSketch: dimension_aligned
import MechanicalSketch: generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: @import_expand, Quantity, @layer
import MechanicalSketch: place_image, set_scale_sketch, scale_to_pt, ∙
import MechanicalSketch: clamped_velocity_matrix, line_integral_convolution_complex, lic_matrix_current, LicSceneOpts
import MechanicalSketch: normalize_datarange, matrix_to_function, noise_for_lic, function_to_interpolated_function
import MechanicalSketch: Movie, Scene, animate, circle, line, Point, arrow, poly, draw_legend
import MechanicalSketch: scale_to_pt, x_y_iterators_at_pixels, Greys_9, BinLegend, scale_pt_to_unit
import Interpolations:   Extrapolation

let

if !@isdefined m²
    @import_expand(m, s)
end

include("test_functions_32.jl")

# Scaling and placement
physwidth = 10.0m
physheight = 4.0m
totheight = physheight * 1.1 * 3
set_scale_sketch(totheight, HE)
totwidth = totheight * WI / HE
Δx = totwidth / 2 - 0.55 * physwidth
Δy = physheight * 1.1

# Reused velocity field from earlier tests
velocity_matrix = clamped_velocity_matrix(ϕ_32; physwidth, physheight, cutoff = 0.5m/s);
# One complex matrix: Phase and amplitude for the visualization. This can generate cyclic movies
complex_convolution_matrix = convolute_image_32(velocity_matrix) # 15.071 s (364518 allocations: 61.17 MiB)
# The distribution of phase angles (complex argument) ought to be flat between -π and π. Let's check that:
n_pixels = round(Int, scale_to_pt(physwidth))
phasetop, binwidth, relfreq = phase_histog(complex_convolution_matrix, n_pixels)
histogrampoints = Point.(collect(phasetop * n_pixels / 2π), -EM .* relfreq);

# We'll also prepare a uniform flow field, constant velocity 0.5 m/s.
fxy_unif(x, y) = (0.5, 0.0)m∙s⁻¹
# Phase and amplitude for the visualization. This can generate cyclic movies
complex_convolution_matrix_uniform = convolute_image_32(fxy_unif; physwidth, physheight)

# We'll also plot a vertically uniform flow field, horizontally increasing velocity from 0 to 0.5 m/s.
fxy_lin(x, y) = (0.5 * (x / physwidth + 0.5), 0.0)m∙s⁻¹
# Phase and amplitude for the visualization. This can generate cyclic movies
complex_convolution_matrix_linear = convolute_image_32(fxy_lin; physwidth, physheight)

legend = BinLegend(;maxlegend = 0.2, minlegend = -1.0, noofbins = 256,
                       colorscheme = reverse(Greys_9),
                       nan_color = color_with_lumin(PALETTE[1], 80), name = Symbol("Value{Float64}"))

# Define scene functions (parts of each image)

# Rectangular flow field plot including a visual frame counter
function plot_flowfield(scene, framenumber)
    @assert scene.opts isa LicSceneOpts string(typeof(scene.opts))
    currentvals = lic_matrix_current(scene, framenumber)
    place_image(scene.opts.O, scene.opts.legend.(currentvals))
end

# Histogram plot, overlay to the flow field
function plot_histogram(scene, framenumber)
    @assert scene.opts isa LicSceneOpts
    histogrampoints = scene.opts.data
    @assert histogrampoints isa Array{Point,1}
    n_pixel_x = size(histogrampoints)[1]
    @layer let
        O = scene.opts.O
        settext("Relative frequency of phase angle",
            O + (-5EM,  -EM ), markup=true)
        settext("Angle values, min = -π, max = π",
            O + (n_pixel_x / 2, 1.5EM), markup=true, halign = "right")
        sethue(PALETTE[3])
        arrow(O, O + (0.0,  -EM))
        arrow(O, O + (n_pixel_x * 1.1 / 2, 0.0))
        poly(O .+ histogrampoints, :stroke)
        settext("π / 2 ", O + (n_pixels / 4,  0.0), markup=true)
     end
end

# Overlay to the flow field
function plot_frame_indicator(scene, framenumber)
    @assert scene.opts isa LicSceneOpts
    @layer let
        O = scene.opts.O
        sethue(PALETTE[9])
        setline(3)
        circle(O + (3EM, -3EM), 1FS, :stroke)
        frames_per_cycle = scene.opts.cycle_duration ∙ scene.opts.framerate
        # Where we are in the repeating cycle, [0, 1 - frame_duration]
        normalized_time = framenumber / frames_per_cycle
        θ = 2π ∙ normalized_time
        line(O + (3EM, -3EM), O + (3.0EM, -3.0EM) + (cos(θ), sin(θ)) .* FS, :stroke)
        settext(string(framenumber), O + (4EM, -4EM))
    end
end

function backdrop(scene, framenumber)
    setfont("Calibri", FS)
    fontsize(FS)
    O = scene.opts.O
    background(color_with_lumin(PALETTE[6], 30))
    sethue(PALETTE[3])
    Δt = scene.opts.cycle_duration
    frames_per_cycle = Δt ∙ scene.opts.framerate
    # Where we are in the repeating cycle, [0, 1 - frame_duration]
    normalized_time = framenumber / frames_per_cycle
    t = Δt ∙ normalized_time
    settext("Uniform velocity 0.5 m/s:\nCheck that one pixel moves \n$(Δt * 0.5m/s) over $Δt",
        O + (-WI / 2 + EM, 0.0), markup=true)
    dimension_aligned(O + (-3.0m, 2.0m), O + (-2.0m, 2.0m),
        fromextension = (0, EM), toextension = (0, EM), offset = -EM)
    dimension_aligned(O + (-4.0m, 2.0m), O + (-3.0m, 2.0m),
        fromextension = (0, EM), toextension = (0, EM), offset = -EM)
    dimension_aligned(O + (-5.0m, 2.0m), O + (-4.0m, 2.0m),
        fromextension = (0, EM), toextension = (0, EM), offset = -EM)

    settext("t = $(round(s,t, digits=2))",
        O + (-WI / 2 + EM, -4EM), markup=true)

    for i = 0:9
        x = t ∙ 0.5m/s + i ∙ 1.0m -5.0m
        line(O + (x, 1.1 * 2.0m), O + (x, 2.0m), :stroke)
    end
    draw_legend(O + (5.0m, 2.0m) + (EM, 0), scene.opts.legend, max_vert_height = 4.0m)
end

Δt = 2.0s
endframe = Int(floor(Δt ∙ (30∙s⁻¹)) - 1)
movie = Movie(WI, HE, "Three flow fields", 15:15)

scenes = [
    Scene(movie, backdrop, 0:endframe,             optarg = LicSceneOpts(Δt, O; legend)),
    Scene(movie, plot_flowfield, 0:endframe,       optarg = LicSceneOpts(Δt,  O + (Δx, -Δy); phase_magnitude_matrix = complex_convolution_matrix, legend)),
    Scene(movie, plot_histogram, 0:endframe,       optarg = LicSceneOpts(Δt, O + (Δx, -Δy ), data = histogrampoints)),
    Scene(movie, plot_frame_indicator, 0:endframe, optarg = LicSceneOpts(Δt, O + (Δx, -Δy ))),
    Scene(movie, plot_flowfield, 0:endframe,       optarg = LicSceneOpts(Δt,  O;  phase_magnitude_matrix = complex_convolution_matrix_uniform, legend)),
    Scene(movie, plot_flowfield, 0:endframe,       optarg = LicSceneOpts(Δt, O + (Δx, Δy); phase_magnitude_matrix = complex_convolution_matrix_linear, legend))
    ];

animate(movie, scenes,
    creategif = true,
    pathname = joinpath(@__DIR__, "test_32.gif"),
    framerate = 30)
set_scale_sketch(m)
end # let
=#


text:
################################################################################
# Toy Example
################################################################################
using Plots
T = Float64

# Define the dynamics of the system
p = 3 # Number of players
# model = DoubleIntegratorGame(p=p) # game with 3 players with unicycle dynamics
model = UnicycleGame(p=p) # game with 3 players with unicycle dynamics
model = BicycleGame(p=p) # game with 3 players with unicycle dynamics
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(1*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2, 0.0,0,0]),
      SVector{model.ni[2],T}([2, 0.0,0,0]),
      SVector{model.ni[3],T}([2, 0.0,0,0]),
      # SVector{model.ni[4],T}([3,+0.8,0,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]

# # Objectives of the game
game_obj = GameObjective(Q,R,xf,uf,N,model)
# radius = 0.5*ones(p)
# μ = 4.0*ones(p)
# μ = 4.0*[1.0, 0.1]
# μ = 4.0*[0.1, 1.0]
# add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = GameConstraintValues(probsize)
# Add collision avoidance
radius = 0.19
add_collision_avoidance!(game_con, radius)

# # Add control bounds
# u_max =  5*ones(SVector{m,T})
# u_min = -5*ones(SVector{m,T})
# add_control_bound!(game_con, u_max, u_min)

# # Add wall constraint
# walls = [Wall([-1.0,-1.75], [1.0,0.25], [1.,-1.]/sqrt(2))]
# add_wall_constraint!(game_con, walls)

# # Add circle constraint
# xc = [1., 2., 3.]
# yc = [1., 2., 3.]
# radius = [0.1, 0.2, 0.3]
# add_circle_constraint!(game_con, xc, yc, radius)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.0, 0.0, 0.0, #0.0,
   -0.4, 0.0, 0.4, #0.6,
    0.0, 0.0, 0.0, #0.0,
    0.0, 0.0, 0.0, #0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-5
opts.α_dual = 1.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-5
opts.ϵ_sta = 1e-5
opts.ϵ_con = 1e-5
opts.ϵ_opt = 1e-5
opts.regularize = true
opts.αx_dual = ones(p)

# s = 0.5
# opts.αx_dual = [sqrt(s), 1/sqrt(s)]
# Define the game problem
prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
@time newton_solve!(prob)
# @profiler newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

function update_λ!(ascore::ActiveSetCore, prob::GameProblem, Δλcol::AbstractVector)
    probsize = prob.probsize
    N = probsize.N
    p = probsize.p
    px = probsize.px

    stamp = CStamp()
    for i = 1:p
        for j ∈ setdiff(1:p,i)
            conval, v = get_collision_conval(prob.game_con, i, j)
            if v
                for (l,k) in enumerate(conval.inds)
                    stampify!(stamp, :h, :col, i, j, k)
                    if conval.active[l][1] == 1
                        conval.λ[l] += Δλcol[horizontal_idx(ascore, stamp) .- probsize.S]
                        conval.λ[l] = max.(0, conval.λ[l])
                    end
                end
            end
        end
    end
    return nothing
end

function subspace_dimension(prob::GameProblem; α::T=2e-3) where {T}
    prob_ = deepcopy(prob)
    probsize = prob.probsize
    ascore = ActiveSetCore(probsize)
    update_nullspace!(ascore, prob, prob.pdtraj)
    ns = length(ascore.null.vec)
    val = 0
    subspace = Vector{PrimalDualTraj}([deepcopy(prob.pdtraj)])
    plt = plot(legend=false)
    anim = @animate for l = 1:ns
        prob_ = deepcopy(prob)
        prob_.opts.shift = 0
        prob_.opts.α_dual = 0.0
        prob_.opts.αx_dual = 1e-2*ones(p)
        # prob_.opts.αx_dual = 0.0*ones(p)
        prob_.opts.dual_reset = false
        prob_.opts.ρ_0 = prob_.pen.ρ[1]
        prob_.opts.ρ_max = prob_.pen.ρ[1]
        prob_.opts.inner_print = false
        prob_.opts.outer_print = false
        prob_.opts.inner_iter = 20
        # prob_.opts.inner_iter = 20
        prob_.opts.outer_iter = 20
        # prob_.opts.outer_iter = 3
        prob_.opts.reg_0 = 1e-5
        # ϵ = 1e-6
        ϵ = 1e-5
        prob_.opts.ϵ_opt = ϵ
        prob_.opts.ϵ_dyn = ϵ
        prob_.opts.ϵ_con = ϵ
        prob_.opts.ϵ_sta = ϵ
        residual!(prob_)
        set_traj!(prob_.core, prob_.Δpdtraj, α*ascore.null.Δtraj[l])
        update_traj!(prob.pdtraj, prob_.pdtraj, 1.0, prob_.Δpdtraj)
        update_λ!(ascore, prob_, α*ascore.null.Δλ[l])
        newton_solve!(prob_)
        plot!(prob_.stats)

        opt = mean(abs.(prob_.stats.opt_vio[end].vio))
        sta = mean(abs.(prob_.stats.sta_vio[end].vio))
        dyn = mean(abs.(prob_.stats.dyn_vio[end].vio))
        con = mean(abs.(prob_.stats.con_vio[end].vio))
        opts = prob_.opts
        if opt < opts.ϵ_opt && sta < opts.ϵ_sta && dyn < opts.ϵ_dyn && con < opts.ϵ_con
            val += 1
            plot!(prob_.model, prob_.pdtraj.pr, plt=plt)
            push!(subspace, deepcopy(prob_.pdtraj))
        end
    end
    @show ns
    @show val
    path = joinpath("~/Documents/nullspace_dim_p3_$(string(typeof(model).name)).gif")
    gif(anim, path, fps = 10)
    return ascore, subspace
end

function pca(prob::GameProblem, subspace::Vector{PrimalDualTraj})
    probsize = prob.probsize
    N = probsize.N
    n = probsize.n
    m = probsize.m
    p = probsize.p
    px = probsize.px
    S = probsize.S
    M = length(subspace)-1

    V = []
    v = zeros(S)
    v_ref = zeros(S)
    get_traj!(prob.core, subspace[1], v_ref)
    for k = 1:M
        get_traj!(prob.core, subspace[k+1], v)
        v = v - v_ref
        push!(V, deepcopy(v))
    end

    # R = (N-1)*(n+m)
    pr_mask = primals_mask(prob.core, probsize)
    R = length(pr_mask)
    W = [v[pr_mask] for v in V]
    # W = [v[1:S] for v in V]
    m = sum(W) ./ M
    Wc = [w .- m for w in W]
    cov = sum([wc*wc' for wc in Wc])/M
    vecs = eigvecs(cov)
    vals = eigvals(cov)
    eigs = [(vecs[:,i], vals[i]) for i=1:R]
    # eigs = [(vecs[:,i], vals[i]) for i=1:S]
    sort!(eigs, by = x -> x[2], rev=true)
    vecs = [e[1] for e in eigs[1:M]]
    vals = [e[2] for e in eigs[1:M]]

    @show round.(log.(10, abs.(vals)), digits=3)

    plt = plot(legend=false)
    pdtraj_ref = subspace[1]
    plot!(prob.model, pdtraj_ref.pr, plt=plt)

    pdtraj_1 = deepcopy(pdtraj_ref)
    v1 = zeros(S)
    v1[pr_mask] = vecs[1]
    set_traj!(prob.core, pdtraj_1, v1)

    pdtraj_2 = deepcopy(pdtraj_ref)
    v2 = zeros(S)
    v2[pr_mask] = vecs[2]
    set_traj!(prob.core, pdtraj_2, v2)

    pdtraj_3 = deepcopy(pdtraj_ref)
    v3 = zeros(S)
    v3[pr_mask] = vecs[3]
    set_traj!(prob.core, pdtraj_3, v3)

    x = [[state(pdtraj_ref.pr[k])[px[i][1]] for k=2:N] for i=1:p]
    y = [[state(pdtraj_ref.pr[k])[px[i][2]] for k=2:N] for i=1:p]
    u1 = [[state(pdtraj_1.pr[k])[px[i][1]]  for k=2:N] for i=1:p]
    v1 = [[state(pdtraj_1.pr[k])[px[i][2]]  for k=2:N] for i=1:p]
    u2 = [[state(pdtraj_2.pr[k])[px[i][1]]  for k=2:N] for i=1:p]
    v2 = [[state(pdtraj_2.pr[k])[px[i][2]]  for k=2:N] for i=1:p]
    u3 = [[state(pdtraj_3.pr[k])[px[i][1]]  for k=2:N] for i=1:p]
    v3 = [[state(pdtraj_3.pr[k])[px[i][2]]  for k=2:N] for i=1:p]
    β = 2.0
    for i = 1:p
        quiver!(x[i],y[i],quiver=(β*u1[i],β*v1[i]), linewidth=1.5, color=:cyan)
        quiver!(x[i],y[i],quiver=(β*u2[i],β*v2[i]), linewidth=1.5, color=:red)
        quiver!(x[i],y[i],quiver=(β*u3[i],β*v3[i]), linewidth=1.5, color=:black)
    end
    # y = rand(1:10,N)
    # u = rand(N)
    # v = rand(N)
    # scatter(x,y)

    display(plt)

    return vals, pdtraj_1
end

function primals_mask(core::NewtonCore, probsize::ProblemSize)
    N = probsize.N
    inds = Vector{Int}()
    stamp = HStamp()
    for k = 1:N-1
        # States
        stampify!(stamp, :x, 1, k+1)
        ind = horizontal_idx(core, stamp)
        inds = vcat(inds, Vector(ind))
        # Controls
        for i = 1:p
            stampify!(stamp, :u, i, k)
            ind = horizontal_idx(core, stamp)
            inds = vcat(inds, Vector(ind))
        end
    end
    return inds
end

function display_eigvals(vals)
    plt = plot(legend=false)
    plot!(vals,
        linetype=:steppost,
        linewidth=3.0,
        title="PCA on the Nash Equilibrium subspace \n around a nominal NE trajectory.",
        xlabel="eigval_index",
        ylabel="eigval_magnitude")
    display(plt)
    return nothing
end

prob_copy = deepcopy(prob)
# ascore, subspace = subspace_dimension(prob_copy, α=4e-3) # Unicycle, # Bicycle
subspace
vals, pdtraj_1 = pca(prob, subspace)
display_eigvals(vals)

pdtraj_1

pr_mask = primals_mask(prob.core, probsize)


text:
@testset "Bytes - Integer Convertions" begin
        tests = ((11, [0x0b]),
                 (18145,[0x46, 0xe1]),
                 (big(18)^256, [0x0b, 0x51, 0xea, 0xa7, 0x75, 0x9a, 0x45, 0x93, 0x57, 0xc0, 0x6f, 0x65, 0x48, 0x7c, 0xe1, 0x8f, 0x23, 0xcf, 0x2b, 0x52, 0xbb, 0xab, 0x34, 0x98, 0x3e, 0x22, 0xb2, 0xa5, 0xca, 0xdd, 0xe3, 0x3f, 0xcb, 0x0f, 0x1a, 0x5a, 0xbf, 0xca, 0x90, 0xc4, 0xf8, 0x13, 0x28, 0x05, 0x2d, 0x31, 0x1f, 0xe7, 0x77, 0x93, 0x41, 0x7d, 0xd7, 0x70, 0x5c, 0x83, 0x09, 0x8e, 0xc0, 0xc7, 0x32, 0xe0, 0xf9, 0xc8, 0x7f, 0xe3, 0xaa, 0x65, 0xac,0x28, 0x94, 0x83, 0x76, 0x0f, 0xbd, 0x2a, 0x7e, 0x53, 0x35, 0xbd, 0x07, 0x5f, 0x4d, 0x13, 0xeb, 0xbd, 0x76, 0x79, 0xef, 0x1b, 0x43, 0x06, 0x89, 0x0c, 0x1d, 0x66, 0x0d, 0x12, 0x76, 0xf1, 0xe8, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]))
        @testset "int2bytes" begin
            for t in tests
                @test int2bytes(t[1]) == t[2]
            end
        end
        @testset "bytes2int" begin
            for t in tests
                @test bytes2int(t[2]) == t[1]
            end
        end
end


text:
"""
    FVectorvariate

Structure that subtypes `VariateForm`, indicating that a single sample
is a vector from a fiber of a vector bundle.
"""
struct FVectorvariate <: VariateForm end

"""
    FVectorSupport(space::AbstractManifold, VectorBundleFibers)

Value support for vector bundle fiber-valued distributions (values from a fiber of a vector
bundle at a `point` from the given manifold).
For example used for tangent vector-valued distributions.
"""
struct FVectorSupport{TSpace<:VectorBundleFibers,T} <: ValueSupport
    space::TSpace
    point::T
end

"""
    FVectorDistribution{TSpace<:VectorBundleFibers, T}

An abstract distribution for vector bundle fiber-valued distributions (values from a fiber
of a vector bundle at point `x` from the given manifold).
For example used for tangent vector-valued distributions.
"""
abstract type FVectorDistribution{TSpace<:VectorBundleFibers,T} <:
              Distribution{FVectorvariate,FVectorSupport{TSpace,T}} end

"""
    MPointvariate

Structure that subtypes `VariateForm`, indicating that a single sample
is a point on a manifold.
"""
struct MPointvariate <: VariateForm end

"""
    MPointSupport(M::AbstractManifold)

Value support for manifold-valued distributions (values from given
[`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M`).
"""
struct MPointSupport{TM<:AbstractManifold} <: ValueSupport
    manifold::TM
end

"""
    MPointDistribution{TM<:AbstractManifold}

An abstract distribution for points on manifold of type `TM`.
"""
abstract type MPointDistribution{TM<:AbstractManifold} <:
              Distribution{MPointvariate,MPointSupport{TM}} end

"""
    support(d::FVectorDistribution)

Get the object of type `FVectorSupport` for the distribution `d`.
"""
function Distributions.support(::T) where {T<:FVectorDistribution}
    return error("support not implemented for type $T")
end

function uniform_distribution(M::AbstractManifold)
    return uniform_distribution(M, allocate_result(M, uniform_distribution))
end

"""
    Random.rand(M::AbstractManifold, [d::Integer]; vector_at=nothing)
    Random.rand(rng::AbstractRNG, M::AbstractManifold, [d::Integer]; vector_at=nothing)

Generate a random point on manifold `M` (when `vector_at` is `nothing`) or a tangent
vector at point `vector_at` (when it is not `nothing`).

Optionally a random number generator `rng` to be used can be specified. An optional integer
`d` indicates that a vector of `d` points or tangent vectors is to be generated.

!!! note

    Usually a uniform distribution should be expected for compact manifolds and a
    Gaussian-like distribution for non-compact manifolds and tangent vectors, although it is
    not guaranteed. The distribution may change between releases.

    `rand` methods for specific manifolds may take additional keyword arguments.

"""
Random.rand(M::AbstractManifold)
function Random.rand(M::AbstractManifold, d::Integer; kwargs...)
    return [rand(M; kwargs...) for _ in 1:d]
end
function Random.rand(rng::AbstractRNG, M::AbstractManifold, d::Integer; kwargs...)
    return [rand(rng, M; kwargs...) for _ in 1:d]
end
function Random.rand(M::AbstractManifold; kwargs...)
    p = allocate_result(M, rand)
    rand!(M, p; kwargs...)
    return p
end
function Random.rand(rng::AbstractRNG, M::AbstractManifold; kwargs...)
    p = allocate_result(M, rand)
    rand!(rng, M, p; kwargs...)
    return p
end


text:
"""

    `comp_smp_ex!(RS)`

    Stores sample average of E[X].

    # Examples
    ```julia-repl
    julia> CpelNano.comp_smp_ex!(rs)
    ```
"""
function comp_smp_ex!(rs::RegStruct)::Nothing

    # Loop over groups
    ex = fill(NaN,rs.N)
    @inbounds for l=1:rs.L
        aux = 0.0
        n_data = 0.0
        # Loop over calls
        @inbounds for m=1:length(rs.calls)
            # Get call
            call = rs.calls[m][l]
            if call.obs
                # If observed add contribution
                aux += call.log_pyx_m>call.log_pyx_u ? 1.0 : -1.0
                n_data += 1.0
            end
        end
        # Store if any data point
        if n_data>0.0
            ex[rs.cpg_grps[l].cpg_ind] .= aux/n_data
        end
    end

    # Store in struct
    rs.exps.ex = ex

    # Return
    return nothing

end
"""

    `comp_smp_exx!(RS)`

    Stores sample average of E[XX].

    # Examples
    ```julia-repl
    julia> CpelNano.comp_smp_exx!(rs)
    ```
"""
function comp_smp_exx!(rs::RegStruct)::Nothing

    # Loop over groups
    exx = fill(NaN,rs.N-1)
    @inbounds for l=1:(rs.L-1)
        aux = 0.0
        n_data = 0.0
        # Loop over calls
        @inbounds for m=1:length(rs.calls)
            # Get call
            call1 = rs.calls[m][l]
            call2 = rs.calls[m][l+1]
            # Interaction between groups
            if call1.obs && call2.obs
                # If observed add contribution
                aux1 = call1.log_pyx_m>call1.log_pyx_u ? 1.0 : -1.0
                aux2 = call2.log_pyx_m>call2.log_pyx_u ? 1.0 : -1.0
                aux += aux1*aux2
                # Increase counter
                n_data += 1.0
            end
        end
        # Store between groups if any data point
        if n_data>0.0
            exx[rs.cpg_grps[l].cpg_ind[end]] = aux/n_data
        end
        # Store within groups if not singleton
        if length(rs.cpg_grps[l].cpg_ind)>1
            exx[rs.cpg_grps[l].cpg_ind[1:(end-1)]] .= 1.0
        end
    end

    # Take care of last group if not singleton
    if length(rs.cpg_grps[end].cpg_ind)>1
        exx[rs.cpg_grps[end].cpg_ind[1:(end-1)]] .= 1.0
    end

    # Store in struct
    rs.exps.exx = exx

    # Return
    return nothing

end
"""

    `comp_smp_exps!(RS)`

    Stores sample averages of E[X] and E[XX].

    # Examples
    ```julia-repl
    julia> CpelNano.comp_smp_exps!(rs)
    ```
"""
function comp_smp_exps!(rs::RegStruct)::Nothing

    # Compute expectations
    comp_smp_ex!(rs)
    comp_smp_exx!(rs)

    # Return
    return nothing

end


text:
using Revise
using Lifting
import Lifting:
    Lifting_Programmes, Lifting_Aux, Lifting_Exercise_Names, Lifting_Progressions
cd(@__DIR__)

name = "nSunsCAP3_OHP_6Day_LP"

# Uncomment to print whole programme.
println(Lifting_Programmes[name])
# Prints the days provided in the second argument.
println(Lifting_Programmes[name], [2, 6, 9])

println(Lifting_Programmes[name], length(Lifting_Programmes[name].days))

write("yugenstronk.csv", Lifting_Programmes["DCG_YugeStr"]; log = true)
#=
# Writes programme to test.cvt
write("test.csv", Lifting_Programmes[name]; log = false)
# Writes days provided in the third argument programme to test.cvt
write("test.csv", Lifting_Programmes[name], [9,15,21]; log = false)
# Write csv file for logging reps, weight and performance for all the lifts. Open the file up in google docs and be log your stuff there.
write("test.csv", Lifting_Programmes[name]; log = true) =#
# Makes a complete copy of the lifting programme with the given name.
prog = deepcopy(Lifting_Programmes[name])

# Loads log file data. Returns dictionaries whose keys are the exercise names.
keyArr, date, day1, Δdays, reps, wght, rpe = loadLogFile(prog)

# Calculates training maxes depending on the weight and reps provided but accounting for intensity defined in the progression.
tm, change = calcTrainingMaxLogs(prog, keyArr, reps, wght)

tm
change
# Plots lifts individually and uses a cubic spline to interpolate values.
figs = plotData(prog, keyArr, Δdays, tm; xlabel = "Days", ylabel = "Weight", lw = 3)
for (i, fig) in enumerate(figs)
    scatterData!(
        fig,
        prog,
        keyArr[i],
        Δdays,
        tm;
        shape = :circle,
        markersize = 5,
        label = "",
    )
end
figs
display.(figs)

# Calculates training maxes depending on the weight and intensities used in the progression and the reps provided.
updateMaxes(prog, keyArr, reps)

# Same as above but also updates the programme with the new values. Run this a couple of times and see how the training maxes change by running the one above.
updateMaxes!(prog, keyArr, reps)

write("nSunsCAP3_RPE.csv", Lifting_Programmes["nSunsCAP3_RPE"])
println(Lifting_Programmes["nSunsCAP3_RPE"], 22)

Lifting_Programmes["DCG_CAP4"]

write("DCG_CAP4.csv", Lifting_Programmes["DCG_CAP4"])
write("DCG_CAP4.csv", Lifting_Programmes["DCG_CAP4"]; log = true)
println(Lifting_Programmes["DCG_CAP4"])
println(Lifting_Programmes["DCG_CAP4"], [1,8,16,23]) # Bench T11 & Bench T21
println(Lifting_Programmes["DCG_CAP4"], [2,9,17,24]) # DL T2 & Row T1
println(Lifting_Programmes["DCG_CAP4"], [3,11,18,26]) # Squat T1 & OHP T2
println(Lifting_Programmes["DCG_CAP4"], [4,12,19,27]) # Bench T12 & Bench T22
println(Lifting_Programmes["DCG_CAP4"], [6,13,21,28]) # DL T1 & Row T2
println(Lifting_Programmes["DCG_CAP4"], [7,14,22,30]) # OHP T1 & Squat T2

max_stars_repo_name: dcelisgarza/Lifting