text:
module InplaceLinalg

using LinearAlgebra

export @inplace, InplaceException

include("declare_types.jl")
include("error_handling.jl")

include("triangle_update.jl")
include("loopy_triangle_updates.jl")

macro inplace(x)
    try
        esc(inplace(x))
    catch e
        :(ip_error($(e.msg)))
    end
end

function inplace(expr::Expr)
    #dump(expr)
    ## parse assignment or err
    lhs, ass, rhs = assignment(expr)
    if ass == :(*=)
        return :(InplaceLinalg.mult_update!($lhs, $rhs, Val(1)))
    elseif ass == :(/=)
        ## accept any RHS
        return :(InplaceLinalg.div_update!($lhs, /, $rhs))
    end
    ## select update expressions involving divisions
    if ass == :(=) && isa(rhs, Expr) && rhs.head == :call && length(rhs.args) == 3
        if rhs.args[1] == :/
            num, den = rhs.args[2:3]
            args = divupdate(lhs, num)
            return :(InplaceLinalg.div_update!($(args...), /, $den))
        elseif rhs.args[1] == :\
            den, num = rhs.args[2:3]
            args = divupdate(lhs, num)
            return :(InplaceLinalg.div_update!($(args...), \, $den))
        end
        ## fall through for other expressions with :=
    end
    term1, pm, term2 = terms(rhs)
    if ass == :(=) && term1 == 0
        pm == :- && ip_error("negated multiplicative update not available.")
        n, facs = multupdate(lhs, rhs)
        n > 0 && return :(InplaceLinalg.mult_update!($lhs, $(facs...), Val($n)))
    end
    if term1 != 0
        β, C = factors(term1, 2)
        #@assert lhs == C "First term must be linear in the LHS" lhs C
        lhs == C || return :(ip_error("First term must be a multiple of the LHS"))
    else
        β = 0
    end
    α, A, B = factors(term2)
    if ass in [:(+=), :(-=)]
        β = dobeta(Val{ass}, β)
        if ass == :(-=)
            #α = negate(α)
            pm = pm == :+ ? :- : :+
        end
        ass = :(=)
    end
    ass in [:(/=), :(*=)] && return :(ip_error("Unexpected assignment operator"))
    ## println((β, α, A, B))
    isa(B, Symbol) ||
        isa(B, Expr) && B.head == Symbol("'") && isa(B.args[1], Symbol) ||
        return :(ip_error("Too complex expression for inplace"))
    return :(InplaceLinalg.C_AB!($lhs, $β, $pm, $α, $A, $B))
end
inplace(x) = x

function assignment(expr::Expr)
    expr.head in [:(=), :(+=), :(*=), :(/=), :(-=)] || error("Unknown assignment operator in " * string(expr))
    isa(expr.args[1], Symbol) || error("LHS should be a symbolic variable")
    return expr.args[1], expr.head, expr.args[2]
end

function divupdate(lhs::Symbol, num::Symbol)
    lhs == num || error("LHS must be equal to numerator in updating divide")
    return tuple(num)
end
function divupdate(lhs::Symbol, num::Expr)
    num.head == :call && length(num.args) == 3 && num.args[1] == :* || error("Numerator must be simple multiplicative expression")
    matches = lhs .== num.args[2:3]
    sum(matches) == 1 || error("LHS must appear exactly once in multiplicative numerator expression")
    return tuple(num.args[2 .+ .!matches]...) ## return args in correct order
end

function terms(expr::Expr)
    if expr.head == :call && length(expr.args) == 3 && expr.args[1] in (:+,:-)
        return expr.args[[2,1,3]]
    else
        return 0, :+, expr
    end
end
terms(x) = 0, :+, x

function multupdate(lhs::Symbol, rhs::Expr)
    if rhs.head == :call && length(rhs.args) ≤ 4 && rhs.args[1] == :*
        factors = checksigns(flatten_mult(rhs.args[2:end]))
        matches = lhs .== factors
        sum(matches) == 1 && return findfirst(matches), factors[.!matches]
    end
    return 0, rhs
end
multupdate(lhs::Symbol, rhs) = 0, rhs

function checksigns(a::Vector)
    n = 0
    ret = []
    for f in a
        if isa(f, Expr) && f.head == :call && length(f.args) == 2 && f.args[1] == :-
            push!(ret, f.args[2])
            n += 1
        else
            push!(ret, f)
        end
    end
    isodd(n) && pushfirst!(ret, -1)
    return ret
end

function flatten_mult(factors::Vector)
    ret = []
    for factor in factors
        expanded = flatten_mult(factor)
        if isa(expanded, Array)
            for fac in expanded
                push!(ret, fac)
            end
        else
            push!(ret, expanded)
        end
    end
    return ret
end
function flatten_mult(expr::Expr)
    if expr.head == :call && length(expr.args) > 2 && expr.args[1] == :*
        return flatten_mult(expr.args[2:end])
    else
        return expr
    end
end
flatten_mult(x) = x

function factors(expr::Expr, n=3)
    ret = []
    scalars = []
    if expr.head == :call && expr.args[1] == :*
        for i in 2:min(length(expr.args), n)
            factor = pop!(expr.args)
            scalar, factor = factors(factor, 2)
            pushfirst!(ret, factor)
            scalar != 1 && pushfirst!(scalars, scalar)
        end
        while length(expr.args) > 1
            scalar, = factors(pop!(expr.args), 1)
            scalar != 1 && pushfirst!(scalars, scalar)
        end
        scalar = length(scalars) == 0 ? 1 : length(scalars) == 1 ? scalars[1] : Expr(:call, :*, scalars...)
        return tuple(scalar, ret...)
    else
        return fill(1, n-1)..., expr
    end
end
factors(x, n=3) = fill(1, n-1)..., x

function factor(expr::Expr)
    nfactors = length(expr.args) - 1
    nfactors == 0 && return 1
    nfactors == 1 && return expr.args[2]
    return expr
end

dobeta(::Type{Val{:(+=)}}, x::Number) = x + 1
dobeta(::Type{Val{:(-=)}}, x::Number) = (1 - x)
dobeta(::Type{Val{:(+=)}}, x) = :($x + 1)
dobeta(::Type{Val{:(-=)}}, x) = :(1 - $x)

negate(x::Number) = -x
negate(x) = :(-$x)

#C_AB!(C, β, α, A, B) = ip_error(": inplace assignment for this combination of types not implemented.")
#C_div!(C, α, B, div, A) = ip_error(": inplace assignment for this combination of types not implemented.")

include("extend_ldiv_and_rdiv.jl")
include("extend_lmul_and_rmul.jl")

#include("C_AB.jl")

# isnumber(x) = false
# isnumber(x::Number) = true
# function reduce_factors(f...)
#     ff = collect(f)
#     num = filter(isnumber,ff)
#     isempty(num) && return f
#     p = prod(num)
#     not_num = filter(!isnumber,ff)
#     p==1 && return tuple(not_num...)
#     return tuple(p, not_num...)
# end

reduce_factors(A, B, C) = (A,B,C)
reduce_factors(A::Number, B, C) = A==1 ? (B,C) : (A,B,C)
reduce_factors(A, B::Number, C) = B==1 ? (A,C) : (B,A,C)
reduce_factors(A, B, C::Number) = C==1 ? (A,B) : (C,A,B)
reduce_factors(A::Number, B::Number, C::Number) = (A*B*C,)
reduce_factors(A, B::Number, C::Number) = (P=B*C; P==1 ? (A,) : (P,A))
reduce_factors(A::Number, B, C::Number) = (P=A*C; P==1 ? (B,) : (P,B))
reduce_factors(A::Number, B::Number, C) = (P=A*B; P==1 ? (C,) : (P,C))

C_AB!(C, β, pm, α, A, B) = add_update!(C, β, pm, reduce_factors(α, A, B)...)
#add_update!(C, β::Number, pm::Function) = add_update(C, β, pm, 1)

include("div_update.jl")
include("mult_update.jl")
include("add_update.jl")

end # module


text:
#
#   File:      feasrepairex1.jl
#
#   Purpose:    To demonstrate how to use the MSK_relaxprimal function to
#               locate the cause of an infeasibility.
#
#   Syntax: On command line
#           julia-basic feasrepairex1.jl feasrepair.lp

using Mosek

printstream(msg::String) = print(msg)

maketask() do task
    putstreamfunc(task,MSK_STREAM_LOG,printstream)

    inputfile = ARGS[1]

    # Read data
    readdata(task,inputfile)

    nvars = getnumvar(task)
    ncons = getnumcon(task)

    putintparam(task,MSK_IPAR_LOG_FEAS_REPAIR,3)

    primalrepair(task,ones(ncons),ones(ncons),ones(nvars),ones(nvars))

    sum_viol = getdouinf(task,MSK_DINF_PRIMAL_REPAIR_PENALTY_OBJ)
    @printf("Minimized sum of violations = %e\n", sum_viol)

    optimize(task)

    solutionsummary(task,MSK_STREAM_MSG)
end

ext: jl


text:
abstract type AlgebraicTrait end
struct CrispTrait <: AlgebraicTrait end
struct FuzzyTrait <: AlgebraicTrait end

"""
Sono in dubbio..
Ho un'idea: Mettere tipo ExistentialIntervalRelation{R,A} dove R è il tipo di relazione (:L, :A, ecc.)
e A è il tipo di algebra (:Boolean, :HeytingChain), però questo non fa emergere l'algebra.. uff...

ho risolto facendo
struct ExistentialIntervalRelation{T}   <: AbstractExistentialIntervalRelation
    algebra::AbstractAlgebra
end
"""

"""
TODO devo cambiare questo? cioè, devo fare i trait crips/fuzzy sulle relazioni/proposizioni? adesso è sulle algebre
"""
# Default behavior is crisp.
AlgebraicTrait(::Type) = CrispTrait()
AlgebraicTrait(::Type{<:BooleanAlgebra}) = CrispTrait()
AlgebraicTrait(::Type{<:HeytingAlgebra}) = FuzzyTrait()

iscrisp(x::T) where {T} = iscrisp(AlgebraicTrait(T), x)
iscrisp(::CrispTrait, x) = true
iscrisp(::FuzzyTrait, x) = false
isfuzzy(x::T) where {T} = !iscrisp(x)

abstract type OntologicalTrait end
struct PropositionalTrait   <: OntologicalTrait end
struct PointTrait           <: OntologicalTrait end
struct IntervalTrait        <: OntologicalTrait end
struct RectangleTrait       <: OntologicalTrait end
struct CubeTrait            <: OntologicalTrait end

# Default behavior is propositional.
OntologicalTrait(::Type) = PropositionalTrait()
# Propositional trait.
# OntologicalTrait(::Type{<:PropositionalState})              = PropositionalTrait()
OntologicalTrait(::Type{<:AbstractPropositionalRelation})   = PropositionalTrait()
# Point trait.
# OntologicalTrait(::Type{<:Point})                           = PointTrait()
OntologicalTrait(::Type{<:AbstractPointRelation})           = PointTrait()
# Interval trait.
OntologicalTrait(::Type{<:Interval})                        = IntervalTrait()
OntologicalTrait(::Type{<:AbstractIntervalRelation})        = IntervalTrait()
# Rectangle trait.
# OntologicalTrait(::Type{<:Rectangle})                       = RectangleTrait()
OntologicalTrait(::Type{<:AbstractRectangleRelation})       = RectangleTrait()
# Cube trait.
# OntologicalTrait(::Type{<:Cube})                            = CubeTrait()
OntologicalTrait(::Type{<:AbstractCubeRelation})            = CubeTrait()

# Is modal relation only if it is not propositional.
ismodal_relation(x::T) where {T}            = ismodal_relation(OntologicalTrait(T), x)
ismodal_relation(::PropositionalTrait, x)   = false
ismodal_relation(::PointTrait, x)           = true
ismodal_relation(::IntervalTrait, x)        = true
ismodal_relation(::RectangleTrait, x)       = true
ismodal_relation(::CubeTrait, x)            = true

# # check_relation(r::T, os::NTuple{2, AbstractOntology}) where {T} = check_relation(OntologicalTrait(T), r, os)
# # function check_relation(::IntervalTrait, r, os)
# #     @show r
# #     @show os
# # end
#
# # function check_relation(r::T, o1::S1, o2::S2) where {T,S1,S2}
# #     check_relation(OntologicalTrait(T), OntologicalTrait(S1), OntologicalTrait(S2), r, o1, o2)
# # end
# # function check_relation(::IntervalTrait, ::IntervalTrait, ::IntervalTrait, r, o1, o2)
# #     @show r
# #     @show o1
# #     @show o2
# #     nothing
# # end
#
# function check_relation(r::RT, os::NTuple{2,OT}) where {RT,OT}
#     check_relation(OntologicalTrait(RT), OntologicalTrait(OT), r, os)
# end
#
# function check_relation(::IntervalTrait, ::IntervalTrait, r::AbstractIntervalRelation, os::NTuple{2,Interval})
#     if iscrisp(r.algebra)
#         check_crisp_relation(IntervalTrait(), IntervalTrait(), r, os)
#     else
#         print("fuzzy relation.. boss..")
#     #     check_fuzzy_relation(::IntervalTrait, ::IntervalTrait, r, os)
#     end
# end
#
# function check_crisp_relation(::IntervalTrait, ::IntervalTrait, r::AbstractIntervalRelation, os::NTuple{2,Interval})
#     check_crisp_relation(r, os)
# end
#
# # TODO qua uso <, in realtà dovrebbe essere dinamico, cioè usare il crisp o fuzzy <
#
# function check_crisp_relation(::ExistentialIntervalRelation{:L}, os::NTuple{2,Interval})
#     return os[1].y < os[2].x ? 0.0 : 1.0
# end
#
# # check_relation(r::T, os::NTuple{2, AbstractOntology}) where {T} = check_relation(OntologicalTrait(T), r, os)
#
# # function check_relation(::IntervalTrait, r, os)
# #     @show r
# #     @show os
# # end


