text:
export TRARC
function TRARC(nlp :: AbstractNLPModel,
               x₀ :: Array{Float64,1},
               TR :: TrustRegion,
			   c :: Combi;
               atol :: Float64 = 1e-8,
               rtol :: Float64 = 1.0e-6,
               max_eval :: Int = 50_000,
               max_iter :: Int = 10_000,
               #max_f :: Int=5000,
               #max_eval :: Int = 5000,
               robust :: Bool = true,
               verbose :: Bool = true
               )
    hessian_rep,PData,solve_model,pre_process,decrease,params = extract(c)

    α = TR.α₀
    x, xnext, d, Df = copy(x₀), copy(x₀), copy(x₀), 0.0
    xopt = x
    λ = 1.0

    n = length(x)
	g = Vector{Float64}(undef, n)
	gnext = Vector{Float64}(undef, n)

    f = obj(nlp,x)
    if isnan(f) | (f==Inf)
        OK = false
	println("f nan or ∞");
	xopt=f
	gopt=f
	fopt=f
        niter = Inf
        calls = [Inf, Inf, Inf, Inf]
	return xopt,gopt,gopt,fopt,niter,calls,OK,:Bad_x0
    end

    fopt = f
    grad!(nlp,x,g)
    norm_g = stop_norm(g)
    norm_g0 = norm(g)
    gopt = g
    H = hessian_rep(nlp,x)

    optimal = (norm_g < atol) || (isinf(f) & (f<0.0))
    tired = false
    stalled = false
    finish = optimal || tired || stalled

    fnext = f
    iter = 0
    verbose && display_header_iterations()
    verbose && display_success(iter,fnext,norm_g0,0.0,α)

    succ, unsucc, verysucc, unsuccinarow = 0, 0, 0, 0

    max_unsuccinarow = TR.max_unsuccinarow

    calls = [0, 0, 0, 0]

    while ~finish

        PData = pre_process(H,g,params,calls,max_eval)

        if ~PData.OK
			return x,fopt,norm_g,Inf,Inf,[Inf,Inf,Inf,Inf],false,:PreFailed;
		end

        success = false

#        while ~success & (iter < max_iter) & (unsuccinarow < TR.max_unsuccinarow)
        while ~success & (iter < max_iter) & (~stalled)
            try
                d,λ = solve_model(PData,α)
            catch
                println(" Problem in solve_model")

                return x,fopt,norm_g,norm(gopt),iter,calls,false,:AscentDir
            end

            Δq = -(g + 0.5 * H * d)⋅d

            if Δq < 0.0 println("*******   Ascent direction in SolveModel: Δq = $Δq")
                println("  g⋅d = $(g⋅d), 0.5 d'Hd = $(0.5*(H*d)⋅d)  α = $α  λ = $λ")
                println("  calls  $calls   iters $verysucc; $succ; $unsucc")

                return x,fopt,norm_g,norm(gopt),iter,calls,false,:AscentDir
            end
            slope = g ⋅ d
            xnext = x + d

            iter += 1

            # Trap illegal values for f, descent will backtrack until a legal value is reached
            try
                fnext = obj(nlp,xnext)
            catch
                fnext = Inf;
            end
            if isnan(fnext)
                fnext=Inf
            end

            Δf = f - fnext

            r, good_grad, gnext = compute_r(nlp,f,Δf,Δq,slope,d,xnext,gnext,robust)

            if r<TR.acceptance_threshold
                verbose && display_failure(iter,fnext,λ,α)
	        unsucc=unsucc+1
	        unsuccinarow = unsuccinarow +1
                stalled = unsuccinarow >= max_unsuccinarow
	        try
                    α = decrease(PData, α, TR)
                catch
                    stalled = true
                end
                #fbidon = obj(nlp,x)
	    else
	        success = true

	        unsuccinarow = 0
	        x = copy(xnext)
	        f = fnext

	        if good_grad  g = copy(gnext)
                else grad!(nlp,x,g);
                end

                norm_g=stop_norm(g)
                H = hessian_rep(nlp,x)
	        if r > TR.increase_threshold
	            α = increase(PData, α, TR)
	            verbose && display_v_success(iter,f,norm_g,λ,α)
                    verysucc += 1
	        else
                    if r < TR.reduce_threshold
                        α = decrease(PData, α, TR)
                    end
	            verbose && display_success(iter,f,norm_g,λ,α)
	            succ += 1
	        end
	    end
        end
        if isnan(f)  OK=false
	    println("f nan");
	    xopt=f
	    gopt=f
	    fopt=f
            niter = Inf
            calls = [Inf, Inf, Inf, Inf]
	    return xopt,gopt,fopt,niter,calls,OK
        end
        calls = [nlp.counters.neval_obj,  nlp.counters.neval_grad, nlp.counters.neval_hess, nlp.counters.neval_hprod]

        optimal = (norm_g < atol)| (norm_g <( rtol * norm_g0)) | (isinf(f) & (f<0.0))
        tired = (iter >= max_iter) | (sum(calls) > max_eval)

        finish = (optimal || tired || stalled)
    end

    xopt = x
    fopt = f
    gopt = g

    niter = verysucc+succ+unsucc;
    verbose && @printf(" %10i unsuccessful iterations \n",unsucc)
    verbose && @printf(" %10i successful iterations \n",succ)
    verbose && @printf(" %10i very successful iterations \n",verysucc)
    verbose && @printf(" %10i total successful iterations \n",verysucc+succ)
    verbose && @printf(" %10i total iterations \n",niter)

    OK = optimal

    calls = [nlp.counters.neval_obj, nlp.counters.neval_grad, nlp.counters.neval_hess, nlp.counters.neval_hprod]

    total_calls = nlp.counters.neval_obj + nlp.counters.neval_grad + n*nlp.counters.neval_hess+ nlp.counters.neval_hprod
     verbose && @printf(" %10i total f %10i g %10i H calls  %10i Hv calls  \n",nlp.counters.neval_obj, nlp.counters.neval_grad, nlp.counters.neval_hess, nlp.counters.neval_hprod)

    if OK && verbose
		@printf("\n     Final cost(x) = %f  ||g|| = %9.2e  number of f and g and H evals = %i  converged\n",fopt, norm(gopt),total_calls)
    elseif verbose
        @printf("\n     Final cost(x) = %f  ||g|| = %9.2e  number of f and g and H evals  = %i  NOT converged\n\n",fopt, norm(gopt),total_calls)
    end;

    if OK status = :Optimal
    elseif stalled status = :stalled
    else status = :UserLimit
    end

    return xopt,fopt,norm_g,norm(gopt),niter,calls,OK,status

end


text:
using Pkg; Pkg.activate(pwd())
using Test, LinearAlgebra, Random, JLD2

include(joinpath(pwd(), "src\\ops_utility.jl"))
include(joinpath(pwd(), "src\\ops_build.jl"))
include(joinpath(pwd(), "src\\ops_methods.jl"))
include(joinpath(pwd(), "src\\cfrops_solve.jl"))