text:
"""
update_state!(model::AbstractModel)

Update the model to the current situation
"""
function update_state!(model)
    update_surface_elevation!(model)
    update_geometry_on_uv_grids!(model)
    update_height_above_floatation!(model)
    update_grounded_fraction_on_huv_grids!(model)
    update_accumulation_rate!(model)
    update_basal_melt!(model)
    update_weertman_c!(model)
    update_dsdh!(model)
    update_model_velocities!(model)
    update_velocities_on_h_grid!(model)
    update_dhdt!(model)
    update_model_wavelets!(model)
    return nothing
end

"""
    update_surface_elevation!(model::AbstractModel)

Adjust surface elevation to hydrostatic equilibrium.
"""
function update_surface_elevation!(model::AbstractModel)
    @unpack params=model
    @unpack gh=model.fields
    gh.s[gh.mask] .= max.(gh.b[gh.mask]+gh.h[gh.mask],
                          params.sea_level_wrt_geoid .+ gh.h[gh.mask]*(1-params.density_ice./params.density_ocean))
    return model
end

"""
    update_geometry_on_uv_grids!(model::AbstractModel)

Interpolate thickness and surface elvation from h-grid to u- and v-grids.

"""
function update_geometry_on_uv_grids!(model::AbstractModel)
    @unpack gh,gu,gv,gc=model.fields
    onesvec=ones(gh.nxh*gh.nyh)
    gu.h[gu.mask].=(gu.samp*(gu.cent'*(gh.crop*gh.h[:])))./(gu.samp*(gu.cent'*(gh.crop*onesvec)))
    gu.s[gu.mask].=(gu.samp*(gu.cent'*(gh.crop*gh.s[:])))./(gu.samp*(gu.cent'*(gh.crop*onesvec)))
    gv.h[gv.mask].=(gv.samp*(gv.cent'*(gh.crop*gh.h[:])))./(gv.samp*(gv.cent'*(gh.crop*onesvec)))
    gv.s[gv.mask].=(gv.samp*(gv.cent'*(gh.crop*gh.s[:])))./(gv.samp*(gv.cent'*(gh.crop*onesvec)))
    return model
end

"""
    update_height_above_floatation!(model::AbstractModel)

Update height above floatation. Zero value is used to define location of grounding line.
"""
function update_height_above_floatation!(model::AbstractModel)
    @unpack params=model
    @unpack gh=model.fields
    gh.haf .= height_above_floatation.(gh.h,gh.b,Ref(params))
    return model
end

"""
    update_grounded_fraction_on_huv_grids!(model::AbstractModel)

Update grounded area fraction on h-, u-, and v-grids for use in subgrid parameterisation.
"""
function update_grounded_fraction_on_huv_grids!(model::AbstractModel)
    @unpack gh,gu,gv = model.fields
    (gfh,gfu,gfv)=pos_fraction(gh.haf;mask=gh.mask)
    gh.grounded_fraction[:] .= gfh[:]
    gu.grounded_fraction[:] .= gfu[:]
    gv.grounded_fraction[:] .= gfv[:]
    return model
end

"""
    update_accumulation_rate!(model::AbstractModel)

Update the accumulation rate.
"""
function update_accumulation_rate!(model::AbstractModel)
    @unpack params = model
    @unpack gh=model.fields
    gh.accumulation .= params.accumulation_rate
    return model
end

"""
    update_basal_melt!(model::AbstractModel)

Update the basal melt rate.
"""
function update_basal_melt!(model::AbstractModel)
    update_melt_rate!(model.melt_rate, model.fields, model.grid)
    return model
end

"""
    update_weertman_c!(model::AbstractModel)

Update coefficient used in the sliding law to account for migration of grounding line.
"""
function update_weertman_c!(model::AbstractModel)
    @unpack gh=model.fields
    @unpack params=model
    gh.weertman_c .= params.weertman_c .* gh.grounded_fraction
    return model
end

"""
    update_dsdh!(model::AbstractModel)

Compute change of surface elevation per unit thickness change, accounting for hydrostatic adjustment.
"""
function update_dsdh!(model::AbstractModel)
    @unpack gh,gu,gv=model.fields
    @unpack params = model
    gh.dsdh .= (1.0 - params.density_ice./params.density_ocean) .+
           (params.density_ice./params.density_ocean).*gh.grounded_fraction;
    return model
end

"""
    update_model_velocity!(model::AbstractModel)

Wrapper function for that which updates the model velocities on the u, v grids (update_velocities in separate file)
"""
function update_model_velocities!(model::AbstractModel)
    update_velocities!(model)
    return model
end

"""
    update_velocities_on_h_grid!(model::AbstractModel)

Update the velocities (depth averaged, surface and bed) on the h grid
"""
function update_velocities_on_h_grid!(model)
    @unpack gh,gu,gv = model.fields
    #depth averaged velocities
    gh.u[:] .= gu.cent*gu.u[:] #(gu.u[1:end-1,:] + gu.u[2:end,:])./2
    gh.v[:] .= gv.cent*gv.v[:] #(gv.v[:,1:end-1] + gv.v[:, 2:end])./2

    #bed velocities
    gh.ub .= gh.u ./ (1 .+ (gh.β .* gh.quad_f2))
    gh.vb .= gh.v ./ (1 .+ (gh.β .* gh.quad_f2))

    #surface velocities
    gh.us .= gh.ub .* (1 .+ (gh.β .* gh.quad_f1))
    gh.vs .= gh.vb .* (1 .+ (gh.β .* gh.quad_f1))
    return model
end
"""
    update_dhdt!(model::AbstractModel)

Evaluate rate of change of thickness using mass conservation.
"""
function update_dhdt!(model::AbstractModel)
    @unpack gh,gu,gv=model.fields
    gh.dhdt[gh.mask].=gh.samp*(gh.accumulation[:] .- gh.basal_melt[:] .-
             (  (gu.∂x*(gu.crop*(gu.h[:].*gu.u[:]))) .+ (gv.∂y*(gv.crop*(gv.h[:].*gv.v[:]))) ) )
    return model
end

"""
    update_model_wavelets(model::AbstractModel)

Wrapper function for that which updates the model wavelets
"""
function update_model_wavelets!(model::AbstractModel)
    update_wavelets!(model)
    return model
end


text:
@assert begin
    bf=get_basefun(
        :(a+x*b+y*c),
        [
            :x, :y
        ],
        [
            :a,
            :b,
            :c
        ]
    )

    isapprox(
        bf(1.0, 2.0, -1.0, 1.0, 2.0), 4.0
    )
end

@assert begin
    bfd=get_basefun_derivative(
        :(a+x*b+y*c),
        [
            :x, :y
        ],
        [
            :a,
            :b,
            :c
        ],
        [1, 0]
    )

    isapprox(
        bfd(
            1.0, -2.0, 4.0, 1.0, -3.0
        ),
        1.0
    )
end

@assert begin
    blpc=get_linpar_contribution(
        :(a+b*x+c*y),
        [
            :x,
            :y
        ],
        [
            :a,
            :b,
            :c
        ],
        :c,
        [0, 0]
    )

    isapprox(
        blpc(
            1.0, -2.0, 4.0, 1.0, -3.0
        ),
        -2.0
    )
end


text:
using Distributed
addprocs(2)

@everywhere begin
using MLJ
using MLJBase
using Test
using Random
end

@testset "ensembles" begin
    @test include("ensembles.jl")
end

@testset "matching models to data" begin
    @test include("model_matching.jl")
end

@testset "scitypes" begin
    @test include("scitypes.jl")
end


text:
# yuck. there were less fucks given making this
"""
    FVMesh{T} <: Mesh{T}

A finite volume discretisation scheme for a DiscretisedFluidQueue.

# Arguments
 - `nodes`: the cell edges (the nodes are at the center of the cell)
 - `order`: the order of the polynomial interpolation used to approximate the flux at the
 cell edges.
"""
struct FVMesh{T} <: Mesh{T}
    nodes::T
    order::Int
    function FVMesh{T}(nodes::T,order::Int) where T
        (order<=0)&&throw(DomainError("order must be positive"))
        (nodes[1]!=0.0)&&throw(DomainError("first node must be 0"))
        return new{T}(nodes,order)
    end
end
# Convenience constructors
FVMesh(nodes::T,order::Int) where T = FVMesh{T}(nodes,order)

"""

    n_bases_per_cell(mesh::FVMesh)

Constant=1
"""
n_bases_per_cell(mesh::FVMesh) = 1
_order(mesh::FVMesh) = mesh.order

"""

    cell_nodes(mesh::FVMesh)

The cell centres
"""
cell_nodes(mesh::FVMesh) = Array(((mesh.nodes[1:end-1] + mesh.nodes[2:end]) / 2 )')

n_bases_per_phase(mesh::FVMesh) = n_intervals(mesh)

function MakeFVFlux(mesh::Mesh)
    order = _order(mesh)
    nNodes = n_bases_per_phase(mesh)
    F = zeros(Float64,nNodes,nNodes)
    ptsLHS = Int(ceil(order/2))
    interiorCoeffs = lagrange_polynomials(cell_nodes(mesh)[1:order],mesh.nodes[ptsLHS+1])
    for n in 2:nNodes
        evalPt = mesh.nodes[n]
        if n-ptsLHS-1 < 0
            nodesIdx = 1:order
            nodes = cell_nodes(mesh)[nodesIdx]
            coeffs = lagrange_polynomials(nodes,evalPt)
        elseif n-ptsLHS-1+order > nNodes
            nodesIdx = (nNodes-order+1):nNodes
            nodes = cell_nodes(mesh)[nodesIdx]
            coeffs = lagrange_polynomials(nodes,evalPt)
        else
            nodesIdx =  (n-ptsLHS-1) .+ (1:order)
            coeffs = interiorCoeffs
        end
        F[nodesIdx,n-1:n] += [-coeffs coeffs]./Δ(mesh)[n-1]
    end
    F[end-order+1:end,end] += -lagrange_polynomials(cell_nodes(mesh)[end-order+1:end],mesh.nodes[end])./Δ(mesh)[end]
    return F
end

function build_full_generator(dq::DiscretisedFluidQueue{FVMesh{T}}; v::Bool=false) where T
    model = dq.model

    order = _order(dq.mesh)
    F = MakeFVFlux(dq.mesh)

    C = rates(model.S)
    n₋ = N₋(model.S)
    n₊ = N₊(model.S)

    B = SparseArrays.spzeros(
        Float64,
        n_phases(model) * n_intervals(dq) + n₋ + n₊,
        n_phases(model) * n_intervals(dq) + n₋ + n₊,
    )
    B[n₋+1:end-n₊,n₋+1:end-n₊] = SparseArrays.kron(
            model.T,
            SparseArrays.I(n_intervals(dq))
        )

    ## Make QBD index
    QBDidx = qbd_idx(dq)

    # Boundary conditions
    T₋₋ = model.T[_has_left_boundary.(model.S),_has_left_boundary.(model.S)]
    T₊₋ = model.T[_has_right_boundary.(model.S),:].*((C.<0)')
    T₋₊ = model.T[_has_left_boundary.(model.S),:].*((C.>0)')
    T₊₊ = model.T[_has_right_boundary.(model.S),_has_right_boundary.(model.S)]
    # yuck
    begin
        nodes = cell_nodes(dq)[1:order]
        coeffs = lagrange_polynomials(nodes,dq.mesh.nodes[1])
        idxdown = ((1:order).+n_bases_per_phase(dq)*(findall(_has_left_boundary.(model.S)) .- 1)')[:] .+ n₋
        down_rates = LinearAlgebra.diagm(0 => C[_has_left_boundary.(model.S)])
        B[idxdown, 1:n₋] = LinearAlgebra.kron(down_rates,-coeffs)
    end
    outLower = [
        T₋₊./Δ(dq,1) SparseArrays.zeros(n₋,n₊+(n_intervals(dq)-1)*n_phases(model))
    ]
    begin
        nodes = cell_nodes(dq)[end-order+1:end]
        coeffs = lagrange_polynomials(nodes,dq.mesh.nodes[end])
        idxup =
            ((1:order).+n_bases_per_phase(dq)*(findall(_has_right_boundary.(model.S)) .- 1)')[:] .+
            (n₋ + n_bases_per_phase(dq) - order)
        B[idxup, (end-n₊+1):end] = LinearAlgebra.kron(
            LinearAlgebra.diagm(0 => C[_has_right_boundary.(model.S)]),
            coeffs,
        )
    end
    outUpper = [
        SparseArrays.zeros(n₊,n₋+(n_intervals(dq)-1)*n_phases(model)) T₊₋./Δ(dq,n_intervals(dq))
    ]

    B[1:n₋,QBDidx] = [T₋₋ outLower]
    B[end-n₊+1:end,QBDidx] = [outUpper T₊₊]
    for i = 1:n_phases(model)
        idx = ((i-1)*n_intervals(dq)+1:i*n_intervals(dq)) .+ n₋
        if C[i] > 0
            B[idx, idx] += C[i] * F
        elseif C[i] < 0
            B[idx, idx] += abs(C[i]) * F[end:-1:1,end:-1:1]
        end
    end

    out = FullGenerator(B[QBDidx,QBDidx],dq)#, mesh.Fil)
    v && println("FullGenerator created with keys ", keys(out))
    return out
end


text:
function get(dr::DataReaderGoogleQuotes, symb::DataSymbol)
    url = "http://www.google.com/finance/info"
    error("ToDo")
end


text:
x = 3
println(2x^2 - 3x + 1)
# 10
println(1.5x^2 - .5x + 1)
# 13.0

println(zero(Float32))
# 0.0f0
println(zero(Int32))
# 0
println(one(Int32))
# 1
println(one(BigFloat))
# 1.000000000000000000000000000000000000000000000000000000000000000000000000000000

max_stars_repo_name: dannypsnl/languages-learn

max_issues_repo_name: dannypsnl/languages-learn

max_forks_repo_name: dannypsnl/languages-learn


text:
function bisect_left(A, p, r, v)
	i = p
	if p < r
	   q = floor(Int, (p + r) / 2)
	   if v <= A[q]
		   i = bisect_left(A, p, q, v)
	   else
		   i = bisect_left(A, q + 1, r, v)
	   end
	end
	return i
end

function bisect_right(A, p, r, v)
	i = p
	if p < r
	   q = floor(Int, (p + r) / 2)
	   if v < A[q]
		   i = bisect_right(A, p, q, v)
	   else
		   i = bisect_right(A, q + 1, r, v)
	   end
	end
	return i
end

function bisect_right(A, v)
	return bisect_right(A, 1, length(A) + 1, v)
end

function bisect_left(A, v)
	return bisect_left(A, 1, length(A) + 1, v)
end


text:
"""
    ChebyshevPoints{T = Float64}(n)

The `n` Chebyshev points on [-1,1] in type `T`.
"""
struct ChebyshevPoints{T <: AbstractFloat} <: AbstractVector{T}
    data::Vector{T}
    function ChebyshevPoints{T}(n::Integer) where {T}
        if n == 1
            return new(T[0])
        else
            data = Vector{T}(undef,n)
            for i = 1:n
                data[i] = cos(T(π)*(n-i)/(n-1))
            end
            return new(data)
        end
    end
end
ChebyshevPoints(n) = ChebyshevPoints{Float64}(n)
Base.size(x::ChebyshevPoints) = size(x.data)
Base.getindex(x::ChebyshevPoints, i::Int) = x.data[i]

function prodpot(x::ChebyshevPoints)
    T = eltype(x); n = length(x)
    n == 1 && return [LogNumber(1,float(T(0)))]
    p = LogNumber.( (-1).^(0:n-1), log(T(n-1)) - (n-2)*log(T(2)) )
    p[1] *= LogNumber(true,log(2)); p[end] *= LogNumber(true,log(2))
    return p
end


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Portfolio Optimization

# **Originally Contributed by**: Arpit Bhatia

# Optimization models play an increasingly important role in financial
# decisions. Many computational finance problems can be solved efficiently using
# modern optimization techniques.
#
# This tutorial solves the famous Markowitz Portfolio Optimization problem with
# data from [lecture notes from a course taught at Georgia Tech by Shabir Ahmed](https://www2.isye.gatech.edu/~sahmed/isye6669/).

# This tutorial uses the following packages
using JuMP
import Ipopt
import Statistics

# Suppose we are considering investing 1000 dollars in three non-dividend paying
# stocks, IBM (IBM), Walmart (WMT), and Southern Electric (SEHI), for a
# one-month period.

# This means we will use the money to buy shares of the three stocks at the
# current market prices, hold these for one month, and sell the shares off at
# the prevailing market prices at the end of the month.
#
# As a rational investor, we hope to make some profit out of this endeavor,
# i.e., the return on our investment should be positive.
#
# Suppose we bought a stock at $p$ dollars per share in the beginning of the
# month, and sold it off at $s$ dollars per share at the end of the month. Then
# the one-month return on a share of the stock is $ \frac{s-p}{p} $.

# Since the stock prices are quite uncertain, so is the end-of-month return on
# our investment. Our goal is to invest in such a way that the expected
# end-of-month return is at least \$50 or 5%. Furthermore, we want to make
# sure that the “risk” of not achieving our desired return is minimum.

# Note that we are solving the problem under the following assumptions:
# 1. We can trade any continuum of shares.
# 2. No short-selling is allowed.
# 3. There are no transaction costs.

# We model this problem by taking decision variables $x_{i}, i=1,2,3,$ denoting
# the dollars invested in each of the 3 stocks.

# Let us denote by $\tilde{r}_{i}$ the random variable corresponding to the
# monthly return (increase in the stock price) per dollar for stock $i$.

# Then, the return (or profit) on $x_{i}$ dollars invested in stock $i$ is
# $\tilde{r}_{i} x_{i},$ and the total (random) return on our investment is
# $\sum_{i=1}^{3} \tilde{r}_{i} x_{i}.$ The expected return on our investment is
# then $\mathbb{E}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right]=\sum_{i=1}^{3} \overline{r}_{i} x_{i},$
# where $\overline{r}_{i}$ is the expected value of the $\tilde{r}_{i}.$

# Now we need to quantify the notion of “risk” in our investment.

# Markowitz, in his Nobel prize winning work, showed that a rational investor’s
# notion of minimizing risk can be closely approximated by minimizing the
# variance of the return of the investment portfolio. This variance is given by:

# ```math
# \operatorname{Var}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right] = \sum_{i=1}^{3} \sum_{j=1}^{3} x_{i} x_{j} \sigma_{i j}
# ```

# where $\sigma_{i j}$ is the covariance of the return of stock $i$ with stock $j$.

# Note that the right hand side of the equation is the most reduced form of the
# expression and we have not shown the intermediate steps involved in getting to
# this form. We can also write this equation as:

# ```math
# \operatorname{Var}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right] =x^{T} Q x
# ```

# Where $Q$ is the covariance matrix for the random vector $\tilde{r}$.

# Finally, we can write the model as:

# ```math
# \begin{aligned}
# \min x^{T} Q x \\
# \text { s.t. } \sum_{i=1}^{3} x_{i} \leq 1000.00 \\
# \overline{r}^{T} x \geq 50.00 \\
# x \geq 0
# \end{aligned}
# ```

# After that long discussion, lets now use JuMP to solve the portfolio
# optimization problem for the data given below.

# | Month        |  IBM     |  WMT    |  SEHI  |
# |--------------|----------|---------|--------|
# | November-00  |  93.043  |  51.826 |  1.063 |
# | December-00  |  84.585  |  52.823 |  0.938 |
# | January-01   |  111.453 |  56.477 |  1.000 |
# | February-01  |  99.525  |  49.805 |  0.938 |
# | March-01     |  95.819  |  50.287 |  1.438 |
# | April-01     |  114.708 |  51.521 |  1.700 |
# | May-01       |  111.515 |  51.531 |  2.540 |
# | June-01      |  113.211 |  48.664 |  2.390 |
# | July-01      |  104.942 |  55.744 |  3.120 |
# | August-01    |  99.827  |  47.916 |  2.980 |
# | September-01 |  91.607  |  49.438 |  1.900 |
# | October-01   |  107.937 |  51.336 |  1.750 |
# | November-01  |  115.590 |  55.081 |  1.800 |

stock_data = [
    93.043 51.826 1.063;
    84.585 52.823 0.938;
    111.453 56.477 1.000;
    99.525 49.805 0.938;
    95.819 50.287 1.438;
    114.708 51.521 1.700;
    111.515 51.531 2.540;
    113.211 48.664 2.390;
    104.942 55.744 3.120;
    99.827 47.916 2.980;
    91.607 49.438 1.900;
    107.937 51.336 1.750;
    115.590 55.081 1.800;
]

# Calculating stock returns

stock_returns = Array{Float64}(undef, 12, 3)
for i in 1:12
    stock_returns[i, :] = (stock_data[i + 1, :] .- stock_data[i, :]) ./ stock_data[i, :]
end
stock_returns

# Calculating the expected value of monthly return:

r = Statistics.mean(stock_returns, dims = 1)

# Calculating the covariance matrix Q

Q = Statistics.cov(stock_returns)

# JuMP Model

portfolio = Model(Ipopt.Optimizer)
set_silent(portfolio)
@variable(portfolio, x[1:3] >= 0)
@objective(portfolio, Min, x' * Q * x)
@constraint(portfolio, sum(x) <= 1000)
@constraint(portfolio, sum(r[i] * x[i] for i = 1:3) >= 50)
optimize!(portfolio)

objective_value(portfolio)

#-

value.(x)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Chordal
using Test

using LinearAlgebra, SparseArrays
using QDLDL
const CD = Chordal

include("chordal_graph.jl")
include("clique_tree.jl")
include("clique_graph.jl")
include("completion.jl")
include("edm.jl")
include("utils.jl")


text:
# Lax-Wendroff Scheme

struct LaxWendroffScheme <: AbstractFVAlgorithm end

function update_flux_value(uold,node_idx,dt,dx,mesh,Flux,alg::LaxWendroffScheme)
    # Local speeds of propagation
    @inbounds ul=cellval_at_left(node_idx,uold,mesh)
    @inbounds ur=cellval_at_right(node_idx,uold,mesh)
    aa = fluxρ(0.5*(ul+ur),Flux)
    # Numerical Fluxes
    return 0.5*(Flux(ul)+Flux(ur))-aa*dt/(2*dx)*(Flux(ur)-Flux(ul))
end
# Ritchmeyer Two-step Lax-Wendroff Method
# Reference:
# R. Leveque. Finite Volume Methods for Hyperbolic Problems.Cambridge University
# Press. New York 2002

struct LaxWendroff2sScheme <: AbstractFVAlgorithm end

function update_flux_value(uold,node_idx,dt,dx,mesh,Flux,alg::LaxWendroff2sScheme)
    # Local speeds of propagation
    @inbounds ul=cellval_at_left(node_idx,uold,mesh)
    @inbounds ur=cellval_at_right(node_idx,uold,mesh)
    # Numerical Fluxes
    return Flux(0.5*(ul+ur)-dx/(2*dt)*(Flux(ur)-Flux(ul)))
end


text:
#=##############################################################################
# DESCRIPTION
Testing models, dynamics, and control of a wind-harvesting aircraft (a.k.a.,
windcraft). This vehicle is a multirotor aircraft tethered to the ground in
similitud to the Makani M600 Project, whose rotors operate both as propellers
and turbines along sections of a circular path in the presence of cross wind.

This example is a snippet of the study published under J. Mehr, E. Alvarez, A.
Cardoza, and A. Ning, "Mixed-fidelity Aerodynamic Analysis of Wind Harvesting
Aircraft," AIAA AVIATION Forum, 2020. See the paper for further details on
geometry, rotor design, and analysis.

# AUTHORSHIP
  * Author    : Judd Mehr and Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Apr 2020
  * License   : MIT
=###############################################################################

# ------------ MODULES ---------------------------------------------------------
# Load simulation engine
# import FLOWUnsteady
reload("FLOWUnsteady")
uns = FLOWUnsteady
vlm = uns.vlm
gt = uns.gt

using PyPlot

# ------------ GLOBAL VARIABLES ------------------------------------------------
# Default path where to save data
extdrive_path = "/media/edoalvar/MyExtDrive/simulationdata7/"
# extdrive_path = "temps/"

# ------------ HEADERS ---------------------------------------------------------
for header_name in ["geometry", "kinematics"]
    include("windcraft_"*header_name*".jl")
end

# ------------ DRIVERS ---------------------------------------------------------

"""
    Visualize kinematic maneuver of the windcraft: saves it as vtk files, and
    calls Paraview visualizing the vehicle path.
"""
function visualize_maneuver_windcraft_kinematic(; save_path=extdrive_path*"windcraft_maneuver11/",
                                                    prompt=true,
                                                    run_name="windcraft",
                                                    verbose=true, v_lvl=0,
                                                    paraview=true)

    # Geometry parameters
    includewing     = true
    includetail     = false
    includepylons   = false
    includerotors   = false
    numrotors       = 8
    inlinerotors    = false
    includecontrols = false

    # Maneuver parameters
    R               = 135/2             # (m) radius of circle
    t_per_rev       = 3.0               # (s) time of one full revolution
    nrevs           = 0.75              # Revolutions to simulate
    RPMref          = 4.0*40.0*30/pi    # Reference RPM: 40m/s with tip speed ratio of 4

    Vmean           = 2*pi*R/(nrevs*t_per_rev) # (m/s) mean velocity along a full circle
    ttot            = nrevs*t_per_rev   # (s) total time to perform maneuver
    nsteps          = 60                # Time steps

    revinit         = 0.25              # Part of revolution where to start the simulation
    tinit           = revinit*t_per_rev # (s) initial time (sim time by end will be tinit+ttot)

    # Circular path parameters
    theta0          = 0.0
    thetan          = 360.0
    omegan          = [4.0/9.0, 20.0/27.0, 4.0/9.0]
    tn              = [0.0, 0.5, 1.0]

    # Generate maneuver
    gt.verbalize("MANEUVER GENERATION", v_lvl, verbose)

    maneuver = generate_maneuver_windcraft_kinematic(nrevs;
                                                       disp_plot        = false,
                                                       includetail      = includetail,
                                                       includewing      = includewing,
                                                       includecontrols  = includecontrols,
                                                       includerotors    = includerotors,
                                                       theta0           = theta0,
                                                       thetan           = thetan,
                                                       omegan           = omegan,
                                                       tn               = tn)

    # Plot maneuver path and controls
    uns.plot_maneuver(maneuver; tstages=[])

    # Generate geometry
    gt.verbalize("VEHICLE GENERATION", v_lvl, verbose)

    (system, vlm_system, rotors,
     tilting_systems, rotor_systems,
        wake_system) = generate_geometry_windcraft(;
                                                    includewing     = includewing,
                                                    includetail     = includetail,
                                                    includepylons   = includepylons,
                                                    includerotors   = includerotors,
                                                    numrotors       = numrotors,
                                                    inlinerotors    = inlinerotors,
                                                    includecontrols = includecontrols,
                                                    xfoil           = false,
                                                    data_path       = uns.def_data_path,
                                                    run_name        = run_name,
                                                    verbose         = verbose,
                                                    v_lvl           = v_lvl+1)

    vehicle = uns.VLMVehicle(   system;
                                tilting_systems = tilting_systems,
                                rotor_systems   = rotor_systems,
                                vlm_system      = vlm_system,
                                wake_system     = wake_system,
                             )

    # Simulation setup
    gt.verbalize("SIMULATION GENERATION", v_lvl, verbose)

    Vref = Vmean                            # (m/s) reference velocity
    Vinit = Vref*maneuver.Vvehicle(tinit/ttot)   # (m/s) initial vehicle velocity
                                            # (rad/s) initial vehicle angular velocity
    angle1 = maneuver.anglevehicle(tinit/ttot)
    angle2 = maneuver.anglevehicle(tinit/ttot + 1e-12)
    Winit = pi/180 * (angle2-angle1)/(ttot*1e-12)

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
                                            Vinit=Vinit, Winit=Winit, t=tinit)

    save_vtk_optsargs = [(:save_horseshoes, false)]

    # Visualize maneuver
    gt.verbalize("STEPPING THROUGH MANEUVER", v_lvl, verbose)

    strn = uns.visualize_kinematics(simulation, nsteps, save_path;
                                    run_name=run_name,
                                    save_vtk_optsargs=save_vtk_optsargs,
                                    prompt=prompt, verbose=verbose, v_lvl=v_lvl,
                                    paraview=false
                                    )

    # # Move landing pad to landing area
    # vlm.vtk.lintransform!(grounds[2], eye(3), Vcruise*telapsed*[-0.25, 0, -0.0025])
    #
    # # Save ground
    # for (i, ground) in enumerate(grounds)
    #     gt.save(ground, run_name*"_Ground$i"; path=save_path)
    #     strn *= run_name*"_Ground$i.vtk;"
    # end

    # Call paraview
    if paraview
        run(`paraview --data="$save_path/$strn"`)
    end

    return strn
end

"""
    Generates geometry of the windcraft, saves it as vtk files, and calls
    Paraview visualizing the VTK geometry.
"""
function visualize_geometry_windcraft(; save_path   = extdrive_path*"windcraft_geometry00/",
                                        run_name    = "windcraft",
                                        prompt      = true,
                                        verbose     = true,
                                        v_lvl       = 0,
                                        paraview    = true,
                                        xfoil       = false,
                                        optargs...)

    gt.verbalize("GEOMETRY VISUALIZATION", v_lvl, verbose)

    # Generate the Geometry in order to visualize
    (system, vlm_system, rotors,
     tilting_systems, rotor_systems,
        wake_system) = generate_geometry_windcraft(; xfoil=xfoil,
                                                  data_path=uns.def_data_path,
                                                  run_name=run_name,
                                                  v_lvl=v_lvl+1,
                                                  verbose=verbose, optargs...)

    # Set up dummy values of RPM and freestream
    Vinf(x,t) = [1,0,0]
    vlm.setVinf(system, Vinf)
    for rotor in rotors
        vlm.setRPM(rotor, 6000)
    end

    # Save visualization of geometry
    gt.create_path(save_path, prompt)
    strn = vlm.save(system, run_name; save_horseshoes=false, path=save_path)

    # Call Paraview
    if paraview
        run(`paraview --data="$save_path/$strn"`)
    end

end


text:
import IterativeSolvers

struct GMRESSolver{L,R}
    linear_operator::L
    maxiter::Int
    restart::Int
    tol::R
end

function GMRESSolver(op; maxiter=0, restart=0, tol=sqrt(eps(real(eltype(op)))))

    m, n = size(op)
    @assert m == n

    maxiter == 0 && (maxiter = div(n, 5))
    restart == 0 && (restart = n)

    GMRESSolver(op, maxiter, restart, tol)
end

operator(solver::GMRESSolver) = solver.linear_operator

function solve(solver::GMRESSolver, b)
    op = operator(solver)
    x, ch = IterativeSolvers.gmres(op, b, log=true,  maxiter=solver.maxiter,
        restart=solver.restart, reltol=solver.tol)
    return x, ch
end

function Base.:*(solver::GMRESSolver, b)
    x, ch = solve(solver, b)
    println("Number of iterations: ", ch.iters)
    ch.isconverged || error("Iterative solver did not converge.")
    return x
end

function gmres(eq::DiscreteEquation; maxiter=0, restart=0, tol=0)

    test_space_dict  = eq.test_space_dict
    trial_space_dict = eq.trial_space_dict

    lhs = eq.equation.lhs
    rhs = eq.equation.rhs

    b = assemble(rhs, test_space_dict)
    Z = assemble(lhs, test_space_dict, trial_space_dict)

    if tol == 0
        invZ = GMRESSolver(Z, maxiter=maxiter, restart=restart)
    else
        invZ = GMRESSolver(Z, maxiter=maxiter, restart=restart, tol=tol)
    end
    x = invZ * b
end


text:
#!/usr/bin/env julia
using XAM,Statistics
bamName = ARGS[1]

function histQuals(bamInName)
	count = 0
    allquals=Int[]
    bamIn = open(BAM.Reader,bamInName)
    for a in bamIn
		count+=1
        push!(allquals,BAM.mappingquality(a))
    end

	println("Mean quality: ",mean(allquals))
	println("Count: ",count)
end

histQuals(bamName)

# time averagequal.jl /Users/james/projects/hla/capture/bam38/PC-HLA-EG4.bam
# Reading BAM...done.
# Average quality: 43.76016208794757
# averagequal.jl /Users/james/projects/hla/capture/bam38/PC-HLA-EG4.bam
# 15.48s user 0.94s system 99% cpu 16.422 total

# time averagequal.jl data/PC-HLA-EG4.bam
# Reading BAM...done.
# Average quality: 56.078356406954285
# averagequal.jl data/PC-HLA-EG4.bam
# 1.85s user 0.42s system 108% cpu 2.105 total

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
function backward_rand_kernel(x, Gf, Gpred, Phi) where T
    xf, Pf = meancov(Gf)
    xpred, Ppred = meancov(Gpred)

    J = Pf*Phi'/Ppred
    xs = xf +  J*(x - (Phi*xf  + b))
    Ps = Pf - J*Ppred*J' # as Ppred = M.Phi*P*M.Phi' + M.Q
    rand(Gaussian(xs, Ps))
end

struct BackwardRand{T,S}
    iter::T # iterate backwards through
    M::S
end

function dyniterate(B::BackwardRand, (x, state))
    ϕ = dyniterate(B.iter, state)
    ϕ === nothing && return nothing
    (n, U), state = ϕ
    Gf, Gpred = U[1], U[2]
    Gs = backward_rand_kernel(x, Gf, Gpred, Phi(RTS.M.sys, (n, U)))
    n=>x, (x, state)
end


text:
export NwtdirectionLDLT, NwtdirectionLDLT_LS, NewtonLDLT

# """
# Computes the Newton direction using the LDLᵗ factorization.
# Use the factorization from LDLFactorizations (JuliaSmoothOptimizers)
# """
# function NwtdirectionLDLT(H, g; verbose::Bool=false)
#     LDLT = ldl(H)
#     d = LDLT \ -g
#     return d
# end

function NwtdirectionLDLT(H, g; verbose :: Bool = false)
   T = eltype(H)
   n, = size(g)
   global L = Matrix{T}(undef, n, n)
   global D = Matrix{T}(undef, n, n)
   global pp = Vector{Int}(undef, n)
   global ρ = nothing
   global ncomp = nothing

   try
       (L, D, pp, ρ, ncomp) = ldlt_symm(H, 'r')
   catch
 	println("*******   Problem in LDLt")
       #res = PDataLDLt()
       #res.OK = false
       #return res
       d = NaN * zeros(n)
       return d
   end

    # A[pp,pp] = P*A*P' =  L*D*L'

    if true in isnan.(D)
 	println("*******   Problem in D from LDLt: NaN")
        println(" cond (H) = $(cond(H))")
        #res = PDataLDLt()
        #res.OK = false
        #return res
        d = NaN * zeros(n)
        return d
    end

   X = eigen(D)
   Δ = X.values
   Q =  X.vectors
   ϵ2 = sqrt(eps(T))
   Γ = max.(abs.(Δ), ϵ2)

    # Ad = P' * L * Q * Δ * Q' * L' *P d = -g
    # replace Δ by Γ to ensure positive definiteness
   d̃ = L \ g[pp]
   d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
   d = -d̂[invperm(pp)]
   return d
end

function NwtdirectionLDLT_LS(H, g; verbose :: Bool = false)
   T = eltype(H)
   n, = size(g)
   global L = Matrix{T}(undef, n, n)
   global D = Matrix{T}(undef, n, n)
   global pp = Vector{Int}(undef, n)
   global ρ = nothing
   global ncomp = nothing
   global u = nothing

   try
       (L, D, pp, ρ, ncomp) = ldlt_symm(H, 'r')
   catch
 	println("*******   Problem in LDLt")
       #res = PDataLDLt()
       #res.OK = false
       #return res
       d = NaN * zeros(n)
       return d
   end

    # A[pp,pp] = P*A*P' =  L*D*L'

    if true in isnan.(D)
 	println("*******   Problem in D from LDLt: NaN")
        println(" cond (H) = $(cond(H))")
        #res = PDataLDLt()
        #res.OK = false
        #return res
        d = NaN * zeros(n)
        return d
    end

   X = eigen(D)
   Δ = X.values
   Q =  X.vectors
   ϵ2 = sqrt(eps(T))
   Γ = max.(abs.(Δ), ϵ2)
   @show Γ
   @show Γ .== ϵ2
   if (true in (Γ .== ϵ2))
      # too_diff = true
      u = 1
   else
      # too_diff = false
      u = 0
   end

    # Ad = P' * L * Q * Δ * Q' * L' *P d = -g
    # replace Δ by Γ to ensure positive definiteness
   d̃ = L \ g[pp]
   d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
   d = -d̂[invperm(pp)]
   return d, u
end

function NewtonLDLT(nlp          :: AbstractNLPModel,
                       nlp_stop     :: NLPStopping;
                       verboseLS    :: Bool = false,
                       kwargs...)
    return  Newton(nlp,
                   nlp_stop;
                   verboseLS = verboseLS,
                   Nwtdirection = NwtdirectionLDLT,
                   hessian_rep = hessian_dense,
                   kwargs...)
end


text:
module TapirSchedulersBenchmarks

using BenchmarkTools: Benchmark, BenchmarkGroup

include("bench_fib.jl")
include("bench_tarai.jl")
include("bench_full.jl")

function setup()
    suite = BenchmarkGroup()
    suite["fib"] = Benchfib.setup()
    return suite
end

function set_smoke_params!(bench)
    bench.params.seconds = 0.001
    bench.params.evals = 1
    bench.params.samples = 1
    bench.params.gctrial = false
    bench.params.gcsample = false
    return bench
end

foreach_benchmark(f!, bench::Benchmark) = f!(bench)
function foreach_benchmark(f!, group::BenchmarkGroup)
    for x in values(group)
        foreach_benchmark(f!, x)
    end
end

function setup_smoke()
    suite = setup()
    foreach_benchmark(set_smoke_params!, suite)
    return suite
end

function clear()
    Benchfib.clear()
end

end  # module TapirSchedulersBenchmarks


text:
function GLHull(points::Array{Float64,2},color::GL.Point4d)::GL.GLMesh
function GLHull(points::Array{Float64,2})::GL.GLMesh
function GLHull2d(points::Array{Float64,2})::GL.GLMesh # points by row
function GLHulls(V::Array{Float64,2},
function GLPolygon(V::Lar.Points,copEV::Lar.ChainOp,copFE::Lar.ChainOp)::GL.GLMesh
function GLPolygon(V::Lar.Points,EV::Lar.Cells)::GL.GLMesh
function GLLar2gl(V::Lar.Points, CV::Lar.Cells)::GL.GLMesh
function GLLines(points::Lar.Points,lines::Lar.Cells,color=COLORS[12])::GL.GLMesh
function GLPoints(points::Lar.Points,color=COLORS[12])::GL.GLMesh # points by row
function GLPolyhedron(V::Lar.Points, FV::Lar.Cells, T::GL.Matrix4=M44)::GL.GLMesh
function GLPolyhedron(V::Lar.Points,
function GLGrid(V::Lar.Points,CV::Lar.Cells,color=GL.COLORS[1])::GL.GLMesh
function GLExplode(V,FVs,sx=1.2,sy=1.2,sz=1.2,colors=1)


text:
@testset "tensorproduct" begin
    rng = MersenneTwister(123456)
    u1 = rand(rng, 10)
    u2 = rand(rng, 10)
    v1 = rand(rng, 5)
    v2 = rand(rng, 5)

    # kernels
    k1 = SqExponentialKernel()
    k2 = ExponentialKernel()
    kernel1 = TensorProduct(k1, k2)
    kernel2 = TensorProduct([k1, k2])

    @test kernel1 == kernel2
    @test kernel1.kernels === (k1, k2) === TensorProduct((k1, k2)).kernels
    @test length(kernel1) == length(kernel2) == 2
    @test_throws DimensionMismatch kernel1(rand(3), rand(3))

    @testset "val" begin
        for (x, y) in (((v1, u1), (v2, u2)), ([v1, u1], [v2, u2]))
            val = k1(x[1], y[1]) * k2(x[2], y[2])

            @test kernel1(x, y) == kernel2(x, y) == val
        end
    end

    # Standardised tests.
    TestUtils.test_interface(kernel1, ColVecs{Float64})
    TestUtils.test_interface(kernel1, RowVecs{Float64})
    test_ADs(
        ()->TensorProduct(SqExponentialKernel(), LinearKernel());
        dims = [2, 2],
    ) # ADs = [:ForwardDiff, :ReverseDiff])
    test_params(TensorProduct(k1, k2), (k1, k2))

    @testset "single kernel" begin
        kernel = TensorProduct(k1)
        @test length(kernel) == 1

        @testset "eval" begin
            for (x, y) in (((v1,), (v2,)), ([v1], [v2]))
                val = k1(x[1], y[1])

                @test kernel(x, y) == val
            end
        end
    end
end


text:
using Documenter, Knet, AutoGrad

# Avoid loading examples and their requirements
# Extract their docstrings instead
# Unfortunately this leads to wrong "source" links :(
# run(`make exampledocs`)
# include("exampledocs.jl")

# Load examples
# load_only = true
# for ex in ("linreg","housing","mnist","lenet","charlm","optimizers","vgg","resnet")
#     println("$ex.jl")
#     include(Knet.dir("examples","$ex.jl"))
# end
# println("Examples loaded")

makedocs(
    # Including modules prevents getting docstrings from Main
    # Including Main in the list leads to too many warnings
    # modules = [Knet,AutoGrad,LinReg,Housing,MNIST,LeNet,CharLM,Optimizers,VGG,ResNet],
    modules = [Knet,AutoGrad],
    clean = false,              # do we clean build dir
    #format = :html,
    sitename = "Knet.jl",
    authors = "Deniz Yuret and contributors.",
    doctest = true,
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        "Manual" => Any[
            "install.md",
            "tutorial.md",
#           "examples.md",
            "reference.md",
        ],
        "Textbook" => Any[
            "backprop.md",
            "softmax.md",
            "mlp.md",
            "cnn.md",
            "rnn.md",
            "rl.md",
            "opt.md",
            "gen.md",
        ],
    ],
#    analytics = "UA-89508993-1",
#    linkcheck = !("skiplinks" in ARGS),
)

deploydocs(
    repo = "github.com/denizyuret/Knet.jl.git",
    julia = "",
    osname = "",
    target = "build",
    make = nothing,
    deps = nothing,
    #deps   = Deps.pip("mkdocs", "python-markdown-math"),
)


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

# It has modified the definition of QR to support general τ

# QR and Hessenberg Factorizations
"""
    QR <: Factorization

A QR matrix factorization stored in a packed format, typically obtained from
[`qr`](@ref). If ``A`` is an `m`×`n` matrix, then

```math
A = Q R
```

where ``Q`` is an orthogonal/unitary matrix and ``R`` is upper triangular.
The matrix ``Q`` is stored as a sequence of Householder reflectors ``v_i``
and coefficients ``\\tau_i`` where:

```math
Q = \\prod_{i=1}^{\\min(m,n)} (I - \\tau_i v_i v_i^T).
```

Iterating the decomposition produces the components `Q` and `R`.

The object has two fields:

* `factors` is an `m`×`n` matrix.

  - The upper triangular part contains the elements of ``R``, that is `R =
    triu(F.factors)` for a `QR` object `F`.

  - The subdiagonal part contains the reflectors ``v_i`` stored in a packed format where
    ``v_i`` is the ``i``th column of the matrix `V = I + tril(F.factors, -1)`.

* `τ` is a vector  of length `min(m,n)` containing the coefficients ``\tau_i``.

"""
struct QR{T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}} <: Factorization{T}
    factors::S
    τ::Tau

    function QR{T,S,Tau}(factors, τ) where {T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}}
        require_one_based_indexing(factors)
        new{T,S,Tau}(factors, τ)
    end
end
QR(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T} = QR{T,typeof(factors), typeof(τ)}(factors, τ)
QR{T}(factors::AbstractMatrix, τ::AbstractVector) where {T} =
    QR(convert(AbstractMatrix{T}, factors), convert(AbstractVector{T}, τ))

QR(F::LinearAlgebra.QR) = QR(F.factors, F.τ)

# iteration for destructuring into components
Base.iterate(S::QR) = (S.Q, Val(:R))
Base.iterate(S::QR, ::Val{:R}) = (S.R, Val(:done))
Base.iterate(S::QR, ::Val{:done}) = nothing

function generic_qrfactUnblocked!(A::AbstractMatrix{T}) where {T}
    require_one_based_indexing(A)
    m, n = size(A)
    τ = zeros(T, min(m,n))
    for k = 1:min(m - 1 + !(T<:Real), n)
        x = view(A, k:m, k)
        τk = reflector!(x)
        τ[k] = τk
        reflectorApply!(x, τk, view(A, k:m, k + 1:n))
    end
    QR(A, τ)
end

qrfactUnblocked!(A::AbstractMatrix) = generic_qrfactUnblocked!(A)
qrfactUnblocked!(A::StridedMatrix{T}) where T<:BlasFloat = QR(LAPACK.geqrf!(A)...)

qrunblocked!(A::AbstractMatrix, ::Val{false}) = qrfactUnblocked!(A)
qrunblocked!(A::AbstractMatrix) = qrunblocked!(A, Val(false))

_qreltype(::Type{T}) where T = typeof(zero(T)/sqrt(abs2(one(T))))

function qrunblocked(A::AbstractMatrix{T}, arg) where T
    require_one_based_indexing(A)
    AA = similar(A, _qleltype(T), size(A))
    copyto!(AA, A)
    return qrunblocked!(AA, arg)
end
function qrunblocked(A::AbstractMatrix{T}) where T
    require_one_based_indexing(A)
    AA = similar(A, _qleltype(T), size(A))
    copyto!(AA, A)
    return qrunblocked!(AA)
end
qrunblocked(x::Number) = qrunblocked(fill(x,1,1))
function qrunblocked(v::AbstractVector)
    require_one_based_indexing(v)
    qrunblocked(reshape(v, (length(v), 1)))
end

# Conversions
QR{T}(A::QR) where {T} = QR(convert(AbstractMatrix{T}, A.factors), convert(AbstractVector{T}, A.τ))
Factorization{T}(A::QR{T}) where {T} = A
Factorization{T}(A::QR) where {T} = QR{T}(A)
AbstractMatrix(F::QR) = F.Q * F.R
AbstractArray(F::QR) = AbstractMatrix(F)
Matrix(F::QR) = Array(AbstractArray(F))
Array(F::QR) = Matrix(F)

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::QR)
    summary(io, F); println(io)
    println(io, "Q factor:")
    show(io, mime, F.Q)
    println(io, "\nR factor:")
    show(io, mime, F.R)
end

@inline function getR(F::QR, _)
    m, n = size(F)
    triu!(getfield(F, :factors)[1:min(m,n), 1:n])
end
@inline getQ(F::QR, _) = QRPackedQ(getfield(F, :factors), F.τ)

getL(F::QR, ::Tuple{AbstractVector,AbstractVector}) = getL(F, size(F))
getQ(F::QR, ::Tuple{AbstractVector,AbstractVector}) = getQ(F, size(F))

getR(F::QR) = getR(F, size(F.factors))
getQ(F::QR) = getQ(F, size(F.factors))

function getproperty(F::QR, d::Symbol)
    if d == :R
        return getR(F)
    elseif d == :Q
        return getQ(F)
    else
        getfield(F, d)
    end
end

Base.propertynames(F::QR, private::Bool=false) =
    (:R, :Q, (private ? fieldnames(typeof(F)) : ())...)

ldiv!(F::QR, B::AbstractVecOrMat) = ArrayLayouts.ldiv!(F, B)
ldiv!(F::QR, B::LayoutVector) = ArrayLayouts.ldiv!(F, B)
ldiv!(F::QR, B::LayoutMatrix) = ArrayLayouts.ldiv!(F, B)

"""
    QRPackedQ <: AbstractMatrix

The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QR`](@ref).
"""
struct QRPackedQ{T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}} <: LayoutQ{T}
    factors::S
    τ::Tau

    function QRPackedQ{T,S,Tau}(factors, τ) where {T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}}
        require_one_based_indexing(factors)
        new{T,S,Tau}(factors, τ)
    end
end
QRPackedQ(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T} = QRPackedQ{T,typeof(factors),typeof(τ)}(factors, τ)
function QRPackedQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T}
    QRPackedQ(convert(AbstractMatrix{T}, factors), convert(AbstractVector{T}, τ))