@testset "Build Ops Game" begin
    # Build game and infosets Random.seed!(579843)
    g = AirDefenseGame()
    A, An, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    Σ, seqn, seqactions = build_Σset(g, A, ns1, ns2)
    Iset, (I1set, I2set), (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    @test length(Iset) == 1695
    gs = GameSet(g)
    # @test seqI2[631] == 7 # sequence [3,4] => [5,8] for Attacker: targetable [1, 4, 5, 8], cyberattack [5, 8], cyberdefend [1, 3]

    # History to infoset
    ndp, nac, ndc = g.ndp, g.nac, g.ndc
    id1 = infoset(SVector(2, 2, 2, 4, 0, 0), 5, g, gs)
    @test Iset[id1] == "Defender: defendable [1, 3, 4, 8], sensed [1, 4]"
    id2 = infoset(SVector(15, 3, 5, 4, 0, 0), 5, g, gs)
    @test Iset[id2] == "Defender: defendable [4, 5, 8, 9], sensed [3, 9]"
    id3 = infoset(SVector(10, 4, 6, 4, 1, 0), 6, g, gs)
    @test Iset[id3] == "Attacker: targetable [1, 3, 5, 8], cyberattack [5, 8], cyberdefend [1, 5]"

    # chance distributions: stage 1
    An = gs.An
    h0 = SVector(0, 0, 0, 0, 0, 0)
    @test all(isapprox(sum(chance(1, a, g, gs) for a in An[1]), 1.0, atol = 0.00001))
    @test all(chance(1, a, g, gs) >= 0.0 for a in An[1])
    # chance distributions: stage 2
    @test all(isapprox(sum(chance(2, a, g, gs) for a in An[2]), 1.0, atol = 0.00001))
    @test all(chance(2, a, g, gs) >= 0.0 for a in An[2], h1 in An[1])
    # chance distributions: stage 4
    @test all(isapprox(sum(chance(4, a, g, gs) for a in An[4]), 1.0, atol = 0.00001))
    @test all(chance(4, a, g, gs) >= 0.0 for a in An[4])
end

@testset "Build Ops Utility Tests (Loads Local Files)" begin
    # Build utility
    g = AirDefenseGame()
    gs = GameSet(g)
    A, An, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    Σ, seqn, seqactions = build_Σset(g, A, ns1, ns2)
    # U, z = build_utility_hist(g, A, An)
    # reward = build_utility_seq(g, gs, seqn, seqactions) # takes an 1 hr to build...
    dir = pwd()
    fn = joinpath(dir, "data\\vars_opt_rewardfunction.jld2")
    @load fn U z reward
    @test length(U) == 3_888_000
    @test length(U) == nnz(reward)

    ind1 = 43
    ind2 = 62491
    h = [1, 6, 5, 1, 6, 120] # h associated with seq1, seq2
    indu = findfirst(x -> x == h, z)
    @test U[indu] == reward[ind1, ind2]
end

@testset "Ops CFR" begin
    # Build game and infosets Random.seed!(579843)
    g = AirDefenseGame()
    gs = GameSet(g)

    # check that strategy builds (still need correctness test)
    T = 100
    u, r, s, σ = cfr(T, g, gs)
    @test all(all(isnan.(σi)) || all(σi .>= 0.0) for σi in σ)
    @test all(all(isnan.(σi)) || isapprox(sum(σi), 1.0, atol = 0.001) for σi in σ)

    gtoy = AirDefenseGame(2, [1.0 1;2 1], [1.0, 1], 0.2, 2, 1, 2, 2, 1, 1, 0.9, 0.5, 0.7, 0.8, [1,2])
    gstoy = GameSet(gtoy)
    T = 15_000
    (utoy, r, s, σtoy), t_cfr = @timed cfr(T, gtoy, gstoy)
    @test isapprox(σtoy[infoset(SVector(1,1,1,1,0,0), 5, gtoy, gstoy)][1], 0.5, atol = 0.1)
    @test isapprox(σtoy[infoset(SVector(1,1,1,2,0,0), 5, gtoy, gstoy)][1], 1.0, atol = 0.001)
end

ext: jl

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
import accrue.cryptoerase.runtime.Condition;

class C {
  static final Condition{H} c = new Condition();
  int x = 7;
  int{L /c H} y = 42;
  void m() {
    c.set();
    this.x = this.y;
  }

  public static void main(String[] args) {
    new C().m();
  }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
const CoeffDict{N,T} = Dict{ SVector{N,Int}, T }

"""Lazy Taylor series with N variables and coefficients of type T
`degree` is the max degree of monomials allowed.
If `degree == -1` then the Tayor series is potentially infinite."""
struct Taylor{N,T,F}
    f::F
    coeffs::CoeffDict{N,T}
    degree::Int   # maximum order of monomial
end

function Taylor(N, T, f::F, degree=-1) where {F}
    t = Taylor{N, T, F}(f, CoeffDict{N,T}(), degree)
    dummy = t[SVector(ntuple(_->0, Val(N)))]  # compile getindex by calculating first coefficient
    return t
end

degree(x::SVector) = sum(x)
degree(f::Taylor) = f.degree

import Base: ^
^(t::Taylor{N,T,F}, n::Integer) where {N,T,F} = Base.power_by_squaring(t, n)

function getindex(t::Taylor{N,T,F}, index::SVector) where {N,T,F}

    if degree(index) > t.degree
        return zero(T)
    end

    coeffs = t.coeffs

    if haskey(coeffs, index)
        return coeffs[index]
    end

    coeffs[index] = (t.f)(t, index)

    return coeffs[index]

end

function getindex(t::Taylor{N,T,F}, index...) where {N,T,F}
    getindex(t, SVector(index))
end

function degree_sum(f, g)
    if f.degree == -1 || g.degree == -1
        return -1  # infinite
    end

    return  f.degree + g.degree
end

function degree_max(f, g)
    if f.degree == -1 || g.degree == -1
        return -1  # infinite
    end

    return max(f.degree, g.degree)
end

# these are memoized, but should look at performance without memoizing perhaps
+(f::Taylor{N,T}, g::Taylor{N,T}) where {N,T} = Taylor(N, T,
                                                (t, i) -> f[i] + g[i],
                                                degree_max(f, g) )

-(f::Taylor{N,T}, g::Taylor{N,T}) where {N,T} = Taylor(N, T,
                                                        (t, i) -> f[i] - g[i],
                                                        degree_max(f, g) )

#
# # formulas from Warwick Tucker, *Validated Numerics*
#
function *(f::Taylor{N,T}, g::Taylor{N,T}) where {N,T}

    return Taylor(N, T,
            (t, index) -> begin

                coeff = zero(T)

                tuples = generate_tuples(index)

                for i in tuples
                    coeff += f[i] * g[index - i]
                end

                coeff
            end,
            degree_sum(f, g)
            )
end

    #     for monomial_f in keys(f.coeffs)
    #         for monomial_g in keys(g.coeffs)
    #
    #             monomial_product = monomial_f + monomial_g
    #
    #
    #

# enumerate tuples less than or equal to t, varying over dimension which

"Generate tuples that are lexicographically less than x"
function generate_tuples(x::SVector{N,Int}) where N

    t = zero(x)

    tuples = [t]

    @inbounds while t != x
        #@show t
        which = N

        while t[which] == x[which] && which > 1
            t = setindex(t, 0, which)
            which -= 1
        end

        t = setindex(t, t[which]+1, which)
        push!(tuples, t)
    end

    return tuples

end

"Generate tuples with a given maximum degree"
function generate_tuples(N, deg)

    t = zero(SVector{N,Int})

    tuples = [t]

    which = N

    @inbounds while degree(t) <= deg && which > 0

        if degree(t) == deg
            t = setindex(t, 0, which)
            which -= 1

            if which == 0
                break
            end

        else
            which = N
        end

        t = setindex(t, t[which]+1, which)
        push!(tuples, t)
    end

    return tuples

end
#
# struct TupleGenerator{N,T}
#     x::SVector{N,T}
# end
#
# function Base.start(g::TupleGenerator{N,T}) where {N,T}
#     return zero(g.x)
# end
#
# function Base.next(g::TupleGenerator{N,T}, state) where {N,T}
#     current = state
#
#     which = N
#
#     while state[which] == g.x[which] && which > 1
#         state = setindex(state, 0, which)
#         which -= 1
#     end
#
#     state = setindex(state, state[which]+1, which)
#
#     return current, state
# end
#
# function Base.done(g::TupleGenerator{N,T}, state) where {N,T}
#     return state == g.x
# end

# enum_tuples(SVector(1,1), 1)

# self is a reference to the object exp(g), that is used recursively
# function exp(g::Taylor)
#     function f(self, k)
#         k == 0 && return exp(g[0])
#         # dummy = g[k]  # preallocate g
#         return sum(i * g[i] * self[k-i] for i in 1:k) / k
#     end
#
#     return Taylor(f)
#
# end
#
# end

const variable_names = ["x", "y", "z", "w"]

function evaluate!(f::Taylor{N,T}) where {N,T}
    tuples = generate_tuples(N, degree(f))

    for t in tuples
        f[t]
    end
end

function Base.show(io::IO, f::Taylor{N,T}) where {N, T}

    if degree(f) == 0
        print(io, f[zero(SVector{N,Int})])
        return
    end

    # todo: order by degree
    for t in sort(collect(keys(f.coeffs)), lt=isless)
        value = f[t]

        if value == 0
            continue
        end

        if (value == 1 && iszero(t))
            print(io, value)
        end

        if value != 1
            print(io, value)
        end

        for i in 1:N

            if t[i] == 1
                print(io, variable_names[i])
            elseif t[i] > 1
                print(io, variable_names[i], "^", t[i])
            end
        end

        print(io, " + ")
    end
end

export x, y, z, o

function variable(N, T, i)
    vec = SVector(ntuple(j->(j==i ? 1 : 0), N))
    return Taylor(N, T, (t,index)->(index==vec ? 1 : 0), 1)
end

function constant(N, T)
    return Taylor(N, T, (t,i)->(i==zero(SVector{N,T}) ? 1 : 0), 0)
end

function variables(N, T)
    vars = [variable(N, T, i) for i in 1:N]

    for var in vars
        evaluate!(var)
    end

    return (constant(N,T), vars ...)
end

export o, x, y, z
o, x, y, z = variables(3, Float64)

max_stars_repo_name: dpsanders/LazyTaylorSeries

max_issues_repo_name: dpsanders/LazyTaylorSeries

max_forks_repo_name: dpsanders/LazyTaylorSeries


text:
import Base.InteractiveUtils
@inline uses_affinemap(::Type{T}) where {T} = supports_affine(T) || supports_affineview(T)

@inline supports_eager(::Type) = true
@inline supports_affine(::Type{<:AffineOperation}) = true
@inline supports_affine(::Type) = false
@inline supports_affineview(::Type{T}) where {T} = supports_affine(T)
@inline supports_permute(::Type) = false
@inline supports_view(::Type) = false
@inline supports_stepview(::Type) = false
@inline supports_lazy(::Type{T}) where {T} = supports_affine(T) || supports_affineview(T) || supports_stepview(T) || supports_view(T) || supports_permute(T)

@inline supports_eager(A)      = supports_eager(typeof(A))
@inline supports_affine(A)     = supports_affine(typeof(A))
@inline supports_affineview(A) = supports_affineview(typeof(A))
@inline supports_permute(A)    = supports_permute(typeof(A))
@inline supports_view(A)       = supports_view(typeof(A))
@inline supports_stepview(A)   = supports_stepview(typeof(A))
@inline supports_lazy(A)       = supports_lazy(typeof(A))

# --------------------------------------------------------------------

"""
    prepareaffine(img)

Make sure `img` is either a `InvWarpedView` or a `SubArray` of
one. If that is already the case, `img` will be returned as is.
Otherwise `invwarpedview` will be called using a `Flat()`
extrapolation scheme.

Doing this will tell subsequent operations that they should also
participate as affine operations (i.e. use `AffineMap` if they
can).
"""
prepareaffine(img) = invwarpedview(img, toaffinemap(NoOp(), img), Flat())
prepareaffine(img::AbstractExtrapolation) = invwarpedview(img, toaffinemap(NoOp(), img))
@inline prepareaffine(img::SubArray{T,N,<:InvWarpedView}) where {T,N} = img
@inline prepareaffine(img::InvWarpedView) = img
prepareaffine(imgs::Tuple) = map(prepareaffine, imgs)

# currently unused
@inline preparelazy(img) = img

# --------------------------------------------------------------------
# Operation and AffineOperation fallbacks

# discard param if operations doesn't need it
@inline randparam(op::Operation, img) = nothing
@inline toaffinemap(op::Operation, img, param) = toaffinemap(op, img)

for FUN in (:applyeager, :applylazy, :applypermute,
            :applyaffine, :applyaffineview,
            :applyaffine_common, :applyaffineview_common,
            :applyview, :applystepview)
    @eval begin
        function ($FUN)(op::Operation, imgs::Tuple)
            param = randparam(op, imgs)
            map(img -> ($FUN)(op, img, param), imgs)
        end
        @inline function ($FUN)(op::Operation, img::AbstractArray)
            ($FUN)(op, img, randparam(op, img))
        end
    end
end

function applyeager(op::Operation, img::AbstractArray, param)
    maybe_copy(applylazy(op, img, param))
end

function applyaffineview(op::Operation, img::AbstractArray, param)
    wv = applyaffine(op, img, param)
    direct_view(wv, indices(wv))
end

function applyaffine(op::AffineOperation, img::AbstractArray, param)
    invwarpedview(img, toaffinemap(op, img, param))
end

# Allow affine operations to omit specifying a custom
# "applylazy". On the other hand this also makes sure that a
# custom implementation of "applylazy" is preferred over
# "applylazy_fallback" which by default just calls "applyaffine".
function applylazy(op::AffineOperation, img::AbstractArray, param)
    _applylazy(op, img, param)
end

# The purpose of having a separate "_applylazy" is to not
# force "applylazy" implementations to specify the type of "img",
# when typeof(img) <: AbstractArray.
function _applylazy(op::AffineOperation, img::InvWarpedView, param)
    applyaffine(op, img, param)
end

function _applylazy(op::AffineOperation, img::SubArray{T,N,<:InvWarpedView}, param) where {T,N}
    applyaffine(op, img, param)
end

function _applylazy(op::AffineOperation, img, param)
    applylazy_fallback(op, img, param)
end

# Defining "applylazy_fallback" instead of "applylazy" will
# make sure that the custom implementation is only used if
# "img" is not already an "InvWarpedView", in which case
# "applyaffine" would be called instead of "applylazy_fallback".
function applylazy_fallback(op::AffineOperation, img, param)
    applyaffine(op, prepareaffine(img), param)
end

# --------------------------------------------------------------------
# Since an "AffineMap" is strongly typed, the resulting type of
# "toaffinemap" can be very different from one operation to the
# next. To allow for type-stable stochastic branches (i.e. "Either")
# we need a way to force a common "AffineMap" type using only
# "SArray" internally (i.e. no "RotMatrix" or other special types).

@generated function toaffinemap_common(op::AffineOperation, img::AbstractArray{T,N}, param) where {T,N}
    quote
        tfm = toaffinemap(op, img, param)
        AffineMap(SMatrix(tfm.linear), SVector(tfm.translation))::AffineMap{SArray{Tuple{$N,$N},Float64,$N,$(N*N)},SVector{$N,Float64}}
    end
end

function applyaffine_common(op::AffineOperation, img::AbstractArray, param)
    invwarpedview(img, toaffinemap_common(op, img, param))
end

function applyaffineview_common(op::AffineOperation, img::AbstractArray, param)
    wv = applyaffine_common(op, img, param)
    direct_view(wv, indices(wv))
end

# We trust that non-affine operations use SArray-only AffineMap.
function applyaffineview_common(op::Operation, img::AbstractArray, param)
    applyaffineview(op, img, param)
end

# --------------------------------------------------------------------
# Functions to unroll sequences of Operation. These are called by
# the pipeline logic to prefer specific behaviour.

# "applyaffine" or "applyaffineview" is accepted
@inline @generated function unroll_applyaffine(op::Operation, img)
    if supports_affine(op)
        :(applyaffine(op, img))
    elseif supports_affineview(op)
        :(applyaffineview(op, img))
    else
        :(throw(MethodError(unroll_applyaffine, (op, img))))
    end
end

@inline unroll_applylazy(op::Operation, img) = applylazy(op, img)

# Recursively unroll the sequentially applied operations.
# This will result in a function that is equivalent to hardcoded
# chaining of operations.
for KIND in (:affine, :lazy)
    FUN = Symbol(:unroll_apply, KIND)
    PRE = Symbol(:prepare, KIND)
    @eval begin
        function ($FUN)(operations::NTuple{N,Operation}, img) where N
            ($FUN)(first(operations), Base.tail(operations), ($PRE)(img))
        end

        @inline function ($FUN)(head::Operation, tail::Tuple, img)
            ($FUN)(first(tail), Base.tail(tail), ($FUN)(head, img))
        end

        @inline function ($FUN)(head::Operation, tail::Tuple{}, img)
            ($FUN)(head, img)
        end
    end
end


text:
using DashBootstrapComponents, DashHtmlComponents

card_content_1 = [
    dbc_cardheader("Card header"),
    dbc_cardbody([
        html_h5("Card title", className="card-title"),
        html_p("This is some card content that we'll reuse", className="card-text"),
    ]),
];

card_content_2 = dbc_cardbody([
    html_blockquote(
        [
            html_p(
                "A learning experience is one of those things that says, " *
                "'You know that thing you just did? Don't do that.'",
            ),
            html_footer(html_small("Douglas Adams", className="text-muted")),
        ],
        className="blockquote",
    ),
]);

card_content_3 = [
    dbc_cardimg(src="/static/images/placeholder286x180.png", top=true),
    dbc_cardbody([
        html_h5("Card with image", className="card-title"),
        html_p(
            "This card has an image on top, and a button below",
            className="card-text",
        ),
        dbc_button("Click me!", color="primary"),
    ]),
];

cards = dbc_cardcolumns([
    dbc_card(card_content_1, color="primary", inverse=true),
    dbc_card(card_content_2, body=true),
    dbc_card(card_content_1, color="secondary", inverse=true),
    dbc_card(card_content_3, color="info", inverse=true),
    dbc_card(card_content_1, color="success", inverse=true),
    dbc_card(card_content_1, color="warning", inverse=true),
    dbc_card(card_content_1, color="danger", inverse=true),
    dbc_card(card_content_3, color="light"),
    dbc_card(card_content_1, color="dark", inverse=true),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# print AST in Lisp form

export pexpression, pexpressions, rmlines
"""
    pexpression(quoted expression)

Print readable form in Lisp format `(header arg1 arg2 ...)`.
"""
function pexpression(io::IO, ex::Expr)
    print(io, "(", ex.head)
    for arg in ex.args
        p0 = position(io)
        print(io, " ")
        pexpression(io, arg)
        if position(io) <= p0 + 1
            seek(io, p0)
        end
    end
    print(io, ")")
end

pexpression(io::IO, s::AbstractString) = print(io, '"', s, '"')

function pexpression(io::IO, s::GlobalRef)
    print(io,"(|.| ", string(s.mod), " ", s.name, ")")
end

function pexpression(io::IO, s::QuoteNode)
    print(io, "(quote "); pexpression(io, s.value); print(io, ")")
end

pexpression(io::IO, s) = print(io, s)

pexpression(io::IO, ::LineNumberNode) = nothing

pexpression(s) = pexpression(stdout, s)

pexpressions(s) = sprint(pexpression, s)

"""
    rmlines(expression)

Return new Expr with all LineiNumberNode's removed
"""
function rmlines(ex::Expr)
    args = filter(x -> !(x isa LineNumberNode), rmlines.(ex.args))
    Expr(ex.head, args...)
end
rmlines(x) = x


text:
"""
    lineplot(x, y; kw...)

# Description

Draws a path through the given points on a new canvas.

The first (optional) vector `x` should contain the horizontal positions for all the points along the path.
The second vector `y` should then contain the corresponding vertical positions respectively.
This means that the two vectors must be of the same length and ordering.

# Usage

    lineplot([x], y; $(keywords(; add = (:canvas,)))

    lineplot(fun, [start], [stop]; kw...)

# Arguments

$(arguments(
    (
        fun = "a unary function ``f: R -> R`` that should be evaluated, and drawn as a path from `start` to `stop` (numbers in the domain)",
        head_tail = "pass the symbol `:head`, `:tail` or `:both` to color the head/tail of the line with the complement of the chosen color",
        x = "horizontal position for each point (can be a real number or of type `TimeType`), if omitted, the axes of `y` will be used as `x`",
    ) ; add = (:x, :y, :canvas)
))

# Returns

A plot object of type `Plot{T<:Canvas}`.

# Author(s)

- Christof Stocker (github.com/Evizero)
- Milktrader (github.com/milktrader)

# Examples

```julia-repl
julia> lineplot([1, 2, 7], [9, -6, 8], title = "My Lineplot")
                      My Lineplot
       ┌────────────────────────────────────────┐
    10 │⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠│
       │⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠁⠀│
       │⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠁⠀⠀⠀⠀│
       │⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠔⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠤⠤⠤⠼⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⠤⠶⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
       │⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⡠⠔⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⢇⡠⠔⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   -10 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       └────────────────────────────────────────┘
       1                                        7
```

# See also

[`Plot`](@ref), [`scatterplot`](@ref), [`stairs`](@ref),
[`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
function lineplot(
    x::AbstractVector,
    y::AbstractVector;
    canvas::Type = KEYWORDS.canvas,
    color::UserColorType = KEYWORDS.color,
    name::AbstractString = KEYWORDS.name,
    head_tail::Union{Nothing,Symbol} = nothing,
    kw...,
)
    idx = map(x, y) do i, j
        isfinite(i) && isfinite(j)
    end
    plot = Plot(x[idx], y[idx], canvas; kw...)
    lineplot!(plot, x, y; color = color, name = name, head_tail = head_tail)
end

lineplot(y::AbstractVector; kw...) = lineplot(axes(y, 1), y; kw...)

function lineplot!(
    plot::Plot{<:Canvas},
    x::AbstractVector,
    y::AbstractVector;
    color::UserColorType = KEYWORDS.color,
    name::AbstractString = KEYWORDS.name,
    head_tail::Union{Nothing,Symbol} = nothing,
)
    color = color == :auto ? next_color!(plot) : color
    name == "" || label!(plot, :r, string(name), color)
    lines!(plot, x, y, color)
    (head_tail === nothing || length(x) == 0 || length(y) == 0) && return plot
    head_tail_color = (col = crayon_256_color(color)) === nothing ? nothing : ~col
    if head_tail in (:head, :both)
        points!(plot, last(x), last(y), head_tail_color)
    end
    if head_tail in (:tail, :both)
        points!(plot, first(x), first(y), head_tail_color)
    end
    plot
end

lineplot!(plot::Plot{<:Canvas}, y::AbstractVector; kw...) =
    lineplot!(plot, axes(y, 1), y; kw...)

# date time

function lineplot(
    x::AbstractVector{D},
    y::AbstractVector;
    xlim = extrema(x),
    kw...,
) where {D<:TimeType}
    d = Dates.value.(x)
    dlim = Dates.value.(D.(xlim))
    plot = lineplot(d, y; xlim = dlim, kw...)
    label!(plot, :bl, string(xlim[1]), color = :light_black)
    label!(plot, :br, string(xlim[2]), color = :light_black)
end

function lineplot!(
    plot::Plot{<:Canvas},
    x::AbstractVector{<:TimeType},
    y::AbstractVector;
    kw...,
)
    d = Dates.value.(x)
    lineplot!(plot, d, y; kw...)
end

# slope and intercept

function lineplot!(plot::Plot{<:Canvas}, intercept::Number, slope::Number; kw...)
    xmin = origin_x(plot.graphics)
    xmax = origin_x(plot.graphics) + width(plot.graphics)
    ymin = origin_y(plot.graphics)
    ymax = origin_y(plot.graphics) + height(plot.graphics)
    lineplot!(
        plot,
        [xmin, xmax],
        [intercept + xmin * slope, intercept + xmax * slope];
        kw...,
    )
end

# plotting a function

function lineplot(
    f::Function,
    x::AbstractVector;
    name = KEYWORDS.name,
    xlabel = "x",
    ylabel = "f(x)",
    kw...,
)
    y = Float64[f(i) for i in x]
    name = name == "" ? string(nameof(f), "(x)") : name
    plot = lineplot(x, y; name = name, xlabel = xlabel, ylabel = ylabel, kw...)
end

function lineplot(
    f::Function,
    startx::Number,
    endx::Number;
    out_stream::Union{Nothing,IO} = nothing,
    width::Int = out_stream_width(out_stream),
    kw...,
)
    diff = abs(endx - startx)
    x = startx:(diff / (3 * width)):endx
    lineplot(f, x; width = width, kw...)
end

lineplot(f::Function; kw...) = lineplot(f, -10, 10; kw...)

function lineplot!(plot::Plot{<:Canvas}, f::Function, x::AbstractVector; name = "", kw...)
    y = Float64[f(i) for i in x]
    name = name == "" ? string(nameof(f), "(x)") : name
    lineplot!(plot, x, y; name = name, kw...)
end

function lineplot!(
    plot::Plot{<:Canvas},
    f::Function,
    startx::Number = origin_x(plot.graphics),
    endx::Number = origin_x(plot.graphics) + width(plot.graphics);
    kw...,
)
    diff = abs(endx - startx)
    x = startx:(diff / (3ncols(plot.graphics))):endx
    lineplot!(plot, f, x; kw...)
end

# plotting vector of functions

lineplot(F::AbstractVector{<:Function}; kw...) = lineplot(F, -10, 10; kw...)

lineplot(F::AbstractVector{<:Function}, startx::Number, endx::Number; kw...) =
    _lineplot(F, startx, endx; kw...)

lineplot(F::AbstractVector{<:Function}, x::AbstractVector; kw...) = _lineplot(F, x; kw...)

function _lineplot(F::AbstractVector{<:Function}, args...; color = :auto, name = "", kw...)
    n = length(F)
    n > 0 || throw(ArgumentError("Can not plot empty array of functions"))
    color_is_vec = color isa AbstractVector
    name_is_vec  = name isa AbstractVector
    color_is_vec && (
        length(color) == n || throw(
            DimensionMismatch(
                "\"color\" must be a symbol or same length as the function vector",
            ),
        )
    )
    name_is_vec && (
        length(name) == n || throw(
            DimensionMismatch(
                "\"name\" must be a string or same length as the function vector",
            ),
        )
    )
    tcolor = color_is_vec ? color[1] : color
    tname  = name_is_vec ? name[1] : name
    plot   = lineplot(F[1], args...; color = tcolor, name = tname, kw...)
    for i in 2:n
        tcolor = color_is_vec ? color[i] : color
        tname  = name_is_vec ? name[i] : name
        lineplot!(plot, F[i], args...; color = tcolor, name = tname)
    end
    plot
end


text:
module Blackbody

using PhysicalConstants

export  PerHertz,
        PerAngstrom,
        planck,
        wien,
        stefan_boltzmann

const c0  = PhysicalConstants.CGS.PlancksConstantH * PhysicalConstants.CGS.SpeedOfLight
const c1h = 2.0 * c0 * PhysicalConstants.CGS.SpeedOfLight
const c1a = c1h / 1.0e8
const c2  = c0 / PhysicalConstants.CGS.Boltzmann

const b1 = 1.0e8 * c2 / 2.821439372
const b2 = 1.0e8 * c2 / 4.965114231744276

const sp = PhysicalConstants.CGS.StefanBoltzmannConstant / pi

abstract PerUnit

immutable PerHertz <: PerUnit
end

immutable PerAngstrom <: PerUnit
end

"""
Planck (blackbody) spectral radiance or specific intensity in CGS units.

    Args:
        wavelength (Array{Number,1} or Number): Wavelength(s) in Angstroms.
        temperature (Number): Temperature in Kelvin.
        output (Type{PerUnit}): Output convention (PerHertz or PerAngstrom).

    Returns:
        Array{Float64,1} or Float64: Planck spectral radiance in CGS units.
"""
function planck{T1<:Number,T2<:PerUnit}( wavelength::Array{T1,1}, temperature::Number, output::Type{T2} )
    [ planck( wavelength[ i ], temperature, output ) for i = 1:length( wavelength ) ]
end

function planck( wavelength::Number, temperature::Number, output::Type{PerHertz} )
    x = 1.0e8 / wavelength
    c1h * x ^ 3 / expm1( c2 * x / temperature )
end

function planck( wavelength::Number, temperature::Number, output::Type{PerAngstrom} )
    x = 1.0e8 / wavelength
    c1a * x ^ 5 / expm1( c2 * x / temperature )
end

"""
Wavelength at the maximum Planck spectral radiance for a given temperature
(Wien's displacement law).

    Args:
        temperature (Number): Temperature in Kelvin.
        output (Type{PerUnit}): Unit convention (PerHertz or PerAngstrom).

    Returns:
        Float64: Wavelength at maximum Planck spectral radiance in Angstroms.
"""
wien( temperature::Number, output::Type{PerHertz} ) = b1 / temperature

wien( temperature::Number, output::Type{PerAngstrom} ) = b2 / temperature

"""
Stefan-Boltzmann law.

    Args:
        temperature (Number): Temperature in Kelvin.

    Returns:
        Float64: Blackbody radiance in CGS units (power over pi).
"""
stefan_boltzmann( temperature::Number ) = sp * temperature ^ 4

end

# using Blackbody
#
# wavelengths = collect( logspace( 3, 4, 100000000 ) )
# for ( wl, bb ) in zip( wavelengths, planck( wavelengths, 10000.0, PerAngstrom ) )
#     println( "$wl $bb" )
# end
# println()
#
# wlmax = wien( 10000.0, PerAngstrom )
# println( wlmax, " ", planck( wlmax, 10000.0, PerAngstrom ) )
# println()
#
# @time planck( wavelengths, 10000.0, PerAngstrom )
# @time planck( wavelengths, 10000.0, PerAngstrom )
# @time planck( wavelengths, 10000.0, PerAngstrom )
# @time planck( wavelengths, 10000.0, PerAngstrom )
# @time planck( wavelengths, 10000.0, PerAngstrom )


text:
export sass_option_get_output_style, sass_option_set_output_style
export sass_option_get_source_comments, sass_option_set_source_comments
export sass_option_get_source_map_file, sass_option_set_source_map_file
export sass_option_get_omit_source_map_url, sass_option_set_omit_source_map_url
export sass_option_get_source_map_embed, sass_option_set_source_map_embed
export sass_option_get_source_map_contents, sass_option_set_source_map_contents
export sass_option_get_source_map_file_urls, sass_option_set_source_map_file_urls
export sass_option_get_source_map_root, sass_option_set_source_map_root
export sass_option_get_is_indented_syntax_src, sass_option_set_is_indented_syntax_src
export sass_option_get_include_path, sass_option_set_include_path
export sass_option_get_plugin_path, sass_option_set_plugin_path
export sass_option_get_indent, sass_option_set_indent
export sass_option_get_linefeed, sass_option_set_linefeed
export sass_option_get_input_path, sass_option_set_input_path
export sass_option_get_output_path, sass_option_set_output_path
export sass_option_get_precision, sass_option_set_precision

@enum Style nested expanded compact compressed

function sass_option_get_output_style(options)
    i = ccall((:sass_option_get_output_style, libsass_so),
        Cint, (Ptr{Cvoid},), options)
    Style(i)
end

sass_option_set_output_style(options, output_style::Style) =
    sass_option_set_output_style(options, Int(output_style))

function sass_option_set_output_style(options, output_style::Integer)
    ccall((:sass_option_set_output_style, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, output_style)
end

function sass_option_get_source_comments(options)
    ccall((:sass_option_get_source_comments, libsass_so),
        Cint, (Ptr{Cvoid},), options)
end

function sass_option_set_source_comments(options, source_comments)
    ccall((:sass_option_set_source_comments, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, source_comments)
end

function sass_option_get_source_map_file(options)
    s = ccall((:sass_option_get_source_map_file, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_option_set_source_map_file(options, map_file)
    ccall((:sass_option_set_source_map_file, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, map_file)
end

function sass_option_get_omit_source_map_url(options)
    ccall((:sass_option_get_omit_source_map_url, libsass_so),
        Cint, (Ptr{Cvoid},), options)
end

function sass_option_set_omit_source_map_url(options, omit_source_map_url)
    ccall((:sass_option_set_omit_source_map_url, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, omit_source_map_url)
end

function sass_option_get_source_map_embed(options)
    ccall((:sass_option_get_source_map_embed, libsass_so),
        Cint, (Ptr{Cvoid},), options)
end

function sass_option_set_source_map_embed(options, source_map_embed)
    ccall((:sass_option_set_source_map_embed, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, source_map_embed)
end

function sass_option_get_source_map_contents(options)
    ccall((:sass_option_get_source_map_contents, libsass_so),
        Cint, (Ptr{Cvoid},), options)
end

function sass_option_set_source_map_contents(options, source_map_contents)
    ccall((:sass_option_set_source_map_contents, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, source_map_contents)
end

function sass_option_get_source_map_file_urls(options)
    ccall((:sass_option_get_source_map_file_urls, libsass_so),
        Cint, (Ptr{Cvoid},), options)
end

function sass_option_set_source_map_file_urls(options, source_map_file_urls)
    ccall((:sass_option_set_source_map_file_urls, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, source_map_file_urls)
end

function sass_option_get_source_map_root(options)
    s = ccall((:sass_option_get_source_map_root, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_option_set_source_map_root(options, source_map_root)
    ccall((:sass_option_set_source_map_root, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, source_map_root)
end

function sass_option_get_is_indented_syntax_src(options)
    ccall((:sass_option_get_is_indented_syntax_src, libsass_so),
        Cint, (Ptr{Cvoid},), options)
end

function sass_option_set_is_indented_syntax_src(options, is_indented_syntax_src)
    ccall((:sass_option_set_is_indented_syntax_src, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, is_indented_syntax_src)
end

function sass_option_get_include_path(options)
    s = ccall((:sass_option_get_include_path, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL && return nothing
    julia_pathlist(unsafe_string(s))
end

function sass_option_set_include_path(options, include_path)
    ccall((:sass_option_set_include_path, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, c_pathlist(include_path))
end

function sass_option_get_plugin_path(options)
    s = ccall((:sass_option_get_plugin_path, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL && return nothing
    julia_pathlist(unsafe_string(s))
end

function sass_option_set_plugin_path(options, plugin_path)
    ccall((:sass_option_set_plugin_path, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, c_pathlist(plugin_path))
end

function sass_option_get_indent(options)
    s = ccall((:sass_option_get_indent, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_option_set_indent(options, indent)
    ccall((:sass_option_set_indent, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, indent)
end

function sass_option_get_linefeed(options)
    s = ccall((:sass_option_get_linefeed, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_option_set_linefeed(options, linefeed)
    ccall((:sass_option_set_linefeed, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, linefeed)
end

function sass_option_get_input_path(options)
    s = ccall((:sass_option_get_input_path, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_option_set_input_path(options, input_path)
    ccall((:sass_option_set_input_path, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, input_path)
end

function sass_option_get_output_path(options)
    s = ccall((:sass_option_get_output_path, libsass_so),
        Cstring, (Ptr{Cvoid},), options)
    s == C_NULL ? nothing : unsafe_string(s)
end

function sass_option_set_output_path(options, output_path)
    ccall((:sass_option_set_output_path, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cstring), options, output_path)
end

function sass_option_get_precision(options)
    ccall((:sass_option_get_precision, libsass_so),
        Cint, (Ptr{Cvoid},), options)
end

function sass_option_set_precision(options, precision)
    ccall((:sass_option_set_precision, libsass_so),
        Cvoid, (Ptr{Cvoid}, Cint), options, precision)
end

const separator = Sys.iswindows() ? ';' : ':'

c_pathlist(s::AbstractString) = s
c_pathlist(s) = join(s, separator)

julia_pathlist(s) = split(s, separator)


text:
@testset "617.merge-two-binary-trees.jl" begin
    @test merge_trees(
        TreeNode{Int}([1, 3, 2, 5]), TreeNode{Int}([2, 1, 3, nothing, 4, nothing, 7])
    ) == TreeNode{Int}([3, 4, 5, 5, 4, nothing, 7])
    @test merge_trees(TreeNode{Int}([1]), TreeNode{Int}([1, 2])) == TreeNode{Int}([2, 2])
    @test merge_trees(nothing, nothing) === nothing
end


text:
@testset "doc string generation" begin
    IterationControl.@create_docs(SquareRooter,
                              header="header",
                              example="example",
                              body="body")

    paragraphs = split(string(@doc SquareRooter), "\n")
    @test paragraphs[2] == "header"
    @test paragraphs[end-1] == "body"
end

true


text:
using Compat

export AbstractCompositeDataFrame, CompositeDataFrame

abstract AbstractCompositeDataFrame <: AbstractDataFrame

function CompositeDataFrame(columns::Vector{Any},
                            cnames::Vector{Symbol} = gennames(length(columns)))
    # TODO: length checks
    typename = symbol("CompositeDF" * string(gensym()))
    e = :(type $(typename) <: AbstractCompositeDataFrame end)
    e.args[3].args = Any[:($(cnames[i]) :: $(typeof(columns[i]))) for i in 1:length(columns)]
    eval(e)   # create the type
    typ = eval(typename)
    return typ(columns...)
end

CompositeDataFrame(; kwargs...) =
    CompositeDataFrame(Any[ v for (k, v) in kwargs ],
                       Symbol[ k for (k, v) in kwargs ])

# CompositeDataFrame(df::DataFrame) = CompositeDataFrame(df.columns, names(df))

CompositeDataFrame(adf::AbstractDataFrame) =
    CompositeDataFrame(values(adf), names(adf))

CompositeDataFrame(adf::AbstractDataFrame, nms::Vector{Symbol}) =
    CompositeDataFrame(values(adf), nms)

DataFrames.DataFrame(cdf::AbstractCompositeDataFrame) = DataFrame(values(cdf), names(cdf))

#########################################
## basic stuff
#########################################

Base.names{T <: AbstractCompositeDataFrame}(cdf::T) = @compat fieldnames(T)

DataFrames.ncol(cdf::AbstractCompositeDataFrame) = length(names(cdf))
DataFrames.nrow(cdf::AbstractCompositeDataFrame) = ncol(cdf) > 0 ? length(cdf.(1))::Int : 0

Base.values(cdf::AbstractCompositeDataFrame) = Any[ cdf.(i) for i in 1:length(cdf) ]

function Base.hcat(df1::AbstractCompositeDataFrame, df2::AbstractCompositeDataFrame)
    nms = DataFrames.make_unique([names(df1); names(df2)])
    columns = Any[values(df1)..., values(df2)...]
    return CompositeDataFrame(columns, nms)
end
Base.hcat(df1::DataFrame, df2::AbstractCompositeDataFrame) = hcat(df1, DataFrame(df2))
Base.hcat(df1::AbstractCompositeDataFrame, df2::AbstractDataFrame) = hcat(DataFrame(df1), DataFrame(df2))
Base.hcat(df1::AbstractDataFrame, df2::AbstractCompositeDataFrame) = hcat(DataFrame(df1), DataFrame(df2))

DataFrames.index(cdf::AbstractCompositeDataFrame) = DataFrames.Index(names(cdf))

#########################################
## getindex
#########################################

Base.getindex(cdf::AbstractCompositeDataFrame, col_inds::DataFrames.ColumnIndex) = cdf.(col_inds)
Base.getindex{T <: DataFrames.ColumnIndex}(cdf::AbstractCompositeDataFrame, col_inds::AbstractVector{T}) = CompositeDataFrame(Any[ cdf.(col_inds[i]) for i = 1:length(col_inds) ], names(cdf)[col_inds])
Base.getindex(cdf::AbstractCompositeDataFrame, row_inds, col_inds::DataFrames.ColumnIndex) = cdf.(col_inds)[row_inds]
Base.getindex(cdf::AbstractCompositeDataFrame, row_inds, col_inds) =
    CompositeDataFrame(Any[ cdf.(col_inds[i])[row_inds] for i = 1:length(col_inds) ],
                       Symbol[ names(cdf)[i] for i = 1:length(col_inds) ])

function Base.getindex(cdf::AbstractCompositeDataFrame, row_inds, col_inds::UnitRange)
    if col_inds.start == 1 && col_inds.stop == length(cdf)
        return typeof(cdf)([ cdf.(i)[row_inds] for i in 1:length(cdf) ]...)
    else
        return CompositeDataFrame(Any[ cdf.(col_inds[i])[row_inds] for i = 1:length(col_inds) ], names(cdf)[col_inds])
    end
end

#########################################
## LINQ-like operations
#########################################

order(d::AbstractCompositeDataFrame; args...) =
    d[sortperm(DataFrame(args...)), :]

transform(d::AbstractCompositeDataFrame; kwargs...) =
    CompositeDataFrame(Any[values(d)..., [ isa(v, Function) ? v(d) : v for (k,v) in kwargs ]...],
                       Symbol[names(d)..., [ k for (k,v) in kwargs ]...])

select(d::AbstractCompositeDataFrame; kwargs...) =
    CompositeDataFrame(Any[ v for (k,v) in kwargs ],
                       Symbol[ k for (k,v) in kwargs ])


text:
# Verify that all task assignment models succeed in solving all examples, and
# that all optimal solvers obtain the same lower bound cost
let

    for (i, f) in enumerate([
        pctapf_problem_1,
        pctapf_problem_2,
        pctapf_problem_3,
        pctapf_problem_4,
        pctapf_problem_5,
        pctapf_problem_6,
        pctapf_problem_7,
        pctapf_problem_8,
    ])
        for cost_model in [MakeSpan(), SumOfMakeSpans()]
            costs = Float64[]
            for solver in [
                TaskGraphsMILPSolver(AssignmentMILP()),
                TaskGraphsMILPSolver(SparseAdjacencyMILP()),
                TaskGraphsMILPSolver(GreedyAssignment()),
            ]
                # MILP formulations alone
                sched, problem_spec, _, _ = f(
                    ;
                    cost_function = cost_model,
                    verbose = false,
                )
                model = formulate_milp(
                    solver,
                    sched,
                    problem_spec;
                    cost_model = cost_model,
                )
                optimize!(model)
                @test termination_status(model) == MOI.OPTIMAL
                cost = Int(round(value(objective_function(model))))
                update_project_schedule!(
                    solver,
                    model,
                    sched,
                    problem_spec,
                    get_assignment_matrix(model),
                )
                if !isa(model, GreedyAssignment)
                    push!(costs, cost)
                end
                @test validate(sched)
                @test cost != Inf
            end
                # @show costs
            @test all(costs .== costs[1]) # error because AssignmentMILP is adding a job shop constraint that shouldn't be there?? On toy_problem_3
        end
    end
end
# Root nodes
let

    for (f, costs) in [
        (pctapf_problem_4, [2, 2]),
        (pctapf_problem_8, [8, 16]),
    ]
        for solver in [
            TaskGraphsMILPSolver(AssignmentMILP()),
            TaskGraphsMILPSolver(SparseAdjacencyMILP()),
        ]
            for (i, cost_model) in enumerate([MakeSpan(), SumOfMakeSpans()])
                let
                    sched, problem_spec, env_graph, _ = f(; cost_function = cost_model, verbose = false)
                    let
                        model = formulate_milp(
                            solver,
                            sched,
                            problem_spec;
                            cost_model = cost_model,
                        )
                        optimize!(model)
                        optimal = (termination_status(model) == MathOptInterface.OPTIMAL)
                        optimal_TA_cost = Int(round(value(objective_function(model))))
                        @test optimal == true
                        @test optimal_TA_cost == costs[i]
                    end
                end
            end
        end
    end
end
# Station Sharing
let
        for solver in [
            TaskGraphsMILPSolver(AssignmentMILP()),
            TaskGraphsMILPSolver(SparseAdjacencyMILP()),
        ]
        for (i, dt) in enumerate([0, 1, 2])
            for (cost_model, costs) in [
                (MakeSpan(), [4, 5, 6])
                (SumOfMakeSpans(), [7, 8, 9])
            ]
                let
                    sched, problem_spec, env_graph, _ = pctapf_problem_9(
                        ;
                        cost_function = cost_model,
                        verbose = false,
                        Δt_op = 0,
                        Δt_collect = [dt, 0],
                        Δt_deposit = [0, 0],
                    )
                    model = formulate_milp(
                        solver,
                        sched,
                        problem_spec;
                        cost_model = cost_model,
                    )
                    optimize!(model)
                    optimal = (termination_status(model) == MathOptInterface.OPTIMAL)
                    optimal_TA_cost = Int(round(value(objective_function(model))))
                    # @show i, cost_model
                    @test optimal == true
                    @test optimal_TA_cost == costs[i]
                end
            end
        end
    end
end
# Job shop constraints
let
    sched, problem_spec, env_graph, _ = pctapf_problem_9(;
        verbose=false,Δt_op=0,Δt_collect=[0,0],Δt_deposit=[0,0]
        );
    cost_model=SumOfMakeSpans()
    solver=TaskGraphsMILPSolver(AssignmentMILP())
    model = formulate_milp(solver,sched,problem_spec;cost_model=cost_model)
    optimize!(model)
    @test termination_status(model) == MathOptInterface.OPTIMAL
    assignment_matrix = get_assignment_matrix(model);
    update_project_schedule!(solver,model,sched,problem_spec,assignment_matrix)
    n_edges = ne(get_graph(sched))
    add_job_shop_constraints!(model,sched,problem_spec)
    @test ne(get_graph(sched)) > n_edges
end


text:
"The game object holds all metadata pertinent to the game and all its contents."
mutable struct Game
    turn_number::Int
    id::Int
    players::Dict{Int, Player}
    game_map::GameMap
    me::Player
end

"Initiates a game object collecting all start-state instances for the contained items for pre-game."
function init_game()
    constants = load_constants(JSON.parse(readline(stdin)))
    num_players, my_id = parse.(Int, split(readline()))
    players = Dict(i => generate_players() for i=0:num_players-1)
    me = players[my_id]
    game_map = generate_game_map()
    constants, Game(0, my_id, players, game_map, me)
end

"Sends a list of commands to the engine."
send_commands(commands::Vector) = println(join(commands, " "))

"""
Indicate that your bot is ready to play.
# Arguments
- `name::String` : The name of your bot.
"""
start_game(name::String) = send_commands([name])

"Updates the game object's state."
function update_frame!(game::Game)
    game.turn_number = parse(Int, readline())
    for _ in keys(game.players)
        player, num_ships, num_dropoffs, halite = parse.(Int, split(readline()))
        update_player!(game.players[player], num_ships, num_dropoffs, halite)
    end
    update_game_map!(game.game_map)
    for player in values(game.players)
        for ship in get_ships(player)
            mark_unsafe!(game.game_map[ship.position], ship)
        end
        game.game_map[player.shipyard.position].structure = player.shipyard
        for dropoff in get_dropoffs(player)
            game.game_map[dropoff.position].structure = dropoff
        end
    end
    game
end

"Method to send all commands to the game engine, effectively ending your turn."
end_turn(commands::Vector) = send_commands(commands)

max_stars_repo_name: ewirkerman/Halite-III

max_issues_repo_name: ewirkerman/Halite-III

max_forks_repo_name: ewirkerman/Halite-III


text:
using Test
using ChemistryFeaturization
using Serialization

@testset "formation energy example" begin
    include(
        joinpath(@__DIR__, "..", "examples", "1_formation_energy", "formation_energy.jl"),
    )
    @test !ismissing(
        train_formation_energy(
            num_pts = 10,
            num_epochs = 1,
            data_dir = joinpath(@__DIR__, "test_data/"),
            verbose = false,
        ),
    )
end


text:
using PtFEM

data = Dict(
  # Solid(ndim, nst, nxe, nye, nze, nip, direction=:r, finite_element(nod, nodof))
  :struc_el => Solid(3, 6, 1, 3, 2, 8, Hexahedron(20, 3)),
  :properties => [
    100.0 0.3;
     50.0 0.3
     ],
  :etype => [1, 2, 1, 2, 1, 2],
  :x_coords => [0.0, 0.5],
  :y_coords => [0.0, 1.0, 2.0, 3.0],
  :z_coords => [0.0, -1.0, -2.0],
  :support => [
    ( 1, [0 0 1]), ( 2, [1 0 1]), ( 3, [1 0 1]), ( 4, [0 0 1]), ( 5, [1 0 1]), ( 6, [0 0 1]),
    ( 7, [1 0 1]), ( 8, [1 0 1]), ( 9, [0 0 1]), (10, [1 0 1]), (11, [0 0 0]), (12, [0 0 0]),
    (13, [0 0 0]), (14, [0 1 1]), (16, [0 1 1]), (18, [0 0 0]), (19, [0 0 0]), (20, [0 1 1]),
    (23, [0 1 1]), (25, [0 1 1]), (28, [0 1 1]), (30, [0 0 0]), (31, [0 0 0]), (32, [0 0 0]),
    (33, [0 1 1]), (35, [0 1 1]), (37, [0 0 0]), (38, [0 0 0]), (39, [0 1 1]), (42, [0 1 1]),
    (44, [0 1 1]), (47, [0 1 1]), (49, [0 0 0]), (50, [0 0 0]), (51, [0 0 0]), (52, [0 1 1]),
    (54, [0 1 1]), (56, [0 0 0]), (57, [0 0 0]), (58, [0 1 1]), (61, [0 1 1]), (63, [0 1 1]),
    (66, [0 1 1]), (68, [0 0 0]), (69, [0 0 0]), (70, [0 0 0])
    ],
  :loaded_nodes => [
    ( 1, [0.0 0.0  0.0417]), ( 2, [0.0 0.0 -0.1667]), ( 3, [0.0 0.0 0.0417]),
    (14, [0.0 0.0 -0.1667]), (15, [0.0 0.0 -0.1667]), (20, [0.0 0.0 0.0417]),
    (21, [0.0 0.0 -0.1667]), (22, [0.0 0.0  0.0417])
    ],
  :cg_tol => 1.0e-6,
  :cg_limit => 200
)

function test_p56_skyline(data::Dict, N::Int)
  for i in 1:N
    m = p56_skyline(data, true)
  end
end

@time m = p56_skyline(data, true)

#=
Profile.clear()  # in case we have any previous profiling data
@profile test_p56_skyline(data, 20)
using ProfileView
ProfileView.view()
=#

@time m = p56_skyline(data)
println()


text:
using Persa
using Test
using DatasetsCF

using ModelBasedCF

dataset = DatasetsCF.MovieLens()

@testset "Random Model Tests" begin
    model = ModelBasedCF.RandomModel(dataset)
    Persa.train!(model, dataset)

    @test !isnan(model[1,1])
end

@testset "Factorization Matrix Models Tests" begin
    @testset "Baseline Tests" begin
        model = ModelBasedCF.Baseline(dataset)
        Persa.train!(model, dataset, max_epochs = 1)

        @test !isnan(model[1,1])
    end

    @testset "RSVD Tests" begin
        model = ModelBasedCF.RSVD(dataset, 1)
        Persa.train!(model, dataset, max_epochs = 1)

        @test !isnan(model[1,1])
    end

    @testset "IRSVD Tests" begin
        model = ModelBasedCF.IRSVD(dataset, 1)
        Persa.train!(model, dataset, max_epochs = 1)

        @test !isnan(model[1,1])
    end
end


text:
using Vinyl: overdub, @overdub, @hook, @primitive
using Test

mutable struct CountCtx
  count
end

@hook c::CountCtx function sin(x)
  c.count += 1
end

function sincount(f, args...)
  ctx = CountCtx(0)
  overdub(ctx, f, args...)
  return ctx.count
end

struct SinCtx end
@primitive SinCtx sin(x) = -1

@testset "Vinyl" begin

@test @overdub(nothing,sum([1,2,3])) == 6

test(x) = sin(x)+cos(x)
test2(x) = test(x) + test(x+1)

@test sincount(test2, 5.0) == 2

@test overdub(SinCtx(), test, 5) == cos(5)-1

include("continuations.jl")

function unused()
  sin(0.)
  return
end

@test sincount(unused) == 1

@test sincount(() -> sin((0.,)...)) == 1

end


text:
#!/bin/bash
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-$JULIA --color=yes --startup-file=no}"
# below gets the directory name of the script, even if there is a symlink involved
# from https://stackoverflow.com/questions/59895/get-the-source-directory-of-a-bash-script-from-within-the-script-itself
SOURCE="${BASH_SOURCE[0]}"
while [ -h "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink
  DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
  SOURCE="$(readlink "$SOURCE")"
  [[ $SOURCE != /* ]] && SOURCE="$DIR/$SOURCE" # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located
done
DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
export JULIA_PROJECT=$DIR/..
export JULIA_LOAD_PATH=@:@stdlib  # exclude default environment
exec $JULIA_CMD -e 'include(popfirst!(ARGS))' "$SOURCE" "$@"
=#
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "basisfile"
        help = "name of the HDF5 containing basis and creation info"
        required = true
    "--outputfile","-o"
        help = "the output pdf file path, otherwise it will make up a name"

end

delete!(ENV, "PYTHONPATH")
using PyCall
using PyPlot
using HDF5
using Pope
using LinearAlgebra
using Statistics
parsed_args = parse_args(ARGS, s)
if parsed_args["outputfile"]==nothing
    parsed_args["outputfile"] = Pope.outputname(parsed_args["basisfile"],"plots","pdf")
end
display(parsed_args);println()
pdf = pyimport("matplotlib.backends.backend_pdf")

function titleme(basisinfo)
    title("Chan $(basisinfo.channel_number)
    Noise file: $(basename(basisinfo.noise_model_file))
    Pulse file: $(basename(basisinfo.pulse_file))")
end

function plot_model_and_residuals(basisinfo)
    figure(figsize=(10,10))
    subplot(311)
    Nsamp, Nb = size(basisinfo.svdbasis.basis)
    artists=plot(basisinfo.svdbasis.basis)
    xlabel("Sample number")
    titleme(basisinfo)
    legend(artists, basisinfo.singular_values,loc="right",title="Singular values")

    subplot(312)
    for i = 1:Nb
        p = basisinfo.svdbasis.projectors[i,:]
        plot(p ./ norm(p))
    end
    xlabel("Sample number")
    ylabel("Projectors (scaled to norm=1)")

    subplot(313)
    noise_std = √(basisinfo.svdbasis.noise_result.autocorr[1])
    normalized_sorted_residuals = sort(basisinfo.std_residuals)/noise_std
    Ntrain = length(normalized_sorted_residuals)
    a=normalized_sorted_residuals[1]
    b=max(1.5*a, normalized_sorted_residuals[ceil(Int, Ntrain*0.9)])
    plot(normalized_sorted_residuals, range(0, stop=1, length=length(normalized_sorted_residuals)))
    xlim(a,b)
    ylim(0,1)
    xlabel("(Residual std)/(Noise std)")
    ylabel("Frac. training pulses w/ lower residual")
    grid(true)
end

function plot_example_pulses(basisinfo)
    for i = 1:ceil(Int,length(basisinfo.percentiles_of_sample_pulses)/9)
        r = (1:9) .+ (i-1)*9
        figure(figsize=(10,10))
        subplot(311)
        artists=plot(basisinfo.example_pulses[:,r])
        legend(artists,basisinfo.percentiles_of_sample_pulses[r],loc="right",title="Percentile of residual std")
        xlabel("Sample number")
        ylabel("Measured pulses")
        titleme(basisinfo)

        subplot(312)
        model_pulses = basisinfo.svdbasis.basis*basisinfo.svdbasis.projectors*basisinfo.example_pulses[:,r]
        residuals = basisinfo.example_pulses[:,r]-model_pulses
        noise_std = √(basisinfo.svdbasis.noise_result.autocorr[1])
        artists=plot(residuals)
        std_residual_normalized = Float32.(std(residuals, dims=1)'[:]/noise_std)
        std_residual_normalized_from_training = basisinfo.std_residuals_of_example_pulses/noise_std
        legend(artists, std_residual_normalized, loc="right", title="std dev/noise std dev")
        xlabel("Sample number")
        ylabel("Residuals (measured-model)")
        ylim(-7*median(std_residual_normalized)*noise_std,7*median(std_residual_normalized)*noise_std)

        subplot(313)
        artists = plot(model_pulses)
        xlabel("Sample number")
        ylabel("Model pulses")
    end
end

function write_all_plots_to_pdf(pdffile)
    for i in plt.get_fignums()
        figure(i)
        pdffile.savefig(i)
        close(i)
    end
end

function main(parsed_args)
    h5 = h5open(parsed_args["basisfile"],"r")
    pdffile = pdf.PdfPages(parsed_args["outputfile"])
    channelnames = string.(sort(parse.(Int,names(h5))))
    for name in channelnames
        println("plotting channel $name")
        channel = parse(Int,name)
        basisinfo = Pope.hdf5load(Pope.SVDBasisWithCreationInfo,h5["$name"])
        plot_model_and_residuals(basisinfo)
        plot_example_pulses(basisinfo)
        write_all_plots_to_pdf(pdffile)
    end
    close(h5)
    pdffile.close()
    plt.close("all")
end

main(parsed_args)


text:
export FastComparer, fn, reset!

type FastComparer{C<:AbstractComparer} <: AbstractComparer
    c::AbstractComparer
    _results::Dict{Tuple, Nullable{Int}}
end

FastComparer() = FastComparer{Comparer}(Comparer())

FastComparer{C}(c::C) = FastComparer{C}(c, Dict{Tuple, Nullable{Int}}())

FastComparer(fn::Function) = FastComparer{Comparer}(Comparer(fn))

fn(fc::FastComparer) = (a, b) -> compare(fc, a, b) < 0

function compare(fc::FastComparer, a, b)
    key1 = (a, b)
    result = get(fc._results, key1, Nullable{Int}())
    if !isnull(result)
        return get(result)
    end

    key2 = (b, a)
    result = get(fc._results, key2, Nullable{Int}())
    if !isnull(result)
        return -get(result)
    end

    baseFn = fn(fc.c)
    result = 0
    if baseFn(a, b)
        # a < b
        result = -1
    elseif baseFn(b, a)
        # b < a
        result = 1
    end
    fc._results[key1] = result
end

function reset!(fc::FastComparer)
    reset!(fc.c)
    empty!(fc._results)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# Test Fishing Problem
# Switched Nonlinear System
# with Uniform Switching Costs
# without Constraints
function fishing( ; solver=Bazinga.ZEROFPR( tol_optim=1e-4,
                                            max_iter=100,
                                            verbose=true,
                                            freq=1 ),
                        objtol = 1e-3,
                        primaltol = 1e-3 )

    println("Test: Fishing, with switching costs")

    # time interval
    t0 = 0.0
    tf = 12.0

    # control sequence
    uvec = Matrix( [repeat([0.0; 1.0], 4, 1); 0.0]' )

    # cost matrix
    C = [1.0 0.0 -1.0 0.0;
         0.0 1.0 0.0  -1.0]
    Q = Matrix( C' * C)

    # initial state
    x0 = [0.5; 0.7; 1; 1]

    # system dynamics
    function nldyn( x::Vector{Float64}, u::Vector{Float64} )
        n = length(x)
        f = zeros(n)
        f[1] = x[1] - x[1]*x[2] - 0.4*x[1]*u[1]
        f[2] = -x[2] + x[1]*x[2] - 0.2*x[2]*u[1]
        return f
    end

    function nldyn_deriv( x::Vector{Float64}, u::Vector{Float64} )
        df = [1.0-x[2]-0.4*u[1]       -x[1]                   0    0;
              x[2]                     -1+x[1]-0.2*u[1]       0    0;
              0                        0                      0    0;
              0                        0                      0    0]
        return df
    end

    # number of points on the fixed grid
    ngrid = 100

    # switching cost
    swc = 1e-3

    # STO problem with switching costs
    prob = scstoproblem(x0, nldyn, nldyn_deriv, uvec, ngrid=ngrid, t0=t0, tf=tf, Q=Q, swc=swc)

    out = solver( prob, L=1.0 )

    print( out )

    # Test Optimal Solution
    #=tauopt = gettau(m)
    @test isapprox(tauopt[1], 2.4434718158206916, atol=primaltol)
    @test isapprox(tauopt[2], 4.122362522221089, atol=primaltol)
    @test isapprox(tauopt[3], 4.433072844726073, atol=primaltol)
    @test isapprox(tauopt[4], 4.68252002334405,  atol=primaltol)
    @test isapprox(tauopt[5], 5.201667827750571, atol=primaltol)
    @test isapprox(tauopt[6], 5.369908894975762, atol=primaltol)
    @test isapprox(tauopt[7], 6.376845190917198, atol=primaltol)
    @test isapprox(tauopt[8], 6.47160340206027, atol=primaltol)

    # Test Optimal Value
    @test isapprox(getobjval(m), 1.3454355602054182, atol=objtol)=#

    println("Passed")
end


text:
"""
    OptSummary

Summary of an `NLopt` optimization

# Members
* `initial`: a copy of the initial parameter values in the optimization
* `final`: a copy of the final parameter values from the optimization
* `fmin`: the final value of the objective
* `feval`: the number of function evaluations
* `optimizer`: the name of the optimizer used, as a `Symbol`
"""
type OptSummary
    initial::Vector{Float64}
    final::Vector{Float64}
    fmin::Float64
    feval::Int
    optimizer::Symbol
end
function OptSummary(initial::Vector{Float64}, optimizer::Symbol)
    OptSummary(initial, initial, Inf, -1, optimizer)
end

"""
    LinearMixedModel

Linear mixed-effects model representation

# Members
* `formula`: the formula for the model
* `mf`: the model frame, mostly used to get the `terms` component for labelling fixed effects
* `wttrms`: a length `nt` vector of weighted model matrices. The last two elements are `X` and `y`.
* `trms`: a vector of unweighted model matrices.  If `isempty(sqrtwts)` the same object as `wttrms`
* `Λ`: a length `nt - 2` vector of lower triangular matrices
* `sqrtwts`: the `Diagonal` matrix of the square roots of the case weights.  Allowed to be size 0
* `A`: an `nt × nt` symmetric matrix of matrices representing `hcat(Z,X,y)'hcat(Z,X,y)`
* `R`: a `nt × nt` matrix of matrices - the upper Cholesky factor of `Λ'AΛ+I`
* `opt`: an [`OptSummary`](@ref) object
"""
type LinearMixedModel{T <: AbstractFloat} <: MixedModel
    formula::Formula
    mf::ModelFrame
    wttrms::Vector
    trms::Vector
    sqrtwts::Diagonal{T}
    Λ::Vector{LowerTriangular{T, Matrix{T}}}
    A::Matrix{Any}        # symmetric cross-product blocks (upper triangle)
    R::Matrix{Any}        # right Cholesky factor in blocks.
    opt::OptSummary
end

"""
    densify(S::SparseMatrix, threshold=0.3)

Convert sparse `S` to `Diagonal` if `S` is diagonal or to `full(S)` if
the proportion of nonzeros exceeds `threshold`.
"""
function densify(S::SparseMatrixCSC, threshold::Real = 0.3)
    m,n = size(S)
    if m == n && isdiag(S)  # convert diagonal sparse to Diagonal
        return Diagonal(diag(S))
    end
    if nnz(S)/(*(size(S)...)) ≤ threshold # very sparse matrices left as is
        return S
    end
    if isbits(eltype(S))
        return full(S)
    end
    # densify a sparse matrix whose elements are arrays of bitstypes
    nzs = nonzeros(S)
    nz1 = nzs[1]
    T = typeof(nz1)
    if !isa(nz1, Array) || !isbits(eltype(nz1)) # branch not tested
        error("Nonzeros must be a bitstype or an array of same")
    end
    sz1 = size(nz1)
    if any(x->typeof(x) ≠ T || size(x) ≠ sz1, nzs) # branch not tested
        error("Inconsistent dimensions or types in array nonzeros")
    end
    M,N = size(S)
    m,n = size(nz1, 1), size(nz1, 2) # this construction allows for nz1 to be a vector
    res = Array(eltype(nz1), M * m, N * n)
    rv = rowvals(S)
    for j in 1:size(S,2)
        for k in nzrange(S, j)
            copy!(Compat.view(res, (rv[k] - 1) * m + (1 : m), (j - 1) * n + (1 : n)), nzs[k])
        end
    end
    res
end
densify(A::AbstractMatrix, threshold::Real = 0.3) = A

function LinearMixedModel(f, mf, trms, Λ, wts)
    n = size(trms[1], 1)
    T = eltype(trms[end])
    optsum = OptSummary(mapreduce(x -> getθ(x), vcat, Λ), :None)
    sqrtwts = Diagonal([sqrt(x) for x in wts])
    wttrms =  isempty(sqrtwts) ? trms :
        size(sqrtwts, 2) == n ? [sqrtwts * t for t in trms] :
        throw(DimensionMismatch("length(wts) must be 0 or length(y)"))
    nt = length(trms)
    A, R = Array{Any}(nt, nt), Array{Any}(nt, nt)
    for j in 1 : nt, i in 1 : j
        A[i, j] = densify(wttrms[i]'wttrms[j])
        R[i, j] = copy(A[i, j])
    end
    for j in 2 : nt
        if isa(R[j, j], Diagonal) || isa(R[j, j], HBlkDiag)
            for i in 1 : (j - 1)     # check for fill-in
                if !isdiag(A[i, j]'A[i, j])
                    for k in j : nt
                        R[j, k] = full(R[j, k])
                    end
                end
            end
        end
    end
    LinearMixedModel(f, mf, wttrms, trms, sqrtwts, Λ, A, R, optsum)
end

"""
    lmm(f::DataFrames.Formula, fr::DataFrames.DataFrame; weights = [])

Create a [`LinearMixedModel`](@ref) from `f`, which contains both fixed-effects terms
and random effects, and `fr`.

The return value is ready to be `fit!` but has not yet been fit.
"""
function lmm(f::Formula, fr::AbstractDataFrame; weights::Vector = [])
    mf = ModelFrame(f, fr)
    X = ModelMatrix(mf)
    T = eltype(X.m)                                       # process the random-effects terms
    retrms = filter(x -> Meta.isexpr(x, :call) && x.args[1] == :|, mf.terms.terms)
    if isempty(retrms)
        throw(ArgumentError("$f has no random-effects terms"))
    end
    trms = Any[remat(e, mf.df) for e in retrms]
    if length(trms) > 1
        nl = [nlevs(t) for t in trms]
        trms = trms[sortperm(nl; rev = true)]
    end
    Λ = LowerTriangular{T, Matrix{T}}[LT(t) for t in trms]
    push!(trms, X.m)
    push!(trms, reshape(convert(Vector{T}, DataFrames.model_response(mf)), (size(X, 1), 1)))
    LinearMixedModel(f, mf, trms, Λ, convert(Vector{T}, weights))
end

function cfactor!(m::LinearMixedModel)
    A, Λ, R = m.A, m.Λ, m.R
    n = size(A, 1)
    for j in 1 : n, i in 1 : j
        inject!(R[i, j], A[i, j])
    end
    for i in eachindex(m.Λ)
        for j in i : n
            tscale!(Λ[i], R[i, j])
        end
        for ii in 1:i
            tscale!(R[ii,i], Λ[i])
        end
        inflate!(R[i, i])
    end
    cfactor!(R)
    m
end

StatsBase.coef(m::LinearMixedModel) = fixef(m)

"""
    fit!(m::LinearMixedModel[, verbose::Bool=false[, optimizer::Symbol=:LN_BOBYQA]])

Optimize the objective of a `LinearMixedModel`.

A value for `optimizer` should be the name of an `NLopt` derivative-free optimizer
allowing for box constraints.
"""
function StatsBase.fit!(m::LinearMixedModel, verbose::Bool=false, optimizer::Symbol=:LN_BOBYQA)
    th = getθ(m)
    k = length(th)
    opt = NLopt.Opt(optimizer, k)
    NLopt.ftol_rel!(opt, 1e-12)   # relative criterion on deviance
    NLopt.ftol_abs!(opt, 1e-8)    # absolute criterion on deviance
    NLopt.xtol_abs!(opt, 1e-10)   # criterion on parameter value changes
    NLopt.lower_bounds!(opt, lowerbd(m))
    feval = 0
    function obj(x::Vector{Float64}, g::Vector{Float64})
        length(g) == 0 || error("gradient not defined")
        feval += 1
        setθ!(m, x) |> cfactor! |> objective
    end
    if verbose
        function vobj(x::Vector{Float64}, g::Vector{Float64})
            length(g) == 0 || error("gradient not defined")
            feval += 1
            val = setθ!(m, x) |> cfactor! |> objective
            print("f_$feval: $(round(val, 5)), [")
            showcompact(x[1])
            for i in 2:length(x) print(","); showcompact(x[i]) end
            println("]")
            val
        end
        NLopt.min_objective!(opt, vobj)
    else
        NLopt.min_objective!(opt, obj)
    end
    fmin, xmin, ret = NLopt.optimize(opt, th)
    ## very small parameter values often should be set to zero
    xmin1 = copy(xmin)
    modified = false
    for i in eachindex(xmin1)
        if 0. < abs(xmin1[i]) < 0.001
            modified = true
            xmin1[i] = 0.
        end
    end
    if modified  # branch not tested
        ff = setθ!(m, xmin1) |> cfactor! |> objective
        if ff ≤ (fmin + 1.e-5)  # zero components if increase in objective is negligible
            fmin = ff
            copy!(xmin, xmin1)
        else
            setθ!(m, xmin) |> cfactor!
        end
    end
    m.opt = OptSummary(th, xmin, fmin, feval, optimizer)
    if verbose
        println(ret)
    end
    m
end

"""
    lowerbd(m::LinearMixedModel)

Return the vector of lower bounds on the covariance parameter vector `θ`
"""
lowerbd(m::LinearMixedModel) = mapreduce(lowerbd, vcat, m.Λ)

"""
    objective(m::LinearMixedModel)

Return negative twice the log-likelihood of model `m`
"""
objective(m::LinearMixedModel) = logdet(m) + nobs(m) * (1 + log2π + log(varest(m)))

"""
    fixef!{T}(v::Vector{T}, m::LinearMixedModel{T})

Overwrite `v` with the fixed-effects coefficients of model `m`
"""
function fixef!{T}(v::AbstractVector{T}, m::LinearMixedModel{T})
    if !isfit(m)
        throw(ArgumentError("Model m has not been fit"))
    end
    A_ldiv_B!(feR(m), copy!(v, m.R[end - 1, end]))
end

"""
    fixef(m::MixedModel)

Returns the fixed-effects parameter vector estimate.
"""
fixef{T}(m::LinearMixedModel{T}) = fixef!(Array(T, (size(m)[2],)), m)

StatsBase.df(m::LinearMixedModel) = size(m.wttrms[end - 1], 2) + sum(A -> nlower(A), m.Λ) + 1

StatsBase.loglikelihood(m::LinearMixedModel) = -deviance(m)/2

StatsBase.nobs(m::LinearMixedModel) = Int(length(m.trms[end]))

function Base.size(m::LinearMixedModel)
    k = length(m.Λ)
    trms = m.wttrms
    q = sum([size(trms[j], 2) for j in 1:k])
    n, p = size(trms[k + 1])
    n, p, q, k
end

"""
    sdest(m::LinearMixedModel)

Return the estimate of σ, the standard deviation of the per-observation noise.
"""
sdest{T}(m::LinearMixedModel{T}) = T(sqrtpwrss(m)/√nobs(m))

"""
    setθ!{T}(m::LinearMixedModel{T}, v::Vector{T})

Install `v` as the θ parameters in `m`.  Changes `m.Λ` only.
"""
function setθ!{T}(m::LinearMixedModel{T}, v::Vector{T})
    Λ = m.Λ
    offset = 0
    for i in eachindex(Λ)
        λ = Λ[i]
        nti = nlower(λ)
        setθ!(λ, Compat.view(v, offset + (1 : nti)))
        offset += nti
    end
    if length(v) ≠ offset
        throw(DimensionMismatch("length(v) = $(length(v)), should be $offset"))
    end
    m
end

"""
    sqrtpwrss(m::LinearMixedModel)

Return the square root of the penalized, weighted residual sum-of-squares (pwrss).

This value is the contents of the `1 × 1` bottom right block of `m.R`
"""
sqrtpwrss{T}(m::LinearMixedModel{T}) = T(m.R[end,end][1])

"""
    varest(m::LinearMixedModel)

Returns the estimate of σ², the variance of the conditional distribution of Y given B.
"""
varest{T}(m::LinearMixedModel{T}) = T(pwrss(m)/nobs(m))

"""
    pwrss(m::LinearMixedModel)

The penalized residual sum-of-squares.
"""
pwrss(m::LinearMixedModel) = abs2(sqrtpwrss(m))

"""
    chol2cor(L::LowerTriangular)

Return the correlation matrix (symmetric, positive definite with unit diagonal)
corresponding to `L * L'`.
"""
function chol2cor(L::LowerTriangular)
    size(L, 1) == 1 && return ones(1, 1)
    res = L * L'
    d = [inv(sqrt(res[i, i])) for i in 1:size(res, 1)]
    scale!(d, scale!(res, d))
end

Base.cor(m::LinearMixedModel) = map(chol2cor, m.Λ)

function StatsBase.deviance(m::LinearMixedModel)
    isfit(m) || error("Model has not been fit")
    objective(m)
end

"""
    isfit(m::LinearMixedModel)

Return a `Bool` indicating if the model has been fit.
"""
isfit(m::LinearMixedModel) = m.opt.fmin < Inf

"""
    lrt(mods::LinearMixedModel...)

Perform sequential likelihood ratio tests on a sequence of models.

The returned value is a `DataFrame` containing information on the likelihood ratio tests.
"""
function lrt(mods::LinearMixedModel...) # not tested
    if (nm = length(mods)) <= 1
        throw(ArgumentError("at least two models are required for a likelihood ratio test"))
    end
    m1 = mods[1]
    n = nobs(m1)
    for i in 2:nm
        if nobs(mods[i]) != n
            throw(ArgumentError("number of observations must be constant across models"))
        end
    end
    mods = mods[sortperm([df(m)::Int for m in mods])]
    degf = Int[df(m) for m in mods]
    dev = [deviance(m)::Float64 for m in mods]
    csqr = unshift!([(dev[i-1]-dev[i])::Float64 for i in 2:nm],NaN)
    pval = unshift!([ccdf(Chisq(degf[i]-degf[i-1]),csqr[i])::Float64 for i in 2:nm],NaN)
    DataFrame(Df = degf, Deviance = dev, Chisq=csqr,pval=pval)
end

"""
    reweight!{T}(m::LinearMixedModel{T}, wts::Vector{T})

Update `m.sqrtwts` from `wts` and `m.wttrms` from `m.trms`.  Recompute `m.A` and `m.R`.
"""
function reweight!{T}(m::LinearMixedModel{T}, weights::Vector{T})
    A, wttrms, trms, sqrtwts = m.A, m.wttrms, m.trms, m.sqrtwts
       # should be able to use map!(sqrt, sqrtwts.diag, weights) but that allocates storage in v0.4
    d = sqrtwts.diag
    if length(weights) ≠ length(d)
        throw(DimensionMismatch("length(weights) = $(length(weights)), should be $(length(d))"))
    end
    for i in eachindex(d)
        d[i] = sqrt(weights[i])
    end
    for j in eachindex(trms)
        A_mul_B!(sqrtwts, copy!(wttrms[j], trms[j]))
    end
    for j in 1 : size(A, 2), i in 1 : j
        Ac_mul_B!(A[i, j], wttrms[i], wttrms[j])
    end
    cfactor!(m)
end

function Base.show(io::IO, m::LinearMixedModel)
    if !isfit(m)
        warn("Model has not been fit")
        return nothing
    end
    n,p,q,k = size(m)
    println(io, "Linear mixed model fit by maximum likelihood")
    println(io, " ", m.formula)
    oo = objective(m)
    nums = showoff([-oo/ 2, oo, aic(m), bic(m)])
    fieldwd = max(maximum(length(nums)) + 2, 11)
    print(' ')
    for label in ["logLik", "-2 logLik", "AIC", "BIC"]
        print(io, Base.cpad(label, fieldwd))
    end
    println(io)
    print(' ')
    for num in nums
        print(io, lpad(num, fieldwd))
    end
    println(io); println(io)

    show(io,VarCorr(m))

    gl = grplevels(m)
    @printf(io," Number of obs: %d; levels of grouping factors: %d", n, gl[1])
    for l in gl[2:end] @printf(io, ", %d", l) end
    println(io)
    @printf(io,"\n  Fixed-effects parameters:\n")
    show(io,coeftable(m))
end

"""
    VarCorr

An encapsulation of information on the fitted random-effects
variance-covariance matrices.

# Members
* `Λ`: the vector of lower triangular matrices from the `MixedModel`
* `fnms`: a `Vector{ASCIIString}` of grouping factor names
* `cnms`: a `Vector{Vector{ASCIIString}}` of column names
* `s`: the estimate of σ, the standard deviation of the per-observation noise

The main purpose of defining this type is to isolate the logic in the show method.
"""
type VarCorr
    Λ::Vector
    fnms::Vector
    cnms::Vector
    s::Float64
    function VarCorr(Λ::Vector, fnms::Vector, cnms::Vector, s::Number)
        length(fnms) == length(cnms) == length(Λ) || throw(DimensionMismatch(""))
        if isfinite(s) && s < 0
            error("s must be non-negative")
        end
        new(Λ, fnms, cnms, s)
    end
end
function VarCorr(m::LinearMixedModel)
    Λ, trms = m.Λ, m.trms
    VarCorr(Λ, [string(trms[i].fnm) for i in eachindex(Λ)],
        [trms[i].cnms for i in eachindex(Λ)], sdest(m))
end

function Base.show(io::IO, vc::VarCorr)
    fnms = isfinite(vc.s) ? vcat(vc.fnms,"Residual") : vc.fnms
    nmwd = maximum(map(strwidth, fnms)) + 1
    write(io, "Variance components:\n")
    stdm = [rowlengths(λ) for λ in vc.Λ]
    cnms = vcat(vc.cnms...)
    if isfinite(vc.s)
        push!(stdm, [1.])
        stdm *= vc.s
        push!(cnms, "")
    end
    cnmwd = max(6, maximum(map(strwidth, cnms))) + 1
    tt = vcat(stdm...)
    vars = showoff(abs2(tt), :plain)
    stds = showoff(tt, :plain)
    varwd = 1 + max(length("Variance"), maximum(map(strwidth, vars)))
    stdwd = 1 + max(length("Std.Dev."), maximum(map(strwidth, stds)))
    write(io, " "^(2+nmwd))
    write(io, Base.cpad("Column", cnmwd))
    write(io, Base.cpad("Variance", varwd))
    write(io, Base.cpad("Std.Dev.", stdwd))
    any(s -> length(s) > 1, stdm) && write(io,"  Corr.")
    println(io)
    cor = [chol2cor(λ) for λ in vc.Λ]
    ind = 1
    for i in 1:length(fnms)
        stdmi = stdm[i]
        write(io, ' ')
        write(io, rpad(fnms[i], nmwd))
        write(io, rpad(cnms[i], cnmwd))
        write(io, lpad(vars[ind], varwd))
        write(io, lpad(stds[ind], stdwd))
        ind += 1
        println(io)
        for j in 2:length(stdmi)
            write(io, " "^(1 + nmwd))
            write(io, rpad(cnms[ind], cnmwd))
            write(io, lpad(vars[ind], varwd))
            write(io, lpad(stds[ind], stdwd))
            ind += 1
            for k in 1:(j-1)
                @printf(io, "%6.2f", cor[i][j, k])
            end
            println(io)
        end
    end
end