end

QRPackedQ{T}(Q::QRPackedQ) where {T} = QRPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(AbstractVector{T}, Q.τ))
QRPackedQ(Q::LinearAlgebra.QRPackedQ) = QRPackedQ(Q.factors, Q.τ)

AbstractMatrix{T}(Q::QRPackedQ{T}) where {T} = Q
AbstractMatrix{T}(Q::QRPackedQ) where {T} = QRPackedQ{T}(Q)

size(F::QR, dim::Integer) = size(getfield(F, :factors), dim)
size(F::QR) = size(getfield(F, :factors))

## Multiplication by Q
### QB

MemoryLayout(::Type{<:QRPackedQ{<:Any,S,T}}) where {S,T} =
    QRPackedQLayout{typeof(MemoryLayout(S)),typeof(MemoryLayout(T))}()

MemoryLayout(::Type{<:QR{<:Any,S,T}}) where {S,T} =
    QRPackedLayout{typeof(MemoryLayout(S)),typeof(MemoryLayout(T))}()

function (\)(A::QR{TA}, B::AbstractVecOrMat{TB}) where {TA,TB}
    require_one_based_indexing(B)
    S = promote_type(TA,TB)
    m, n = size(A)
    m == size(B,1) || throw(DimensionMismatch("Both inputs should have the same number of rows"))

    AA = Factorization{S}(A)

    X = _zeros(S, B, n)
    X[1:size(B, 1), :] = B
    ldiv!(AA, X)
    return _cut_B(X, 1:n)
end

function (\)(A::QR{T}, BIn::VecOrMat{Complex{T}}) where T<:BlasReal
    require_one_based_indexing(BIn)
    m, n = size(A)
    m == size(BIn, 1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(BIn,1)) rows"))

# |z1|z3|  reinterpret  |x1|x2|x3|x4|  transpose  |x1|y1|  reshape  |x1|y1|x3|y3|
# |z2|z4|      ->       |y1|y2|y3|y4|     ->      |x2|y2|     ->    |x2|y2|x4|y4|
#                                                 |x3|y3|
#                                                 |x4|y4|
    B = reshape(copy(transpose(reinterpret(T, reshape(BIn, (1, length(BIn)))))), size(BIn, 1), 2*size(BIn, 2))

    X = _zeros(T, B, n)
    X[1:size(B, 1), :] = B

    ldiv!(A, X)

# |z1|z3|  reinterpret  |x1|x2|x3|x4|  transpose  |x1|y1|  reshape  |x1|y1|x3|y3|
# |z2|z4|      <-       |y1|y2|y3|y4|     <-      |x2|y2|     <-    |x2|y2|x4|y4|
#                                                 |x3|y3|
#                                                 |x4|y4|
    XX = reshape(collect(reinterpret(Complex{T}, copy(transpose(reshape(X, div(length(X), 2), 2))))), _ret_size(A, BIn))
    return _cut_B(XX, 1:n)
end


text:
dev = device()

@test name(dev) isa String
@test uuid(dev) isa Base.UUID
totalmem(dev)
attribute(dev, CUDA.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
@test warpsize(dev) == 32
capability(dev)
@grab_output show(stdout, "text/plain", dev)

@test eval(Meta.parse(repr(dev))) == dev

@test eltype(devices()) == CuDevice
@grab_output show(stdout, "text/plain", CUDA.DEVICE_CPU)
@grab_output show(stdout, "text/plain", CUDA.DEVICE_INVALID)


text:
include("Scanner.jl")

struct SyntaxException <: Exception
  msg::String
end

abstract type Node end
abstract type Statement <: Node end
abstract type Value <: Node end

struct Statements <: Node
  statements::Array{Statement,1}
  line::Int
end

struct Declaration <: Statement
  variable::Token
  type::Token
  line::Int
end

struct DecAssignment <: Statement
  variable::Token
  type::Token
  value::Value
  line::Int
end

struct Assignment <: Statement
  variable::Token
  value::Value
  line::Int
end

struct For <: Statement
  variable::Token
  from::Value
  to::Value
  body::Statements
  line::Int
end

struct Read <: Statement
  variable::Token
  line::Int
end

struct Print <: Statement
  argument::Value
  line::Int
end

struct Assert <: Statement
  argument::Value
  line::Int
end

struct Literal <: Value
  token::Token
  line::Int
end

struct Operator <: Node
  token::Token
  line::Int
end

struct BinaryOperation <: Value
  leftOperand::Value
  operator::Token
  rightOperand::Value
  line::Int
end

struct UnaryOperation <: Value
  operator::Token
  operand::Value
  line::Int
end

struct Var <: Value
  variable::Token
  line::Int
end

struct LeftVal <: Node
  token::Token
  line::Int
end

"""
The main function of the parser. Maintains an index ("next") to point
at the token to be processed next. Builds the AST by calling mutually
recursive parsing functions, starting with statements().
"""
function parseInput(input::Array{Token,1})

  # Gets the next token.
  function nxtok()
    if next > length(input)
      throw(SyntaxException(
        "Failed to parse program. Do all your statements have a terminating semicolon?"))
    end
    return input[next]
  end

  nxtclass() = nxtok().class

  # Convenience function. Returns the token, its class and its line number.
  function tok_class_line()
    tok = nxtok()
    return tok, tok.class, tok.line
  end

  """
  Consumes the next token, checking that the token class matches.
  """
  function match_term(terminal::TokenClass)
    DEBUG && println("match_term called in $(currentUnit) with terminal $(terminal), next is $(nxtclass())")
    token, class, line = tok_class_line()
    does_match = terminal == class
    if !does_match
      if terminal == semicolon && class == eoi
        throw(SyntaxException(
          "Unexpected end of input. Did you forget a semicolon?"))
      elseif terminal == ident && class ∈ values(keywords)
        throw(SyntaxException(
          "Expected a variable identifier but got keyword \"$(token.lexeme)\" on line $(line)."))
      else
        throw(SyntaxException(
          "Failed to parse $(currentUnit) on or around line $(line)."))
      end
    end
    next += 1
    return does_match
  end

  function statements()
    currentUnit = "group of statements"
    token, class, line = tok_class_line()
    DEBUG && println("this is statements, nxtype is ", class)
    stmts = Array{Node,1}()
    while nxtclass() ∉ [eoi, kw_end]
      push!(stmts, statement())
      match_term(semicolon)
    end
    return Statements(stmts, line)
  end

  function statement()
    currentUnit = "statement"
    token, class, line = tok_class_line()
    DEBUG && println("this is statement, nxtype is ", class)
    if class == kw_var
      # Parse a declaration or a declaration-assignment
      match_term(kw_var)
      leftval_tok = leftval_ident().token
      match_term(colon)
      token, class, line = tok_class_line()
      if class ∉ [kw_bool, kw_int, kw_string]
        throw(SyntaxException(
          "Expected a type keyword on line $(line), got '$(token.lexeme)'."))
      end
      var_type_tok = token
      match_term(class)
      if nxtclass() == assign
        match_term(assign)
        return DecAssignment(leftval_tok, var_type_tok, expr(), line)
      end
      return Declaration(leftval_tok, var_type_tok, line)
    end
    if class == ident
      # Parse an assignment
      leftval_tok = leftval_ident().token
      line = nxtok().line
      match_term(assign)
      value = expr()
      return Assignment(leftval_tok, value, line)
    end
    if class == kw_for
      # Parse a for-loop statement
      match_term(kw_for)
      iter_var_token = leftval_ident().token
      match_term(kw_in)
      from = expr()
      match_term(rng)
      to = expr()
      match_term(kw_do)
      body = statements()
      match_term(kw_end)
      match_term(kw_for)
      return For(iter_var_token, from, to, body, line)
    end
    if class == kw_read
      # Parse a read statement
      match_term(kw_read)
      return Read(leftval_ident().token, line)
    end
    if class == kw_print
      # Parse a print statement
      match_term(kw_print)
      return Print(expr(), line)
    end
    if class == kw_assert
      # Parse an assert statement
      match_term(kw_assert)
      if nxtclass() != open_paren
        throw(SyntaxException(
          "The argument of an assert statement must be in parentheses (line $(line))."))
      end
      match_term(open_paren)
      argument = expr()
      match_term(close_paren)
      return Assert(argument, line)
    end
    throw(SyntaxException(
      "Failed to parse statement on line $(line)"))
  end

  function expr()
    currentUnit = "expression"
    token, class, line = tok_class_line()
    DEBUG && println("this is expr, nxtclass is ", nxtclass())
    if class ∈ keys(unary_ops)
      return UnaryOperation(unary_op().token, operand(), line)
    end
    leftOperand = operand()
    token, class, line = tok_class_line()
    if class ∈ keys(binary_ops)
      currentUnit = "binary operation"
      DEBUG && println("this is expr parsing binary op, nxtclass is ", class)
      operatorToken = operator().token
      rightOperand = operand()
      return BinaryOperation(leftOperand, operatorToken, rightOperand, line)
    end
    return leftOperand
  end

  function unary_op()
    currentUnit = "unary operation"
    token, class, line = tok_class_line()
    DEBUG && println("this is unary_op, nxtclass is ", class)
    match_term(class)
    return Operator(token, line)
  end

  function operand()
    currentUnit = "operand"
    token, class, line = tok_class_line()
    DEBUG && println("this is operand, nxtclass is ", class)
    if class == open_paren
      match_term(open_paren)
      expression = expr()
      match_term(close_paren)
      return expression
    end
    if class == ident
      match_term(class)
      return Var(token, line)
    end
    if class ∈ [int_literal, string_literal, kw_true, kw_false]
      match_term(class)
      return Literal(token, line)
    end
    throw(SyntaxException(
      "Expected an operand on line $(line), got '$(token.lexeme)'."
    ))
  end

  function operator()
    currentUnit = "operator"
    token, class, line = tok_class_line()
    DEBUG && println("this is operator, nxtclass is ", class)
    if class ∈ keys(binary_ops)
      match_term(class)
      return Operator(token, line)
    end
    throw(SyntaxException(
      "Expected a binary operator on line $(line), got '$(token.lexeme)'."
    ))
  end

  function leftval_ident()
    currentUnit = "identifier"
    token, class, line = tok_class_line()
    DEBUG && println("this is leftval_ident, nxtclass is ", class)
    match_term(ident)
    return LeftVal(token, line)
  end

  if length(input) < 2
    throw(SyntaxException(
      "The empty string is not a valid program."
    ))
  end
  next = 1
  currentUnit = "program"
  statements()
end

parseInput(source::String) = parseInput(scanInput(source))

max_stars_repo_name: juhamyllari/mini-pl-julia

max_issues_repo_name: juhamyllari/mini-pl-julia

max_forks_repo_name: juhamyllari/mini-pl-julia


text:
# p12.jl - accuracy of Chebyshev spectral differentiation
#         (compare p7.jl)

# Compute derivatives for various values of N:
Nmax = 50; allN = 6:2:Nmax; E = zeros(4,Nmax);
for N = 1:Nmax
    (D,x) = cheb(N);
    v = @. abs(x)^3; vprime = @. 3*x*abs(x);   # 3rd deriv in BV
    E[1,N] = norm(D*v-vprime,Inf);
    v = @. exp(-x^(-2)); vprime = @. 2*v/x^3;  # C-infinity
    E[2,N] = norm(D*v-vprime,Inf);
    v = @. 1/(1+x^2); vprime = @. -2*x*v^2;    # analytic in [-1,1]
    E[3,N] = norm(D*v-vprime,Inf);
    v = x.^10; vprime = 10*x.^9;               # polynomial
    E[4,N] = norm(D*v-vprime,Inf);
end

# Plot results:
titles = [L"|x|^3",L"\exp(-x^2)",L"1/(1+x^2)",L"x^{10}"];  clf();
for iplot = 1:4
    subplot(2,2,iplot);
    semilogy(1:Nmax,E[iplot,:],".-",markersize=6);
    axis([0,Nmax,1e-16,1e3]); grid(true);
    xticks(0:10:Nmax); yticks(10.0.^(-15:5:0));
    title(titles[iplot]);
    iplot > 2 ? xlabel("N") : nothing;
    iplot%2 > 0? ylabel("error") : nothing;
end


text:
module PyTSP

# Write your package code here.

using PyCall

include("tsp_solvers.jl")

export  solve_TSP_Concorde, solve_TSP_LKH,
        distance_matrix, tour_length, distance_matrix_float

end


text:
module PSSFSS

module GSMs
export GSM

mutable struct GSM{T1 <: AbstractMatrix, T2 <: AbstractMatrix}
    s11::T1
    s12::T2
    s21::T2
    s22::T1
end

end # module GSMs
##############################################
module Outputs
using ..GSMs: GSM
using JLD2
using StaticArrays
export runit

const SV2 = SVector{2,Float64}

SteerType = Union{NamedTuple{(:ψ₁, :ψ₂), Tuple{Float64, Float64}},
                  NamedTuple{(:θ, :ϕ), Tuple{Float64, Float64}}}

struct Result
    gsm::GSM
    steering::SteerType
    β⃗₀₀::SV2
    FGHz::Float64
    ϵᵣin::ComplexF64
    μᵣin::ComplexF64
    β₁in::SV2
    β₂in::SV2
    ϵᵣout::ComplexF64
    μᵣout::ComplexF64
    β₁out::SV2
    β₂out::SV2
end

function append_result_data(fname, gname, result)
    jldopen(fname, "a") do fid
        group = JLD2.Group(fid, gname)
        group["result"] = result
    end
    return
end

function read_result_file(fname::AbstractString)::Vector{Result}
    dat = load(fname) # a Dict
    ks = collect(keys(dat))
    sort!(ks, by = x -> parse(Int,split(x, '/')[1]))
    Result[dat[k] for k in ks]
end

function runit(n1=1)
    infile = "pssfss.res"
    outfile = "test.res"
    isfile(outfile) && rm(outfile)
    results = read_result_file(infile)
    for k = n1:length(results)
        println(k)
        gname = string(k)
        append_result_data(outfile, gname, results[k])
    end
end

end # module Outputs
##############################################

end # module PSSFSS


text:
using Test
using LinearAlgebra
using BQCESubroutine
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using BQCESubroutine: ctrl_offset, bcomspace, bsubspace, broutine2x2!, broutine!
using ThreadingUtilities

function naive_broutine!(st::AbstractVector, U, locs)
    n = log2dim(st)
    subspace = bsubspace(n, sort(locs))
    comspace = bcomspace(n, sort(locs))
    indices = [idx + 1 for idx in comspace]
    for k in subspace
        idx = indices .+ k
        st[idx] = U * st[idx]
    end
    return st
end

function naive_broutine!(st::AbstractMatrix, U, locs)
    n = log2dim(st)
    subspace = bsubspace(n, sort(locs))
    comspace = bcomspace(n, sort(locs))
    indices = [idx + 1 for idx in comspace]
    for b in 1:size(st, 1)
        for k in subspace
            idx = indices .+ k
            st[b, idx] = U * st[b, idx]
        end
    end
    return st
end

function naive_broutine!(st::AbstractVector, U, locs, ctrl)
    n = log2dim(st)
    subspace = bsubspace(n, sort(merge_locations(locs, ctrl.storage)))
    comspace = bcomspace(n, sort(locs))
    offset = ctrl_offset(ctrl)
    indices = [idx + 1 for idx in comspace]
    for k in subspace
        idx = indices .+ k .+ offset
        st[idx] = U * st[idx]
    end
    return st
end

function naive_broutine!(st::AbstractMatrix, U, locs, ctrl)
    n = log2dim(st)
    subspace = bsubspace(n, sort(merge_locations(locs, ctrl.storage)))
    comspace = bcomspace(n, sort(locs))
    offset = ctrl_offset(ctrl)
    indices = [idx + 1 for idx in comspace]
    for b in 1:size(st, 1)
        for k in subspace
            idx = indices .+ k .+ offset
            st[b, idx] = U * st[b, idx]
        end
    end
    return st
end

function naive_broutine_multi!(st, U, locs)
    for loc in locs
        naive_broutine!(st, U, loc)
    end
    return st
end

const N = 10

@testset "generic broutine" begin
    @testset "size(st)=$(size(st))" for st in [rand(ComplexF64, 1 << N), rand(ComplexF64, 5, 1 << N)]
        @testset "$(1<<M)x$(1<<M)" for M in [1,2,3,4]
            U = rand(ComplexF64, 1<<M, 1<<M)
            @testset "i=$i" for i in 1:N
                locs = Locations([mod1(i+k, N) for k in 0:M-1])
                locs = sort(locs)
                ctrl = CtrlLocations((mod1(i+M, N), mod1(i+M+2, N)))
                @test broutine!(copy(st), U, locs) ≈ naive_broutine!(copy(st), U, locs)
                @test broutine!(copy(st), U, locs, ctrl) ≈ naive_broutine!(copy(st), U, locs, ctrl)
            end
        end
    end
end

@testset "basic broutine $G" for G in [:X, :Y, :Z, :H, :T, :S, :Sdag, :Tdag]
    @testset "size(st)=$(size(st))" for st in [rand(ComplexF64, 1 << N), rand(ComplexF64, 5, 1 << N)]
        U = getfield(BQCESubroutine, Symbol(:B, G))

        @testset "$i=>$G" for i in 1:N
            locs = Locations(i)
            @test isapprox(broutine!(copy(st), Val(G), locs), naive_broutine!(copy(st), U, locs); rtol=1e-8)
        end

        @testset "($i, $(mod1(i+1, N)))=>$G" for i in 1:N
            locs = sort(Locations((i, mod1(i+2, N))))
            S1 = broutine!(copy(st), Val(G), locs)
            S2 = naive_broutine_multi!(copy(st), U, locs)
            @test isapprox(S1, S2; rtol=1e-8)
        end

        @testset "@ctrl $(mod1(i+1, N)) $i=>$G" for i in 1:N
            locs = Locations(i)
            ctrl = CtrlLocations(mod1(i+1, N))
            @test isapprox(broutine!(copy(st), Val(G), locs, ctrl), naive_broutine!(copy(st), U, locs, ctrl); rtol=1e-8)
        end
    end
end

@testset "parameterize gate $G" for G in [:Rx, :Ry, :Rz]
    @testset "size(st)=$(size(st))" for st in [rand(ComplexF64, 1 << N), rand(ComplexF64, 5, 1 << N)]
        U = getfield(BQCESubroutine, Symbol(:B, G))(0.2)

        @testset "$i=>$G" for i in 1:N
            locs = Locations(i)
            @test broutine!(copy(st), Val(G), locs, 0.2) ≈ naive_broutine!(copy(st), U, locs)
        end

        @testset "($i, $(mod1(i+1, N)))=>$G" for i in 1:N
            locs = sort(Locations((i, mod1(i+2, N))))
            @test broutine!(copy(st), Val(G), locs, 0.2) ≈ naive_broutine_multi!(copy(st), U, locs)
        end

        @testset "@ctrl $(mod1(i+1, N)) $i=>$G" for i in 1:N
            locs = Locations(i)
            ctrl = CtrlLocations(mod1(i+1, N))
            @test broutine!(copy(st), Val(G), locs, ctrl, 0.2) ≈ naive_broutine!(copy(st), U, locs, ctrl)
        end
    end
end

@testset "parameterize gate PSWAP" begin
    U = BQCESubroutine.BPSWAP(0.2)
    @testset "size(st)=$(size(st))" for st in [rand(ComplexF64, 1 << N), rand(ComplexF64, 5, 1 << N)]
        @testset "$((i, mod1(i+1, N)))=>PSWAP" for i in 1:N
            locs = Locations((i, mod1(i+1, N)))
            @test broutine!(copy(st), Val(:PSWAP), locs, 0.2) ≈ naive_broutine!(copy(st), U, locs)
        end

        @testset "@ctrl $(mod1(i+1, N)) $i=>PSWAP" for i in 1:N
            locs = sort(Locations((i, mod1(i+1, N))))
            ctrl = CtrlLocations(mod1(i+2, N))
            @test broutine!(copy(st), Val(:PSWAP), locs, ctrl, 0.2) ≈ naive_broutine!(copy(st), U, locs, ctrl)
        end
    end
end

@testset "gate SWAP" begin
    U = BQCESubroutine.BSWAP
    @testset "size(st)=$(size(st))" for st in [rand(ComplexF64, 1 << N), rand(ComplexF64, 5, 1 << N)]
        @testset "$((i, mod1(i+1, N)))=>SWAP" for i in 1:N
            locs = Locations((i, mod1(i+1, N)))
            @test broutine!(copy(st), Val(:SWAP), locs) ≈ naive_broutine!(copy(st), U, locs)
        end

        @testset "@ctrl $(mod1(i+1, N)) $i=>SWAP" for i in 1:N
            locs = sort(Locations((i, mod1(i+1, N))))
            ctrl = CtrlLocations(mod1(i+2, N))
            @test broutine!(copy(st), Val(:SWAP), locs, ctrl) ≈ naive_broutine!(copy(st), U, locs, ctrl)
        end
    end
end

@testset "diagonal" begin
    @testset "dtype_st=$(dtype_st), size(st)=$(size(st))" for dtype_st in [Float64, ComplexF64], st in [rand(dtype_st, 1 << N), rand(dtype_st, 5, 1 << N)]
        @testset "dtype_U=$(dtype_U), $(1<<M) x $(1<<M)" for dtype_U in [Float64, ComplexF64], M in [1,2,3,4]
            dtype_st == Float64 && dtype_U == ComplexF64 && break
            U = Diagonal(rand(dtype_U, 1<<M))
            U_dense = Matrix(U)
            @testset "i=$i" for i in 1:N
                locs = Locations([mod1(i+k, N) for k in 0:M-1])
                locs = sort(locs)
                ctrl = CtrlLocations((mod1(i+M, N), mod1(i+M+2, N)))
                @test broutine!(copy(st), U, locs) ≈ broutine!(copy(st), U_dense, locs)
                @test broutine!(copy(st), U, locs, ctrl) ≈ broutine!(copy(st), U_dense, locs, ctrl)
            end
        end
    end
end


text:
using GigaFV

configfile = "./src/advection.yaml"
GigaFV.main(configfile)


text:
module SomeModule
import TestSubgraphTrackingPath
using Test

@test TestSubgraphTrackingPath.addfive(1) == 6

end


text:
ProjDir = joinpath(dirname(@__FILE__), "..", "examples", "Bernoulli")
cd(ProjDir) do

  isdir("tmp") &&
    rm("tmp", recursive=true);

  include(joinpath(ProjDir, "bernoulli_nochains.jl"))

  isdir("tmp") &&
    rm("tmp", recursive=true);

end # cd


text:
# %% Sensititivy to trajectory length
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

# %%
N = 1*10^5
systems = [:koch, :henon_chaotic]
slabels = ["Koch", "Hénon"]
qs = 2:4
Cmethod = "standard" # bueno or standard. Decides εmax for correlation sum.

eHs, eCs, Hs, Cs = [Vector{Float64}[] for i in 1:4]
for data in systems
    for q in qs
        qH = qC = q
        @show (q, data)
        # Here we simply pack all parameters into a dictionary
        # (other parameters are (probably) globals)
        params = @strdict N qH qC data
        if data == :standardmap_chaotic
            params["k"] = 1.0
        # elseif data == :henon_chaotic
            # params["z"] = -4
        end
        params["theiler"] = 0
        if Cmethod ≠ "standard"
            params["Cmethod"] = Cmethod
        end

        # This is the main call that calculates everything
        output, s = produce_or_load(
            datadir("main"), params, make_C_H;
            prefix = string(data), suffix = "jld2", force = false,
            ignores = ["data"], storepatch = false
        )

        @unpack eH, eC, H, C = output
        push!(eHs, eH); push!(Hs, H); push!(eCs, eC); push!(Cs, C)
    end
end

legendtitle = "impact of order \$q\$"
labels = [s*" \$q=$(q)\$" for s in slabels for q in qs]

fig, axs = mainplot(
    Hs, Cs, eHs, eCs, labels, legendtitle;
    qH = "q", qC = "q", tol = 0.25,
    offsets = range(0; length = 6, step = 1.5),
    dimension_fit_C = FractalDimension.linear_regression_fit_glm,
)

wsave(plotsdir("paper", "orderq"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using Distributions
export LazyΩ
using Random
using ..IDS, ..Util, ..Tagging, ..RNG, ..Space
import ..Basis: replacetags
import ..Var

"Lazily constructs randp, values as they are needed"
struct LazyΩ{TAGS <: Tags, T} <: AbstractΩ
  data::T
  tags::TAGS
end

# AbstractDict interface

Base.keys(ω::LazyΩ) = Base.keys(ω.data)
Base.values(ω::LazyΩ) = Base.values(ω.data)
# Move to specific omega types
Base.merge!(ω::LazyΩ, ω_) = error("Unimplemented")
Basis.like(ω::LazyΩ{Tags, T}, kv::Pair) where {Tags, T} = T(kv)

const EmptyTags = Tags{(),Tuple{}}
LazyΩ{EmptyTags, T}() where T = LazyΩ(T(), Tags())

"Construct `LazyΩ` from `rng` -- `ω.data` will be generated from `rng`"
LazyΩ{T}(rng::AbstractRNG) where T = tagrng(LazyΩ{T}(), rng)

Basis.idtype(ω::LazyΩ{TAGS, Dict{T, V}}) where {TAGS, T, V} = T
ids(ω::LazyΩ) = keys(ω.data)

replacetags(ω::LazyΩ, tags) = LazyΩ(ω.data, tags)
traits(::Type{LazyΩ{TAGS, T}}) where {TAGS, T} = traits(TAGS)

Base.setindex!(ω::LazyΩ, value, id) =
  ω.data[convertid(idtype(ω), id)] = value

function (exo::Var.ExoRandVar)(ω::LazyΩ)
  result = get(ω.data, exo, 0)
  if result === 0
    ω.data[exo] = rand(rng(ω), exo.class)
  else
    result
  end::eltype(exo.class)
end

# function resolve(dist, id, ω::LazyΩ)
#   id_ = convertid(idtype(ω), id)
#   # id_ = id
#   if haskey(ω.data,  id_)
#     d, val = ω.data[id_]
#   else
#     ω.data[id_] = (dist, rand(rng(ω), dist))
#     d, val = ω.data[id_]
#   end
#   val::eltype(dist)
# end


text:
#using NamedTupleTools

function add_tuple_sum(x::NT, y::NT) where { NT<:NamedTuple }
    @assert all(fieldtypes(x) .== fieldtypes(y))
    names_of_values = propertynames(x)
    @assert propertynames(y)  == names_of_values
    values_of_x = fieldvalues(x)
    values_of_y = fieldvalues(y)
    #nt = namedtuple(names_of_values, values_of_x .+ values_of_y )
    values = map(t->t[1].+t[2], zip(values_of_x,values_of_y))
    nt = namedtuple(names_of_values, values )
    return nt
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
@recipe function f(det::SolidStateDetector{T}; n = 30, φ = missing, seriescolor = missing, label = missing) where {T}

    clabel = (ismissing(label) ? map(c -> (c.name == "" ? c.id : c.name), det.contacts) : label)
    if !(typeof(clabel) <: AbstractArray) clabel = [clabel] end
    ccolor = (ismissing(seriescolor) ? map(c -> c.id, det.contacts) : seriescolor)
    if !(typeof(ccolor) <: AbstractArray) ccolor = [ccolor] end

    if ismissing(φ)
        xguide --> "x / m"
        yguide --> "y / m"
        zguide --> "z / m"

        for (cn,contact) in enumerate(det.contacts)
            linewidth --> 2
            seriescolor := ccolor[(cn-1)%length(ccolor)+1]
            @series begin
                label := ""
                n --> n
                contact
            end
            @series begin
                label := clabel[(cn-1)%length(clabel)+1]
                seriescolor := ccolor[(cn-1)%length(ccolor)+1]
                []
            end
        end
    else
        @info "2D plotting of the detector not yet implemented."
    end
end

@recipe function f(contact::Contact{T}; n = 30, seriescolor = missing) where {T}
    ccolor = (ismissing(seriescolor) ? contact.id : seriescolor)
    @series begin
        seriescolor := ccolor
        label --> (contact.name == "" ? contact.id : contact.name)
        []
    end
    for (cn,c) in enumerate(contact.geometry_positive)
        @series begin
            seriescolor := ccolor
            label := ""
            n --> n
            c
        end
    end
end

#=
@recipe function f(geometry::AbstractGeometry)
    pos, neg = get_decomposed_volumes(geometry)
    for g in pos
        @series begin
            g
        end
    end
end
@recipe function f(object::AbstractObject)
    @series begin
        object.geometry
    end
end

## 2D from Felix

function line_2d(r1,r2,z1,z2)
   x1::Vector = [r1,z1]
   x2::Vector = [r2,z2]
   return [x1[1],x2[1]],[x1[2],x2[2]]
end

function polar_circle_2d(r::T, phiStart::T, phiStop::T;nSteps=360) where T <: SSDFloat
    if r == 0; return [],[]; end
    phiRange = collect(phiStart:(phiStop-phiStart)/nSteps:phiStop)
    rRange = [r for a in 1:length(phiRange)]
    return phiRange, rRange
end

@recipe function f(d::SolidStateDetector{T}, dim::Symbol; φ = missing, z = missing) where{T <: SSDFloat}
    if ismissing(z) && ismissing(φ)
        print("Please specify φ or z.")
    elseif ismissing(z)
        for c in d.contacts
            for g in c.geometry_positive
                @series begin
                    if d.name == "Public Inverted Coax"
                        if typeof(c) == Contact{T,:N}; color --> :orange
                        elseif typeof(c) == Contact{T,:P}; color --> :blue; end
                    else color --> :black
                    end
                    lw --> 2
                    label = ""
                    g, :φ, T(deg2rad(mod(φ,2π)))
                end
            end
        end
    elseif ismissing(φ)
        proj --> :polar
        for c in d.contacts
            for g in c.geometry_positive
                @series begin
                    if d.name == "Public Inverted Coax"
                        if typeof(c) == Contact{T,:N}; color --> :orange
                        elseif typeof(c) == Contact{T,:P}; color --> :blue; end
                    else color --> :black
                    end
                    lw --> 2
                    label = ""
                    g, :z, T(z)
                end
            end
        end
    else
        print("Please specify only one of either φ or z.")
    end
end

@recipe function f(Vol::Tube{T}, dim::Symbol, parameter::T) where{T <: SSDFloat}
    if dim == :φ
        if parameter in Vol.φ_interval
           rStart = Vol.r_interval.left
           rStop = Vol.r_interval.right
           zStart = Vol.z_interval.left
           zStop = Vol.z_interval.right

           @series begin
                   label --> ""
                   line_2d(rStop,rStop,zStart,zStop)
               end

            @series begin
                   label --> ""
                   line_2d(rStart,rStop,zStart,zStart)
            end

           if rStart != rStop
                @series begin
                   label --> ""
                   line_2d(rStart,rStop,zStop,zStop)
               end
            end

            if zStart != zStop && rStart != 0
                @series begin
                   label --> ""
                   line_2d(rStart,rStart,zStart,zStop)
               end
            end
        end

    elseif dim == :z
        if parameter in Vol.z_interval
            rStart = Vol.r_interval.left
            rStop = Vol.r_interval.right
            phiStart = Vol.φ_interval.left
            phiStop = Vol.φ_interval.right
            zStart = Vol.z_interval.left
            zStop = Vol.z_interval.right

            if rStart != 0
                @series begin
                    label --> ""
                    polar_circle_2d(rStart,phiStart,phiStop)
                end
            end
            @series begin
                label --> ""
                polar_circle_2d(rStop,phiStart,phiStop)
            end

            if !isapprox(mod(phiStart,2π),mod(phiStop,2π),atol = 1e-5)
                @series begin
                    label --> ""
                    line_2d(phiStart,phiStart,rStart,rStop)
                end
                @series begin
                    label --> ""
                    line_2d(phiStop,phiStop,rStart,rStop)
                end
            end
        end
    end
end

# @recipe function f(Vol::ConeMantle{T}, dim::Symbol, parameter::T) where{T <: SSDFloat}
#     newVol = Vol.cone
#     @series begin
#         newVol, dim, parameter
#     end
# end

@recipe function f(Vol::Cone{T}, dim::Symbol, parameter::T) where{T <: SSDFloat}
    if dim == :φ
        if parameter in Vol.φ_interval
            rStart = Vol.r_interval.left
            rStop = Vol.r_interval.right
            zStart = Vol.z_interval.left
            zStop = Vol.z_interval.right
            orientation = Vol.orientation

            if orientation in [:top_right, :bottom_left]
                @series begin
                    label --> ""
                   line_2d(rStop,rStart,zStart,zStop)
                end
            else
                @series begin
                    label --> ""
                   line_2d(rStart,rStop,zStart,zStop)
                end
            end
        end

    elseif dim == :z
        if parameter in Vol.z_interval
            rStart = Vol.r_interval.left
            rStop = Vol.r_interval.right
            phiStart = Vol.φ_interval.left
            phiStop = Vol.φ_interval.right
            zStart = Vol.z_interval.left
            zStop = Vol.z_interval.right
            orientation = Vol.orientation

            if orientation in [:top_right, :bottom_left]
                @series begin
                    label --> ""
                    r = (rStop-rStart)*(parameter-zStop)/(zStart-zStop) + rStart
                    polar_circle_2d(r,phiStart,phiStop)
                end
            else
                @series begin
                    label --> ""
                    r = (rStop-rStart)*(parameter-zStart)/(zStop-zStart) + rStart
                    polar_circle_2d(r,phiStart,phiStop)
                end
            end
        end

    end
end
=#


text:
using FixedSizeArrays

import Base.*

export M22, M23, M24, M32, M33, M34, M42, M43, M44, M2, M3, M4, V2, V3, V4

typealias M22 Mat{2, 2, Float32}
typealias M23 Mat{2, 3, Float32}
typealias M24 Mat{2, 4, Float32}
typealias M32 Mat{3, 2, Float32}
typealias M33 Mat{3, 3, Float32}
typealias M34 Mat{3, 4, Float32}
typealias M42 Mat{4, 2, Float32}
typealias M43 Mat{4, 3, Float32}
typealias M44 Mat{4, 4, Float32}

typealias M2 M22
typealias M3 M33
typealias M4 M44

typealias V2 Vec{2, Float32}
typealias V3 Vec{3, Float32}
typealias V4 Vec{4, Float32}

max_stars_repo_name: twadleigh/aminda

max_issues_repo_name: twadleigh/aminda

max_forks_repo_name: twadleigh/aminda


text:
#
# This file is part of the Guards.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

#
# Indexing
#
function Base.getindex(G::Guard{T}, i::Int) where T<:AbstractArray
    @grd getindex(G, i)
end

function Base.setindex!(G::Guard{T}, v, i::Int) where T<:AbstractArray
    @grd setindex!(G, v, i)
end

function Base.firstindex(G::Guard{T}) where T<:AbstractArray
    @grd firstindex(G)
end

function Base.lastindex(G::Guard{T}) where T<: AbstractArray
    @grd lastindex(G)
end

#
# AbstractArrays
#
Base.supertype(::Type{Guard{T}}) where T<: AbstractArray = T

function Base.size(G::Guard{T}) where T<: AbstractArray
    @grd size(G)
end

function Base.IndexStyle(::Type{Guard{T}}) where T<: AbstractArray
    IndexStyle(T)
end

function Base.getindex(G::Guard{T}, I::Vararg{Int, N}) where {T<: AbstractArray, N}
    @grd getindex(G, I)
end

function Base.setindex!(G::Guard{T}, v, I::Vararg{Int, N}) where {T<: AbstractArray, N}
    @grd setindex!(G, v, I)
end

function Base.length(G::Guard{T}) where T <: AbstractArray
    @grd length(G)
end

#
# Dicts
#
function Base.getindex(G::Guard{T}, key) where T <: AbstractDict
    @grd getindex(G, key)
end

function Base.setindex!(G::Guard{T}, value, key) where T <: AbstractDict
    @grd setindex!(G, value, key)
end


text:
module HackerNews

# Module Imports
using HTTP
using JSON

# Abstract Type
abstract type HNApiRoute end

# HNApiRoute Subtypes
struct MaxItem <: HNApiRoute end
struct TopStories <: HNApiRoute end
struct NewStories <: HNApiRoute end
struct BestStories <: HNApiRoute end
struct AskStories <: HNApiRoute end
struct ShowStories <: HNApiRoute end
struct JobStories <: HNApiRoute end
struct Updates <: HNApiRoute end

# Declaring the values
hn_route_string(::Type{MaxItem}) = "maxitem"
hn_route_string(::Type{TopStories}) = "topstories"
hn_route_string(::Type{NewStories}) = "newstories"
hn_route_string(::Type{BestStories}) = "beststories"
hn_route_string(::Type{AskStories}) = "askstories"
hn_route_string(::Type{ShowStories}) = "showstories"
hn_route_string(::Type{JobStories}) = "jobstories"
hn_route_string(::Type{Updates}) = "updates"

"""
HackerNews initializer composite type
"""
struct HN
    story
    nposts::Int64
    user_related::Bool

    # Inner constructor
    HN(story, nposts=1, user_related=false) = new(story, nposts, user_related)
end

"""
HackerNews post composite type
"""
struct HNPost
    data
    id
    deleted
    itemtype
    by
    time
    text
    dead
    parent
    poll
    kids
    url
    score
    title
    parts
    descendants

    # Inner constructor
    HNPost(data,
            id      =   get(data,"id",nothing),
            deleted =   get(data,"deleted",nothing),
            itemtype=   get(data, "type", nothing),
            by      =   get(data,"by",nothing),
            time    =   get(data,"time",nothing),
            text    =   get(data,"text",nothing),
            dead    =   get(data,"dead",nothing),
            parent  =   get(data,"parent",nothing),
            poll    =   get(data,"poll",nothing),
            kids    =   get(data, "kids",nothing),
            url     =   get(data,"url",nothing),
            score   =   get(data,"score",nothing),
            title   =   get(data,"title",nothing),
            parts   =   get(data,"parts",nothing),
            descendants=get(data,"descendants",nothing)) = new(data,id, deleted, itemtype, by, time,text, dead,
                                                            parent, poll, kids, url, score, title, parts, descendants)
end

"""
Hackernews user composite type
"""
struct HNUser
    data
    id
    delay
    created
    karma
    about
    submitted

    # Inner constructor
    HNUser(data,
            id      =   get(data, "id", nothing),
            delay   =   get(data, "delay", nothing),
            created =   get(data, "created", nothing),
            karma   =   get(data, "karma", nothing),
            about   =   get(data, "about", nothing),
            submitted=  get(data, "submitted", nothing)) = new(data, id, delay, created, karma, about, submitted)
end

"""
Parse response for updates
"""
function parse_hn_response(::Type{Updates}, response, hn)
    @info("fetching updates..")
    JSON.parse(String(response.body))[hn.user_related ? "profiles" : "items"]
end

"""
Parse response for maxitem
"""
function parse_hn_response(::Type{MaxItem}, response, hn)
    @info("fetching maxitem..")
    JSON.parse(String(response.body))
end

"""
Parse response
"""
function parse_hn_response(story, response, hn)
    @info("fetching $(story)...")
    eval(JSON.parse(String(response.body)))[1:hn.nposts]
end

"""
function to get the ids for a story
"""
function getids(hn::HN)
    response = HTTP.request("GET","https://hacker-news.firebaseio.com/v0/" * hn_route_string(hn.story) * ".json?print=pretty")
    parse_hn_response(hn.story, response, hn)
end

"""
function to get the actual posts for given array of IDs
"""
function getinfo(hn::HN)
    ids = HackerNews.getids(hn)
    posts = HNPost[]
    users = HNUser[]

    if hn.user_related == false
        @info("generating post data...")
        for i in 1:hn.nposts
            response = HTTP.request("GET","https://hacker-news.firebaseio.com/v0/item/" * string(ids[i]) * ".json?print=pretty")
            response_json = JSON.parse(String(response.body))
            hnpost = HNPost(response_json)
            push!(posts, hnpost)
        end
    else
        for i in 1:hn.nposts
            @info("generating user data...")
            response = HTTP.request("GET","https://hacker-news.firebaseio.com/v0/user/" * string(ids[i]) * ".json?print=pretty")
            response_json = JSON.parse(String(response.body))
            hnuser = HNUser(response_json)
            push!(users, hnuser)
        end
    end

    # return whats required
    if hn.user_related users else posts end
end

"""
function to get a single user info
"""
function getuser(id::String)
    response = HTTP.request("GET","https://hacker-news.firebaseio.com/v0/user/" * id * ".json?print=pretty")
    response_json = JSON.parse(String(response.body))
    HNUser(response_json)
end

# Module exports
export getinfo, getuser, getids
export HN, HNUser, HNPost
export HNApiRoute, MaxItem, TopStories, NewStories, BestStories
export AskStories, ShowStories, JobStories, Updates
export hn_route_string
end