text:
function calc_Vandermonde_inv_formula(x,n)

    detX=BigInt(1)
    for j=1:n
        for i=1:j-1
            detX*=(x[j]-x[i])
        end
    end

    L_1=zeros(BigFloat, n, n)
    for i=1:n
        for j=1:n
            if i>=j
                lij=BigInt(1)
                for k=1:i
                    if k!=j
                        lij/=(x[j]-x[k])
                    end
                end
                L_1[i,j]=lij
            end
        end
    end
    L_1[1,1]=BigInt(1)

    U_1=zeros(BigFloat, n, n)
    for i=1:n
        for j=1:n
            if i==j
                U_1[i,j]=BigInt(1)
            elseif j>i
                if i==1
                    U_1[i,j]=BigInt(0)-U_1[i,j-1]*x[j-1]
                else
                    U_1[i,j]=U_1[i-1,j-1]-U_1[i,j-1]*x[j-1]
                end
            end
        end
    end

    U_1L_1=U_1*L_1

    return U_1L_1,detX

end


text:
#######################################################################################
#                                 Load Packages
#######################################################################################
# set the working directory to the directory in which this file is contained
cd(@__DIR__)
# load the package manager
using Pkg
# activate the project environment
Pkg.activate("../../")
using Revise, AdaptiveDesignOptimization, Random, Distributions
includet("Delay_Discounting.jl")
#######################################################################################
#                               Define Experiment Design
#######################################################################################
Random.seed!(120341)
prior = [Uniform(-5, 5), Uniform(-5, 50)]

model = Model(;prior, loglike)

parm_list = (κ = range(-5, 0, length=50) .|> x->10^x,
   τ = range(0, 5, length=11)[2:end])

# parm_list = (κ = [.1],
#    τ = [.2,.5])

design_list = (
    t_ss = [0.0],
    t_ll = [0.43, 0.714, 1, 2, 3,
        4.3, 6.44, 8.6, 10.8, 12.9,
        17.2, 21.5, 26, 52, 104,
        156, 260, 520],
    r_ss = 12.5:12.5:787.5,
    r_ll = [800.0]
)

# design_list = (
#     t_ss = [0.0],
#     t_ll = [5.0, 10.0],
#     r_ss = [12.0, 20.0],
#     r_ll = [80.0]
# )

data_list = (choice=[true, false],)
#######################################################################################
#                              Simulate Experiment
#######################################################################################
using DataFrames
true_parms = (κ=.12, τ=1.5)
n_trials = 100
optimizer = Optimizer(;design_list, parm_list, data_list, model);
design = optimizer.best_design
df = DataFrame(design=Symbol[], trial=Int[], mean_κ=Float64[], mean_τ=Float64[],
    std_κ=Float64[], std_τ=Float64[])
new_data = [:optimal, 0, mean_post(optimizer)..., std_post(optimizer)...]
push!(df, new_data)

for trial in 1:n_trials
    data = simulate(true_parms..., design...)
    design = update!(optimizer, data)
    new_data = [:optimal, trial, mean_post(optimizer)..., std_post(optimizer)...]
    push!(df, new_data)
end
#######################################################################################
#                              Random Experiment
#######################################################################################
randomizer = Optimizer(;design_list, parm_list, data_list, model, design_type=Randomize);
design = randomizer.best_design
new_data = [:random, 0, mean_post(randomizer)..., std_post(randomizer)...]
push!(df, new_data)

for trial in 1:n_trials
    data = simulate(true_parms..., design...)
    design = update!(randomizer, data)
    new_data = [:random, trial, mean_post(randomizer)..., std_post(randomizer)...]
    push!(df, new_data)
end
#######################################################################################
#                                 Plot Results
#######################################################################################
using StatsPlots
@df df plot(:trial, :mean_κ, xlabel="trial", ylabel="mean κ", group=:design, grid=false)
hline!([true_parms.κ], label="true")

@df df plot(:trial, :mean_τ, xlabel="trial", ylabel="mean τ", group=:design, grid=false)
hline!([true_parms.τ], label="true")

@df df plot(:trial, :std_κ, xlabel="trial", ylabel="σ of κ", grid=false, group=:design, ylims=(0,.3))

@df df plot(:trial, :std_τ, xlabel="trial", ylabel="σ of τ", grid=false, group=:design, ylims=(0,2))

ext: jl


text:
@testset "`Log forwarding`" begin
    mktempdir() do logdir
        vector = rand(10)
        data = Dict("string" => "asdf", "int" => 42, "float64" => 42.0, "float32" => 42.0f0,
                    "vector" => vector, "dict" => Dict("a" => 0, 42 => identity))
        logger = LearnLogger(logdir, "test_run")
        channel = Channel()
        forwarding_task = Lighthouse.forward_logs(channel, logger)
        N = 5
        for _ in 1:N
            for (field, value) in data
                put!(channel, field => value)
            end
            put!(channel, "events" => "happens")
        end
        put!(channel, "plotted" => vector)
        close(channel)
        wait(forwarding_task)
        loaded = logger.logged
        for (k, v) in data
            if typeof(v) <: Union{Real,Complex}
                @test length(loaded[k]) == N
                @test first(loaded[k]) == v
            end
        end
    end
end


text:
using RedPitayaDAQServer
using PyPlot

include("config.jl")

# Alternative RedPitayaClusterView(URLs, [1, 3])
rpcv = RedPitayaClusterView(URLs, [true, false, true])

dec = 64
modulus = 4800
base_frequency = 125000000
samples_per_period = div(modulus, dec)
periods_per_frame = 3

decimation(rpcv, dec)
samplesPerPeriod(rpcv, samples_per_period)
periodsPerFrame(rpcv, periods_per_frame)

# Cluster is configured as before, channel indices refer directly to channels in the cluster
modeDAC(rpcv, "STANDARD")
frequencyDAC(rpcv, 1, 1, base_frequency / modulus)
frequencyDAC(rpcv, 3, 1, base_frequency / modulus)
frequencyDAC(rpcv, 5, 1, base_frequency / modulus)

println(" frequency = $(frequencyDAC(rpcv,1,1))")
amplitudeDAC(rpcv, 1, 1, 0.5)
phaseDAC(rpcv, 1, 1, 0.0 ) #
signalTypeDAC(rpcv, 1, "SINE")

amplitudeDAC(rpcv, 3, 1, 0.5)
phaseDAC(rpcv, 3, 1, 0.0 )
signalTypeDAC(rpcv, 3, "SINE")

amplitudeDAC(rpcv, 5, 1, 0.5)
phaseDAC(rpcv, 5, 1, 0.0 )
signalTypeDAC(rpcv, 5, "SINE")

triggerMode(rpcv, "EXTERNAL")
ramWriterMode(rpcv, "TRIGGERED")
masterTrigger(rpcv, false)

startADC(rpcv)
masterTrigger(rpcv, true)

sleep(1.0)

# ClusterView only reads from selected RedPitayas
uCurrentPeriod = readData(rpcv, currentFrame(rpcv), 1)

fig = figure(1)
clf()
subplot(1,2,1)
plot(vec(uCurrentPeriod[:,1,:,:]))
# Channels from result can be mapped to channels in the cluster
PyPlot.title("Cluster channel $(translateViewToClusterChannel(rpcv,1))")
subplot(1,2,2)
plot(vec(uCurrentPeriod[:,3,:,:]))
PyPlot.title("Cluster channel $(translateViewToClusterChannel(rpcv,3))")

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
module TestBasicTimestamps

using Timestamps
using Dates
using Test

# c = convert(TimestampArray, cl)

@testset "Timestamp construction" begin
    @test Timestamp(Date(2000,1,1),12).ts  == Date(2000,1,1)
    @test Timestamp(Date(2000,1,1),12).val == 12
end

# facts("Getindex works on TimeType") do

#     context("getindex on single Date, DateTime") do
#         @fact length(c[(Date(2000,1,1))])    => 0
#         @fact c[(Date(2000,1,3))][1].value   => 111.94
#     end

#     context("getindex on Vector Date, DateTime") do
#         @fact length(c[[(Date(2000,1,1)), Date(2000,1,3)]])  => 1
#         @fact c[[(Date(2000,1,1)), Date(2000,1,3)]][1].value => 111.94
#     end

#     context("getindex on StepRange Date, DateTime") do
#         @fact length(c[[(Date(2000,1,1)):Date(2000,1,3)]])  => 1
#         @fact c[[(Date(2000,1,1)):Date(2000,1,3)]][1].value => 111.94
#     end
# end

end # module


text:
################################################################################
#                              DocStringExtension                              #
################################################################################

############################################################
#                        Attributes                        #
############################################################

struct DocThemer <: DocStringExtensions.Abbreviation end

const ATTRIBUTES = DocThemer()

function DocStringExtensions.format(::DocThemer, buf, doc)
	local binding = doc.data[:binding] |> Docs.resolve
	help_attributes(buf, binding; extended=true)
end

############################################################
#                        Instances                         #
############################################################

struct DocInstances <: DocStringExtensions.Abbreviation end

const INSTANCES = DocInstances()

function DocStringExtensions.format(::DocInstances, buf, doc)
	local binding = doc.data[:binding] |> Docs.resolve

	# @assert binding isa Enum "Binding $binding must be an `Enum`!"

	insts = instances(binding) # get the instances of the enum

	# initialize a vector of rows for the table
	rows = Vector{Vector{String}}(undef, length(insts) + 1)

	rows[1] = ["Instance", "Value"] # set the header

	# iterate through the instances and create a row for each
	for (i, inst) in enumerate(insts)
		rows[i + 1] = ["`$(inst)`", "`$(Int(inst))`"]
	end

	# print the Markdown table into the buffer
	show(buf, Markdown.MD(Markdown.Table(rows, [:l, :l])))
end

# """
# A lolly instead of a lol!
#
# $(INSTANCES)
# """
# @enum Lolly Pop Bang Snap Crackle Jerk
#
# @doc Lolly
# A lolly instead of a lol!
#
# Instance Value
# –––––––– –––––
# Pop      0
# Bang     1
# Snap     2
# Crackle  3
# Jerk     4


text:
# common facilities

typealias DenseVector{T} DenseArray{T,1}
typealias DenseMatrix{T} DenseArray{T,2}
typealias DenseVecOrMat{T} Union(DenseVector{T},DenseMatrix{T})

typealias NumericArray{T<:Number,N} DenseArray{T,N}
typealias DenseArrOrNum{T<:Number} Union(DenseArray{T}, T)
typealias ArrayOrNum{T<:Number} Union(Array{T}, T)

typealias NumericVector{T<:Number} DenseVector{T}
typealias NumericMatrix{T<:Number} DenseMatrix{T}

typealias ContiguousArrOrNum{T<:Number} Union(ContiguousArray{T}, T)
typealias ContiguousNumericArray{T<:Number} ContiguousArray{T}
typealias ContiguousRealArray{T<:Real} ContiguousArray{T}

typealias DimSpec Union(Int,Dims,Vector{Int})

parent(a::Number) = a
stride(a::Number, d::Int) = 1
ellipview(a::Number, i::Int) = a

getvalue(a::Number, i::Integer) = a
getvalue(a::DenseArray, i::Integer) = a[i]

gt_or_nan(s::Number, x) = (s > x)
lt_or_nan(s::Number, x) = (s < x)

gt_or_nan(s::FloatingPoint, x) = (s > x || s != s)
lt_or_nan(s::FloatingPoint, x) = (s < x || s != s)


text:
export TRARC

struct TRARCWorkspace{T,S,Hess}
    xt::S
    xtnext::S
    d::S
    ∇f::S
    ∇fnext::S
    Hstruct::Hess
    function TRARCWorkspace(::Type{T}, ::Type{S}, ::Type{Hess}, n) where {T,S,Hess}
        return new{T,S,Hess}(
            S(undef, n), # xt
            S(undef, n), # xtnext
            S(undef, n), # d
            S(undef, n), # ∇f
            S(undef, n), # ∇fnext
            Hess(T, S, n),
        )
    end
end

function TRARC(
    nlp_stop::NLPStopping{Pb,M,SRC,NLPAtX{T,S},MStp,LoS};
    TR::TrustRegion = TrustRegion(T(10.0)),
    hess_type::Type{Hess} = HessDense,
    pdata_type::Type{ParamData} = PDataLDLt,
    solve_model::Function = solve_modelTRDiag,
    robust::Bool = true,
    verbose::Bool = false,
    kwargs...,
) where {Pb,M,SRC,MStp,LoS,S,T,Hess,ParamData}
    nlp, nlp_at_x = nlp_stop.pb, nlp_stop.current_state

    PData = ParamData(S, T, nlp.meta.nvar; kwargs...)
    workspace = TRARCWorkspace(T, S, Hess, nlp.meta.nvar)
    xt, xtnext, d, ∇f, ∇fnext =
        workspace.xt, workspace.xtnext, workspace.d, workspace.∇f, workspace.∇fnext

    α = TR.α₀  # initial Trust Region size
    max_unsuccinarow = TR.max_unsuccinarow
    acceptance_threshold = TR.acceptance_threshold
    increase_threshold = TR.increase_threshold
    reduce_threshold = TR.reduce_threshold

    xt .= nlp_at_x.x
    ft, ∇f = objgrad!(nlp, xt, ∇f)
    ftnext = ft
    norm_∇f = norm(∇f)
    nlp_stop.meta.optimality0 = norm_∇f

    OK = update_and_start!(nlp_stop, x = xt, fx = ft, gx = ∇f)
    !OK && Stopping.update!(
        nlp_at_x,
        Hx = hessian!(workspace.Hstruct, nlp, xt),
        convert = true,
    )

    iter = 0 # counter different than stop count
    succ, unsucc, verysucc, unsuccinarow = 0, 0, 0, 0

    verbose && @info log_header(
        [:iter, :f, :nrm_g, :λ, :status, :α, :nrm_dtr, :f_p_dTR, :ΔqN],
        [Int64, T, T, T, String, T, T, T],
    )
    verbose && @info log_row(Any[iter, ft, norm_∇f, 0.0, "First iteration", α])

    while !OK
        max_hprod = nlp_stop.meta.max_cntrs[:neval_hprod]
        PData = preprocess(PData, nlp_at_x.Hx, ∇f, neval_hprod(nlp), max_hprod)

        if ~PData.OK
            @warn("Something wrong with PData")
            return nlp_stop
        end

        success = false
        while !success & (unsuccinarow < max_unsuccinarow)
            d, λ = solve_model(nlp_at_x.Hx, ∇f, nlp_stop, PData, α) # Est-ce que d et λ ne sont pas dans PData ?

            Δq = -(∇f + 0.5 * nlp_at_x.Hx * d) ⋅ d

            if Δq < 0.0
                println("*******   Ascent direction in SolveModel: Δq = $Δq")
                println(
                    "  g⋅d = $(∇f⋅d), 0.5 d'Hd = $(0.5*(nlp_at_x.Hx*d)⋅d)  α = $α  λ = $λ",
                )
                # cond issue with H?
                return nlp_stop
            end
            slope = ∇f ⋅ d
            xtnext .= xt .+ d
            ftnext = obj(nlp, xtnext)
            Δf = ft - ftnext

            iter += 1

            r, good_grad, ∇fnext =
                compute_r(nlp, ft, Δf, Δq, slope, d, xtnext, ∇fnext, robust)

            if r < acceptance_threshold # unsucessful
                verbose && @info log_row(Any[iter, ft, norm_∇f, λ, "U", α, norm(d), Δq])
                unsucc += 1
                unsuccinarow += 1
                α = decrease(PData, α, TR)
            else
                success = true
                unsuccinarow = 0

                xt .= xtnext
                ft = ftnext
                if good_grad
                    ∇f .= ∇fnext
                else
                    grad!(nlp, xt, ∇f)
                end

                verysucc += 1
                if r > increase_threshold # very sucessful
                    α = increase(PData, α, TR)
                    verbose &&
                        @info log_row(Any[iter, ft, stop_norm(∇f), λ, "V", α, norm(d), Δq])
                else # sucessful
                    if r < reduce_threshold
                        α = decrease(PData, α, TR)
                    end
                    verbose &&
                        @info log_row(Any[iter, ft, stop_norm(∇f), λ, "S", α, norm(d), Δq])
                    succ += 1
                end
            end
        end # while !success

        nlp_stop.meta.nb_of_stop = iter
        OK = update_and_stop!(nlp_stop, x = xt, fx = ft, gx = ∇f)
        success && Stopping.update!(nlp_at_x, Hx = hessian!(workspace.Hstruct, nlp, xt))
    end # while !OK

    return nlp_stop
end


text:
# Adapted from:
# SOSDEMO3 --- Bound on Global Extremum
# Section 3.3 of SOSTOOLS User's Manual

@testset "SOSDEMO3 with $(factory.optimizer_constructor)" for factory in sdp_factories
    @polyvar x1 x2

    m = SOSModel(factory)

    @variable m γ

    # Constraint : r(x)*(f(x) - gam) >= 0
    # f(x) is the Goldstein-Price function
    f1 = x1+x2+1
    f2 = 19-14*x1+3*x1^2-14*x2+6*x1*x2+3*x2^2
    f3 = 2*x1-3*x2
    f4 = 18-32*x1+12*x1^2+48*x2-36*x1*x2+27*x2^2

    f = (1+f1^2*f2)*(30+f3^2*f4)

    @constraint m f >= γ

    @objective m Max γ

    JuMP.optimize!(m)

    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT || JuMP.primal_status(m) == MOI.NEARLY_FEASIBLE_POINT

    @test isapprox(JuMP.objective_value(m), 3; rtol=1e-2)
end


text:
visa = r"^(?:4[0-9]{12}(?:[0-9]{3})?)$";
input = "4457418557635128";

ismatch(visa, input)  #=> true
if ismatch(visa, input)
    println("credit card found")
    m = match(visa, input)
    println(m.match) #=> 4457418557635128
    println(m.captures) #=> nothing
    println(m.offset) #=> 1
    println(m.offsets) #=> []
end

m = match(r"(ju|l)(i)?(a)", "Julia")
println(m.match) #=> "lia"
println(m.captures) #=> l - i - a
println(m.offset) #=> 3
println(m.offsets) #=> 3 - 4 - 5

str = "The sky is blue"
reg = r"[\w]{3,}"
# matchall returns an array with a RegexMatche for each match
r = matchall(reg, str)
show(r) #=> ["The","sky","blue"]
# eachmatch returns an iterator over all the matches
iter = eachmatch(reg, str)
for i in iter
    println("\"$(i.match)\" ")
end

max_stars_repo_path: Module 1/Chapter02/Old/regexp.jl

max_issues_repo_path: Module 1/Chapter02/Old/regexp.jl

max_forks_repo_path: Module 1/Chapter02/Old/regexp.jl


text:
using Comonicon
using FromFile
using Test

@testset "compat" begin
    include("compat.jl")
end

@testset "options" begin
    include("options.jl")
end

@testset "ast" begin
    include("ast/ast.jl")
end

@testset "frontend" begin
    include("frontend/cast.jl")
    include("frontend/markdown.jl")
end

@testset "codegen" begin
    include("codegen/codegen.jl")
end

@testset "builder" begin
    include("builder/install.jl")
    include("builder/cli.jl")
end

@testset "scripts" begin
    include("scripts.jl")
end

# the script workaround doesn't work on MacOS
# https://github.com/actions/runner/issues/241
# @testset "tools" begin
#     include("tools.jl")
# end


text:
using Oceananigans

include("utils.jl")

Nx, Ny, Nz = 256, 1, 128
Lx, Ly, Lz = 2000, 1, 1000
Nt, Δt = 1000, 10*86400  # Time step of 10 days.

model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), κ=1e-2)

# Set a +1 g/kg salinity anomaly in a square located at the center of the
# vertical slice with width 0.2Lx and height 0.2Lz.
i1, i2 = floor(Int, 4Nx/10), floor(Int, 6Nx/10)
k1, k2 = floor(Int, 4Nz/10), floor(Int, 6Nz/10)
@. model.tracers.S.data[i1:i2, 1, k1:k2] += 1

# Add a NetCDF output writer that saves NetCDF files to the current directory
# "." with a filename prefix of "diffusion_2d_" every 100 iterations.
nc_writer = NetCDFOutputWriter(dir=".", prefix="isotropic_diffusion_2d_", frequency=100)
push!(model.output_writers, nc_writer)

time_step!(model, Nt, Δt)

make_vertical_slice_movie(model, nc_writer, "S", Nt, Δt, model.eos.S₀)


text:
@testset "Lattice" begin
  @testset "Graph" begin
    # first test the graph stuff
    G = Hecke.Graph{Int, Int}()
    append!(G, 1)
    @test haskey(G.degrees, 1)
    @test haskey(G.new_low_degrees, 1)
    @test G.degrees[1] == 0
    @test 1 in keys(G.edges)

    append!(G, 2)
    @test 2 in keys(G.edges)
    append!(G, (2, 3), 5)
    @test haskey(G.degrees, 2)
    @test haskey(G.new_low_degrees, 2)
    @test haskey(G.new_low_degrees, 3)
    @test G.degrees[2] == 1
    @test G.degrees[3] == 1
    @test 3 in keys(G.edges)
    @test 3 in keys(G.edges[2])
    @test 5 == G.edges[2][3]

    append!(G, (3, 4), 5)
    @test !(3 in keys(G.new_low_degrees))
    @test G.degrees[3] == 2

    append!(G, (5, 4), 5)

    # Now G looks like this:
    #
    #       4
    #      / \
    #     /   \
    # 1   3    5
    #     |
    #     2

    @test_throws ErrorException append!(G, 1)

    # Now delete 3
    delete!(G, 3)
    @test !(3 in keys(G.edges))
    @test !(3 in keys(G.edges[2]))
    @test G.degrees[2] == 0
    @test G.degrees[4] == 1
    @test 4 in keys(G.new_low_degrees)

    # Now add it again and restore G
    append!(G, (2, 3), 5)
    append!(G, (3, 4), 5)

    append!(G, (2, 4), 5)
    append!(G, (6, 5), 5)
    append!(G, (5, 7), 5)

    # Now G looks like this:
    #
    #       4      7
    #      / |\   /
    #     /  | \ /
    # 1   3  |  5
    #     |  /   \
    #     2       6

    # Find shortest path
    b, P = @inferred Hecke.find_shortest(G, 2, 4)
    @test b
    @test P == [4, 2]
    b, P = @inferred Hecke.find_shortest(G, 6, 4)
    @test b
    @test P == [4, 5, 6]
    b, P = @inferred Hecke.find_shortest(G, 1, 1)
    @test b
    @test P == [1]

    b, P = @inferred Hecke.find_shortest(G, 1, 2)
    @test !b

    b, P = @inferred Hecke.find_shortest(G, 3, 5)
    @test !b

    # Find something in the intersection
    b, P1, P2 = @inferred Hecke.find_common(G, 3, 5)
    @test b
    @test P1 == [4, 3]
    @test P2 == [4, 5]

    b, P1, P2 = @inferred Hecke.find_common(G, 2, 3)
    @test b
    @test P1 == [3, 2]
    @test P2 == [3]

    b, P1, P2 = @inferred Hecke.find_common(G, 1, 2)
    @test !b

    b, P1, P2 = @inferred Hecke.find_common(G, 2, 7)
    @test !b
  end

  @testset "Lattice of groups" begin
    # For testing don't use the global lattice
    L = Hecke.GrpAbLatticeCreate()
    G = DiagonalGroup([3,3,4,5])
    H, mH = sub(G, [G[1]], true, L)
    Q, mQ = quo(G, [G[2]], true, L)

    b, M = @inferred Hecke.can_map_into(L, H, G)
    @test b
    @test M == mH.map
    @test_throws ErrorException H[1] + G[1] == 2*G[1]
    @test +(H[1], G[1], L) == 2*G[1]

    b, M = @inferred Hecke.can_map_into(L, H, Q)
    @test b
    @test M == mH.map * mQ.map
    @test_throws ErrorException H[1] + Q[1]
    @test +(H[1], Q[1], L) == 2*Q[1]

    b, M = @inferred Hecke.can_map_into(L, Q, H)
    @test !b

    HH, mHH = sub(G, [G[2]], true, L)
    b, GG, MH, MHH = @inferred Hecke.can_map_into_overstructure(L, H, HH)
    @test b
    @test GG == G
    @test MH == mH.map
    @test MHH == mHH.map
    @test_throws ErrorException H[1] + HH[1]
    @test +(H[1], HH[1], L) == G[1] + G[2]

    HHH, mHHH = sub(H, [H[1]], true, L)
    b, GG, MHHH, MHH = @inferred Hecke.can_map_into_overstructure(L, HHH, HH)
    @test b
    @test GG == G
    @test_throws ErrorException HHH[1] + HH[1]
    @test +(HHH[1], HH[1], L) == G[1] + G[2]

    Q2, mQ2 = quo(G, [G[1]], true, L)
    b, GG, MHHH, MHH = @inferred Hecke.can_map_into_overstructure(L, Q, Q2)
    @test !b
  end
end


text:
function _normal_transition_params(model, scheme::ImplicitEulerMaruyama, t0, x0, x1)
  t1 = t0 + scheme.Δt
  μ = x0 + corrected_drift(model, t1, x1, 1.0) * scheme.Δt
  σ = diffusion(model, t1, x1)
  Σ = scheme.Δt * σ * σ'
  μ, Σ
end


text:
module GaussianMixtureRegressors

# http://cs229.stanford.edu/section/more_on_gaussians.pdf
# http://pypr.sourceforge.net/mog.html

using AutomotiveDrivingModels
using GaussianMixtures
using Distributions
using PyCall

@pyimport sklearn.mixture as mixture

import AutomotiveDrivingModels:
    AbstractVehicleBehavior,
    AbstractVehicleBehaviorPreallocatedData,
    AbstractVehicleBehaviorTrainParams,

    preallocate_learning_data,
    select_action,
    calc_action_loglikelihood,
    train

export
    GMRBehavior,
    GMR_TrainParams,

    select_action,
    calc_action_loglikelihood,
    train

const DEFAULT_MIN_EIGENVALUE = 1e-6

function ensure_posdef!(Σ::Matrix{Float64}, min_eigenvalue::Float64)

    @assert(size(Σ, 1) == size(Σ, 2))

    Λ, V = eig(Σ)
    for (j,λ) in enumerate(Λ)
        Λ[j] = max(DEFAULT_MIN_EIGENVALUE, λ)
    end
    Σ[:] = V*diagm(Λ)/V

    Σ
end

function nearestSPD(A::Matrix{Float64})

    # see http://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd

    # output:
    #  α, β ≥ 0.0 such that
    #         α ≤ δ₂(A) ≤ β ≤ α + 2 max(fα, tol)
    #  and a PSD matrix X such that |A - X|₂ = β

    n = size(A, 1)
    @assert(n == size(A, 2)) # ensure it is square

    I = eye(n)

    # symmetrize A into B
    B = (A+A')./2

    # Compute the symmetric polar factor of B. Call it H.
    # Clearly H is itself SPD.
    U, σ, V = svd(B)
    H = V*diagm(σ)*V'

    # get Ahat in the above formula
    Ahat = (B+H)/2

    # ensure symmetry
    Ahat = (Ahat + Ahat')/2;

    # test that Ahat is in fact PD. if it is not so, then tweak it just a bit.
    worked = false
    iteration_count = 0
    while !worked && iteration_count < 100

        iteration_count += 1

        try
            chol(Ahat)
            worked = true
        catch
            # do nothing
        end

        if !worked
            # Ahat failed the chol test. It must have been just a hair off,
            # due to floating point trash, so it is simplest now just to
            # tweak by adding a tiny multiple of an identity matrix.

            min_eig = minimum(eigvals(Ahat))
            Ahat = Ahat + (-min_eig*iteration_count.^2 + eps(Float32))*I
        end
    end

    Ahat
end

type GMR_TrainParams <: AbstractVehicleBehaviorTrainParams

    targets::ModelTargets
    indicators::Vector{AbstractFeature}

    # params for sci-kit learn GMM fitting
    n_components::Int
    n_gmm_iter::Int # number of iterations taken during GMM fitting
    n_init::Int # number of GMMs to attempt to fit
    verbose::Int # defaults to 0, with 1 prints some stuff, with 2 prints more
    tol::Float64 # tolerance on convergence, if Δlogl < tol, we have converged
    min_covar::Float64 # minimum size of covariance matrix diagonal entry

    # params for julia fitting and training
    max_n_indicators::Int
    use_PCA::Bool
    unlearned_component_weight::Float64 # the weighting given to the unlearned component

    function GMR_TrainParams(;
        targets::ModelTargets = ModelTargets(FUTUREDESIREDANGLE, FUTUREACCELERATION),
        indicators::Vector{AbstractFeature} = AbstractFeature[],

        # params for sci-kit learn GMM fitting
        n_components::Int = 2,
        n_gmm_iter::Int = 100,
        n_init::Int = 3,
        verbose::Int = 0,
        tol::Float64 = 1e-3,
        min_covar::Float64 = 1e-6,

        # params for julia fitting and training
        max_n_indicators::Int = 3,
        use_PCA::Bool = false,
        unlearned_component_weight::Float64 = 0.01, # ∈ [0,1]
        )

        retval = new()

        retval.targets = targets
        retval.indicators = indicators

        retval.n_components = n_components
        retval.n_gmm_iter = n_gmm_iter
        retval.n_init = n_init
        retval.verbose = verbose
        retval.tol = tol
        retval.min_covar = min_covar

        retval.max_n_indicators = max_n_indicators
        retval.use_PCA = use_PCA
        retval.unlearned_component_weight = unlearned_component_weight

        retval
    end
end
function Base.print(io::IO, θ::GMR_TrainParams)
    println(io, "GMR Train Params")
    println(io, "\ttargets: ", θ.targets)
    println(io, "\tindicators: ", map(f->symbol(f), θ.indicators))
    println(io, "\tn_components:     ", θ.n_components)
    println(io, "\tmax_n_indicators: ", θ.max_n_indicators)
    println(io, "\tuse_PCA:   ", θ.use_PCA)
    println(io, "\tmin_covar:        ", θ.min_covar)
end

type GMR_PreallocatedData <: AbstractVehicleBehaviorPreallocatedData

    YX::Matrix{Float64} # [nframes × 2+p]
                        # this includes both targets as 1st two columns (lat, lon)
                        # and then the predictors
    YZ::Matrix{Float64} # [nframes × 2+p]
                        # this uses PCA'd predictors

    extractor::FeatureSubsetExtractor
    preprocess::DataPreprocessor
    preprocess_pca::DataPreprocessor
    preprocess_target::DataPreprocessor

    function GMR_PreallocatedData(dset::ModelTrainingData2, params::GMR_TrainParams)

        retval = new()

        targets = params.targets
        indicators = params.indicators
        trainingframes = dset.dataframe
        nframes = nrow(trainingframes)
        nindicators = length(indicators)

        X = Array(Float64, nindicators, nframes)
        Y = Array(Float64, 2, nframes)

        pull_design_and_target_matrices!(X, Y, trainingframes, targets, indicators)

        ###########################

        preprocess_target = ChainedDataProcessor(2)
        push!(preprocess_target, Y, DataClamper)
        process!(Y, preprocess_target.processors[end]) # process in place
        # push!(preprocess_target, Y, DataScaler) # ensure data is standardized before running PCA
        # process!(Y, preprocess_target.processors[end]) # process in place

        @assert(findfirst(v->isnan(v), Y) == 0)
        @assert(findfirst(v->isinf(v), Y) == 0)

        extractor = FeatureSubsetExtractor(Array(Float64, nindicators), indicators)
        preprocess = ChainedDataProcessor(extractor)
        push!(preprocess, X, DataNaReplacer, extractor.indicators)
        process!(X, preprocess.processors[end]) # process in place

        @assert(findfirst(v->isnan(v), X) == 0)
        @assert(findfirst(v->isinf(v), X) == 0)

        push!(preprocess, X, DataClamper)
        process!(X, preprocess.processors[end]) # process in place

        @assert(findfirst(v->isnan(v), X) == 0)
        @assert(findfirst(v->isinf(v), X) == 0)

        push!(preprocess, X, DataScaler) # ensure data is standardized before running PCA
        process!(X, preprocess.processors[end]) # process in place

        @assert(findfirst(v->isnan(v), X) == 0)
        @assert(findfirst(v->isinf(v), X) == 0)

        retval.YX = vcat(Y, copy(X))
        @assert(findfirst(v->isnan(v), retval.YX) == 0)
        @assert(findfirst(v->isinf(v), retval.YX) == 0)

        # Add a PCA scaler
        preprocess_pca = deepcopy(preprocess)
        n_PCA_features = min(nindicators, 20)
        push!(preprocess_pca, X, DataLinearTransform, n_PCA_features)

        Z = Array(Float64, n_PCA_features, nframes)
        Z = process!(Z, X, preprocess_pca.processors[end])

        retval.YZ = vcat(Y, Z)
        @assert(findfirst(v->isnan(v), retval.YZ) == 0)
        @assert(findfirst(v->isinf(v), retval.YZ) == 0)

        retval.extractor = extractor
        retval.preprocess = preprocess
        retval.preprocess_pca = preprocess_pca
        retval.preprocess_target = preprocess_target

        retval
    end
end
function preallocate_learning_data(
    dset::ModelTrainingData2,
    params::GMR_TrainParams)

    GMR_PreallocatedData(dset, params)
end

######################################################################
#
#
#
######################################################################

type GMR
    converged::Bool

    # μ₁₋₂ = μ₁ + Σ₁₂ * Σ₂₂⁻¹ * (x₂ - μ₂) = A*x₂ + b
    vec_A::Vector{Matrix{Float64}} # [n_components [2×nindicators]]
    vec_b::Vector{Vector{Float64}} # [n_components [2]]

    # pdf(a|p), means μⱼ_ₚ and weights βⱼ(p) are functions of p and must be updated every time, covariance is constant
    # the last component is the extra, unlearned component, and its weight is static
    mixture_Act_given_Obs::MixtureModel{Multivariate,Continuous,MvNormal}

    # pdf(p), all pre-computed. Used to compute βⱼ(p)
    mixture_Obs::MixtureModel{Multivariate,Continuous,MvNormal} # p(obs), all pre-computed, should never be edited

    function GMR(gmm::PyObject, params::GMR_TrainParams)

        retval = new()

        if !gmm[:converged_]
            retval.converged = false
        else

            weights = gmm[:weights_]::Vector{Float64} # [n_components]
            means = gmm[:means_]::Matrix{Float64}     # [n_components, n_features]
            covars = gmm[:covars_]::Array{Float64, 3} # shape depends on covariance_type
            @assert(gmm[:covariance_type] == "full")
            # (n_components, n_features)             if 'spherical',
            # (n_features, n_features)               if 'tied',
            # (n_components, n_features)             if 'diag',
            # (n_components, n_features, n_features) if 'full'

            @assert(0.0 ≤ params.unlearned_component_weight ≤ 1.0)

            n_targets = 2
            n_indicators = size(means, 2) - n_targets
            n_components = length(weights)

            # println("n_targets:    ", n_targets)
            # println("n_indicators: ", n_indicators)
            # println("n_components: ", n_components)

            vec_A = Array(Matrix{Float64}, n_components) # μ₁₋₂ = μ₁ + Σ₁₂ * Σ₂₂⁻¹ * (x₂ - μ₂) = A*x₂ + b
            vec_b = Array(Vector{Float64}, n_components)
            vec_G = Array(MvNormal, n_components)
            vec_H = Array(MvNormal, n_components)

            for i = 1 : n_components

                μₐ = vec(means[i,1:n_targets])
                μₚ = vec(means[i,n_targets+1:end])

                Σ = reshape(covars[i, :, :], (n_targets + n_indicators,n_targets + n_indicators))

                Σₐₐ = Σ[1:n_targets,1:n_targets]
                Σₐₚ = Σ[1:n_targets,n_targets+1:end]
                Σₚₚ = nearestSPD(Σ[n_targets+1:end,n_targets+1:end])
                iΣₚₚ = inv(Σₚₚ)

                A = Σₐₚ * iΣₚₚ
                vec_A[i] = A
                vec_b[i] = vec(μₐ - A*μₚ)
                C = nearestSPD(Σₐₐ - Σₐₚ * iΣₚₚ * ((Σₐₚ)'))

                vec_G[i] = MvNormal(Array(Float64, n_targets), C) # p(action|obs), mean and weighting must be updated with each observation, cov is pre-computed
                vec_H[i] = MvNormal(μₚ, Σₚₚ) # p(obs), all pre-computed, should never be edited
            end

            # add the un-learned component (needed because direct learning will NaN with outliers)
            push!(vec_G, MvNormal([0.0,0.0], [0.01,0.1])) # TODO(tim): tune this
            wᵤ = params.unlearned_component_weight
            Act_given_Obs_prior = push!(fill((1-wᵤ)/n_components, n_components), wᵤ)

            retval.converged = true
            retval.vec_A = vec_A
            retval.vec_b = vec_b
            retval.mixture_Act_given_Obs = MixtureModel(vec_G, Act_given_Obs_prior) # p(action|obs), mean and weighting must be updated with each observation, cov is pre-computed
            retval.mixture_Obs = MixtureModel(vec_H, weights) # p(obs), all pre-computed, should never be edited
        end

        retval
    end
end
function Base.print(model::GMR)
    if !model.converged
        println("GMR: <not converged>")
    else
        println("GMR:")
        for (i, mat) in enumerate(model.vec_A)
            println(i)
            print("\t[")
            for j in 1:size(mat,2)
                @printf("  %10.6f", mat[1,j])
            end
            @printf("] + [  %10.6f]\n", model.vec_b[i][1])
            print("\t[")
            for j in 1:size(mat,2)
                @printf("  %10.6f", mat[2,j])
            end
            @printf("] + [  %10.6f]\n", model.vec_b[i][2])
        end
        println("\tmixture_Obs: ")
        println("\t\tprior: ", model.mixture_Obs.prior)
    end
end

_n_learned_components(gmr::GMR) = length(gmr.vec_A)
function nsuffstats(gmr::GMR)

    dimA = length(gmr.vec_A[1])
    n_learned_components = _n_learned_components(gmr)

    n_learned_components * (2*dimA + 2 # bias
                                   + 3 # covariance in mixture_Act_given_Obs
                                   + div(dimA*dimA,2)) # covariance for mixture_Obs
end
function _calc_action_loglikelihood(gmr::GMR, a::Vector{Float64}, f::Vector{Float64})

    if !gmr.converged
        println("converged: ", gmr.converged)
    end

    mixture_Act_given_Obs = gmr.mixture_Act_given_Obs
    mixture_Obs = gmr.mixture_Obs
    n_learned_components = _n_learned_components(gmr)

    for j in 1 : n_learned_components

        # compute the β value, unweighted
        # βⱼ(f) ∝ wⱼ Nⱼ(μₚ, Σₚ)
        wⱼ = mixture_Obs.prior.p[j]
        Nⱼ = mixture_Obs.components[j]

        # println("a   ", a)
        # println("f   ", f)
        # println("wⱼ   ", wⱼ)
        # println("Nⱼ\n", Nⱼ)

        mixture_Act_given_Obs.prior.p[j] = wⱼ * pdf(Nⱼ, f)

        # compute the conditional mean
        #  μₐ_ₚ = A⋅f + b
        A = gmr.vec_A[j]
        b = gmr.vec_b[j]

        # println("A:\n", A)
        # println("b:\n", b)

        mixture_Act_given_Obs.components[j].μ[:] = A*f + b
    end

    # normalize the β values
    sum_β = sum(mixture_Act_given_Obs.prior.p)
    for i in 1 : n_learned_components
        mixture_Act_given_Obs.prior.p[i] /= sum_β
    end

    if findfirst(v->isinf(v), mixture_Act_given_Obs.prior.p) != 0 || findfirst(v->isnan(v), mixture_Act_given_Obs.prior.p) != 0

        @printf("j  %10s  %10s\n", "wⱼ", "pdf(Nⱼ, f)")
        for j in 1 : n_learned_components
            wⱼ = mixture_Obs.prior.p[j]
            Nⱼ = mixture_Obs.components[j]
            @printf("%2d  %10.6f  %10.6f\n", j, wⱼ, pdf(Nⱼ, f))
        end
        # sleep(100)

        println("INF1: ")
        println("\t a|p prior:  ", mixture_Act_given_Obs.prior.p)
        println("\t sum_β:      ", sum_β)
        println("\t indicators: ", f)
        println("\t targets:    ",    a)

        sleep(0.01)

        return -Inf
    end

    # compute logl using the built-in mixture
    logl = logpdf(mixture_Act_given_Obs, a)

    if isinf(logl) || isnan(logl)
        println(mixture_Act_given_Obs.prior.p)
        println(sum_β)
        println("INF2: obs: ", f, "  action: ",  a)
        sleep(0.01)
    end

    logl
end
function _train_model(gmm::PyObject, YX::Matrix{Float64}, chosen_indicators::Vector{Int}, params::GMR_TrainParams)
    # println("GMM trainging ", chosen_indicators)
    # tic()
    gmm[:fit](YX[:,[1;2;chosen_indicators+2]])
    # toc()
    GMR(gmm, params)
end
function _calc_bic_score(gmr::GMR, YX::Matrix{Float64}, chosen_indicators::Vector{Int})

    bic = -Inf

    if gmr.converged

        a = Array(Float64, 2)
        f = Array(Float64, length(chosen_indicators))

        m = size(YX, 1)
        logl = 0.0
        for i in 1 : m

            a[1] = YX[i,1]
            a[2] = YX[i,2]

            for (j,p) in enumerate(chosen_indicators)
                f[j] = YX[i, 2+p]
            end

            logl += _calc_action_loglikelihood(gmr, a, f)
        end

        # # DEBUG
        # println("\tlogl:     ", logl)
        # println("\t2nd term: ", log(m)*nsuffstats(gmr)/2)
        # println("\ttotal:    ", logl - log(m)*nsuffstats(gmr)/2)

        bic = logl - log(m)*nsuffstats(gmr)/2
    end

    bic
end

######################################################################
#
#
#
######################################################################

type GMRBehavior <: AbstractVehicleBehavior

    #=
    The GMR behavior does Gaussian Mixture Regression

    1 - fit a joint GMM over (lat,lon,p₁,p₁,...,pₙ)

        pdf({a,P}) = pdf(x) = ∑ wⱼ N(μⱼ,Σⱼ)

    2 - convert to a conditional: (lat,lon | p₁,p₁,...,pₙ)

        pdf(a|P) = ∑ βⱼ(p) N(μⱼ_ₚ, Σⱼ_ₚ)

            μⱼ_ₚ = μₐ + Σₐₚ/Σₚₚ (p - μₚ)
            Σⱼ_ₚ = Σₐₐ - (Σₐₚ/Σₚₚ)Σₚₐ
            βⱼ(p) ∝ wⱼ N(μₚ, Σₚₚ)
    =#

    gmr::GMR
    targets::ModelTargets
    extractor::FeatureSubsetExtractor # - note that the extractor and processor are linked internally
    processor::DataPreprocessor
    processor_target::DataPreprocessor

    function GMRBehavior(
        gmr::GMR,
        chosen_indicators::Vector{Int},
        params::GMR_TrainParams,
        pre::GMR_PreallocatedData,
        )

        indicators = pre.extractor.indicators[chosen_indicators]
        preprocess = params.use_PCA > 0 ? pre.preprocess_pca : pre.preprocess

        retval = new()
        retval.gmr = gmr
        retval.targets = params.targets
        retval.extractor = FeatureSubsetExtractor(deepcopy(pre.extractor.indicators))
        retval.processor = deepcopy(preprocess, retval.extractor)
        push!(retval.processor, DataSubset, chosen_indicators)
        retval.processor_target = deepcopy(pre.preprocess_target)
        retval
    end
end
function Base.print(io::IO, GM::GMRBehavior)

    num_components = ncomponents(GM.gmr.mixture_Obs)
    n_indicators = length(GM.processor.z)

    println(io, "Gaussian Mixture Regression Model")
    println(io, "\ttargets: ", GM.targets)
    println(io, "\tncomponents: ", num_components)
    println(io, "\tnindicators: ", n_indicators)
    print(GM.gmr)
end
AutomotiveDrivingModels.get_targets(behavior::GMRBehavior) = behavior.targets
function AutomotiveDrivingModels.get_indicators(behavior::GMRBehavior)
    subset_index = findfirst(p->isa(p, DataSubset), behavior.processor.processors)
    subset = behavior.processor.processors[subset_index]
    behavior.extractor.indicators[subset.indeces]
end
_n_learned_components(GM::GMRBehavior) = length(GM.gmr.vec_A)

######################################################################
#
#
#
######################################################################

const TOLERANCE_ZERO_PROB_WEIGHT = 1e-20

function select_action(
    behavior::GMRBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int
    )

    observe!(behavior.extractor, runlog, sn, colset, frame)
    process!(behavior.processor)

    f = behavior.processor.z

    # condition on the mixture
    vec_A = behavior.gmr.vec_A
    vec_b = behavior.gmr.vec_b
    n_components = _n_learned_components(behavior)
    mixture_Obs = behavior.gmr.mixture_Obs
    mixture_Act_given_Obs = behavior.gmr.mixture_Act_given_Obs

    total_prior_sum = mixture_Act_given_Obs.prior.p[end] # include weight from unlearned component
    for i = 1 : n_components
        mixture_Act_given_Obs.prior.p[i] = mixture_Obs.prior.p[i] * pdf(mixture_Obs.components[i], f)
        total_prior_sum += mixture_Act_given_Obs.prior.p[i]
    end
    if abs(total_prior_sum) < TOLERANCE_ZERO_PROB_WEIGHT
        return (0.0,0.0)
    end

    for i = 1 : n_components
        mixture_Act_given_Obs.prior.p[i] /= total_prior_sum
    end
    @assert(!isnan(mixture_Act_given_Obs.prior.p[1]))

    # pick a component to sample from
    component_index = rand(mixture_Act_given_Obs.prior)
    component = mixture_Act_given_Obs.components[component_index]

    if component_index != n_components+1
        # not the unlearned component
        # modify said component with observation
        component.μ[:] = vec_b[component_index] + vec_A[component_index]*f
    end

    # sample from said component
    Distributions._rand!(component, behavior.processor_target.x)
    process!(behavior.processor_target)

    # return the selected action
    action_lat = behavior.processor_target.z[1]
    action_lon = behavior.processor_target.z[2]

    (action_lat, action_lon)
end

function _set_and_process_action!(behavior::GMRBehavior, action_lat::Float64, action_lon::Float64)
    behavior.processor_target.x[1] = action_lat
    behavior.processor_target.x[2] = action_lon
    process!(behavior.processor_target)
    behavior.processor_target
end
function _calc_action_loglikelihood(behavior::GMRBehavior)

    # NOTE(tim): this assumes that:
    #   - observations have already been pulled and are in behavior.processor.z
    #   - action_lat and action_lon are already processed and in behavior.processor_target.z

    a = behavior.processor_target.z
    f = behavior.processor.z
    _calc_action_loglikelihood(behavior.gmr, a, f)
end

function calc_action_loglikelihood(
    behavior::GMRBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int,
    action_lat::Float64,
    action_lon::Float64,
    )

    #=
    Compute the log-likelihood of the action taken during a single frame
    given the VehicleBehaviorGaussian.
    =#

    observe!(behavior.extractor, runlog, sn, colset, frame)
    process!(behavior.processor)

    _set_and_process_action!(behavior, action_lat, action_lon)

    _calc_action_loglikelihood(behavior)
end
function calc_action_loglikelihood(
    behavior::GMRBehavior,
    features::DataFrame,
    frameind::Integer,
    )

    # println("before: behavior.extractor.x: ", behavior.extractor.x)

    observe!(behavior.extractor, features, frameind)
    process!(behavior.processor)

    # println("after: behavior.extractor.x: ", behavior.extractor.x)
    # println("after: behavior.processor.z: ", behavior.processor.z)
    # println("before: behavior.processor_target.z: ", behavior.processor_target.z)
    # println("before: behavior.processor_target.x: ", behavior.processor_target.x)

    action_lat = features[frameind, symbol(behavior.targets.lat)]::Float64
    action_lon = features[frameind, symbol(behavior.targets.lon)]::Float64
    # println("action_lat: ", action_lat)
    # println("action_lon: ", action_lon)
    _set_and_process_action!(behavior, action_lat, action_lon)

    # println("after: behavior.processor_target.x: ", behavior.processor_target.x)
    # println("after: behavior.processor_target.z: ", behavior.processor_target.z)

    _calc_action_loglikelihood(behavior)
end

######################################################################
#
#
#
######################################################################

function train(
    training_data::ModelTrainingData2, # TODO: do we even need this?
    preallocated_data::GMR_PreallocatedData,
    params::GMR_TrainParams,
    foldset::FoldSet,
    )

    YX = params.use_PCA ?
            copy_matrix_fold(preallocated_data.YZ, foldset) :
            copy_matrix_fold(preallocated_data.YX, foldset)

    @assert(findfirst(v->isnan(v), YX) == 0)
    @assert(findfirst(v->isinf(v), YX) == 0)

    # -----------------------------------------
    # run greedy ascent
    #  - always add the next best parent

    gmm = mixture.GMM(n_components    = params.n_components, # int, optional. Number of mixture components. Defaults to 1.
                      covariance_type = "full", # string, optional. ‘spherical’, ‘tied’, ‘diag’, ‘full’. Defaults to ‘diag’
                      random_state    = 1, # RandomState or an int seed (None by default)
                      min_covar       = params.min_covar, # float, optional
                      n_iter          = params.n_gmm_iter,
                      n_init          = params.n_init,
                      tol             = params.tol,
                      verbose         = params.verbose
                      )

    n_indicators = size(YX, 2) - 2
    max_n_indicators = min(params.max_n_indicators, n_indicators, div(size(YX, 1), 2))
    chosen_indicators = Int[] # start with no parents
    best_model = _train_model(gmm, YX, chosen_indicators, params)
    best_score = _calc_bic_score(best_model, YX, chosen_indicators)

    finished = false
    while !finished && length(chosen_indicators) < max_n_indicators

        finished = true
        test_indicators = [deepcopy(chosen_indicators); 0]

        for i in 1:n_indicators
            if !in(i, chosen_indicators)
                test_indicators[end] = i
                test_model = _train_model(gmm, YX, test_indicators, params)
                test_score = _calc_bic_score(test_model, YX, test_indicators)

                if test_score > best_score
                    best_score = test_score
                    best_model = test_model
                    finished = false
                    chosen_indicators = deepcopy(test_indicators)
                end
            end
        end
    end

    GMRBehavior(best_model, chosen_indicators, params, preallocated_data)
end

end # end module


text:
@kernel function polar_kernel(n::Int, nlines::Int, line_start::Int, scale::Float64,
                     u_curr, v_curr,
                     l_curr, rho,
                     shift_lines::Int, param,
                     _YffR, _YffI,
                     _YftR, _YftI,
                     _YttR, _YttI,
                     _YtfR, _YtfI,
                     frBound, toBound)

    I = @index(Group, Linear)
    J = @index(Local, Linear)
    id_line = I + shift_lines

    x = @localmem Float64 (4,)
    xl = @localmem Float64 (4,)
    xu = @localmem Float64 (4,)

    @inbounds begin
        YffR = _YffR[id_line]; YffI = _YffI[id_line]
        YftR = _YftR[id_line]; YftI = _YftI[id_line]
        YttR = _YttR[id_line]; YttI = _YttI[id_line]
        YtfR = _YtfR[id_line]; YtfI = _YtfI[id_line]

        pij_idx = line_start + 8*(I-1)

        xl[1] = sqrt(frBound[2*(id_line-1)+1])
        xu[1] = sqrt(frBound[2*id_line])
        xl[2] = sqrt(toBound[2*(id_line-1)+1])
        xu[2] = sqrt(toBound[2*id_line])
        xl[3] = -2*pi
        xu[3] = 2*pi
        xl[4] = -2*pi
        xu[4] = 2*pi

        x[1] = min(xu[1], max(xl[1], sqrt(u_curr[pij_idx+4])))
        x[2] = min(xu[2], max(xl[2], sqrt(u_curr[pij_idx+5])))
        x[3] = min(xu[3], max(xl[3], u_curr[pij_idx+6]))
        x[4] = min(xu[4], max(xl[4], u_curr[pij_idx+7]))

        param[1,id_line] = l_curr[pij_idx]
        param[2,id_line] = l_curr[pij_idx+1]
        param[3,id_line] = l_curr[pij_idx+2]
        param[4,id_line] = l_curr[pij_idx+3]
        param[5,id_line] = l_curr[pij_idx+4]
        param[6,id_line] = l_curr[pij_idx+5]
        param[7,id_line] = l_curr[pij_idx+6]
        param[8,id_line] = l_curr[pij_idx+7]
        param[9,id_line] = rho[pij_idx]
        param[10,id_line] = rho[pij_idx+1]
        param[11,id_line] = rho[pij_idx+2]
        param[12,id_line] = rho[pij_idx+3]
        param[13,id_line] = rho[pij_idx+4]
        param[14,id_line] = rho[pij_idx+5]
        param[15,id_line] = rho[pij_idx+6]
        param[16,id_line] = rho[pij_idx+7]
        param[17,id_line] = v_curr[pij_idx]
        param[18,id_line] = v_curr[pij_idx+1]
        param[19,id_line] = v_curr[pij_idx+2]
        param[20,id_line] = v_curr[pij_idx+3]
        param[21,id_line] = v_curr[pij_idx+4]
        param[22,id_line] = v_curr[pij_idx+5]
        param[23,id_line] = v_curr[pij_idx+6]
        param[24,id_line] = v_curr[pij_idx+7]

        @synchronize

        status, minor_iter = tron_kernel(n, shift_lines, 500, 200, 1e-6, scale, true, x, xl, xu,
                                         param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, I, J)

        vi_vj_cos = x[1]*x[2]*cos(x[3] - x[4])
        vi_vj_sin = x[1]*x[2]*sin(x[3] - x[4])

        u_curr[pij_idx] = YffR*x[1]^2 + YftR*vi_vj_cos + YftI*vi_vj_sin
        u_curr[pij_idx+1] = -YffI*x[1]^2 - YftI*vi_vj_cos + YftR*vi_vj_sin
        u_curr[pij_idx+2] = YttR*x[2]^2 + YtfR*vi_vj_cos - YtfI*vi_vj_sin
        u_curr[pij_idx+3] = -YttI*x[2]^2 - YtfI*vi_vj_cos - YtfR*vi_vj_sin
        u_curr[pij_idx+4] = x[1]^2
        u_curr[pij_idx+5] = x[2]^2
        u_curr[pij_idx+6] = x[3]
        u_curr[pij_idx+7] = x[4]
    end
end

function polar_kernel_cpu(n::Int, nline::Int, line_start::Int, scale::Float64,
                          u_curr::AbstractVector{Float64}, v_curr::AbstractVector{Float64},
                          l_curr::AbstractVector{Float64}, rho::AbstractVector{Float64},
                          shift::Int,
                          param::Array{Float64},
                          YffR::Array{Float64}, YffI::Array{Float64},
                          YftR::Array{Float64}, YftI::Array{Float64},
                          YttR::Array{Float64}, YttI::Array{Float64},
                          YtfR::Array{Float64}, YtfI::Array{Float64},
                          frBound::Array{Float64}, toBound::Array{Float64})
    avg_minor_it = 0

    x = zeros(n)
    xl = zeros(n)
    xu = zeros(n)

    xl[3] = -2*pi
    xu[3] = 2*pi
    xl[4] = -2*pi
    xu[4] = 2*pi

    @inbounds for I=1:nline
        pij_idx = line_start + 8*(I-1)
        id_line = shift + I

        xl[1] = sqrt(frBound[2*(id_line-1)+1])
        xu[1] = sqrt(frBound[2*id_line])
        xl[2] = sqrt(toBound[2*(id_line-1)+1])
        xu[2] = sqrt(toBound[2*id_line])

        x[1] = min(xu[1], max(xl[1], sqrt(u_curr[pij_idx+4])))
        x[2] = min(xu[2], max(xl[2], sqrt(u_curr[pij_idx+5])))
        x[3] = min(xu[3], max(xl[3], u_curr[pij_idx+6]))
        x[4] = min(xu[4], max(xl[4], u_curr[pij_idx+7]))

        param[1,id_line] = l_curr[pij_idx]
        param[2,id_line] = l_curr[pij_idx+1]
        param[3,id_line] = l_curr[pij_idx+2]
        param[4,id_line] = l_curr[pij_idx+3]
        param[5,id_line] = l_curr[pij_idx+4]
        param[6,id_line] = l_curr[pij_idx+5]
        param[7,id_line] = l_curr[pij_idx+6]
        param[8,id_line] = l_curr[pij_idx+7]
        param[9,id_line] = rho[pij_idx]
        param[10,id_line] = rho[pij_idx+1]
        param[11,id_line] = rho[pij_idx+2]
        param[12,id_line] = rho[pij_idx+3]
        param[13,id_line] = rho[pij_idx+4]
        param[14,id_line] = rho[pij_idx+5]
        param[15,id_line] = rho[pij_idx+6]
        param[16,id_line] = rho[pij_idx+7]
        param[17,id_line] = v_curr[pij_idx]
        param[18,id_line] = v_curr[pij_idx+1]
        param[19,id_line] = v_curr[pij_idx+2]
        param[20,id_line] = v_curr[pij_idx+3]
        param[21,id_line] = v_curr[pij_idx+4]
        param[22,id_line] = v_curr[pij_idx+5]
        param[23,id_line] = v_curr[pij_idx+6]
        param[24,id_line] = v_curr[pij_idx+7]

        function eval_f_cb(x)
            f = eval_f_polar_kernel_cpu(id_line, scale, x, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            return f
        end

        function eval_g_cb(x, g)
            eval_grad_f_polar_kernel_cpu(id_line, scale, x, g, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            return
        end

        function eval_h_cb(x, mode, rows, cols, _scale, lambda, values)
            eval_h_polar_kernel_cpu(id_line, x, mode, scale, rows, cols, lambda, values,
                                    param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            return
        end

        nele_hess = 10
        tron = ExaTronKernels.createProblem(4, xl, xu, nele_hess, eval_f_cb, eval_g_cb, eval_h_cb;
                                     :tol => 1e-6, :matrix_type => :Dense, :max_minor => 200,
                                     :frtol => 1e-12)

        tron.x .= x
        status = ExaTronKernels.solveProblem(tron)
        x .= tron.x
        avg_minor_it += tron.minor_iter

        cos_ij = cos(x[3] - x[4])
        sin_ij = sin(x[3] - x[4])
        vi_vj_cos = x[1]*x[2]*cos_ij
        vi_vj_sin = x[1]*x[2]*sin_ij

        u_curr[pij_idx] = YffR[id_line]*x[1]^2 + YftR[id_line]*vi_vj_cos + YftI[id_line]*vi_vj_sin
        u_curr[pij_idx+1] = -YffI[id_line]*x[1]^2 - YftI[id_line]*vi_vj_cos + YftR[id_line]*vi_vj_sin
        u_curr[pij_idx+2] = YttR[id_line]*x[2]^2 + YtfR[id_line]*vi_vj_cos - YtfI[id_line]*vi_vj_sin
        u_curr[pij_idx+3] = -YttI[id_line]*x[2]^2 - YtfI[id_line]*vi_vj_cos - YtfR[id_line]*vi_vj_sin
        u_curr[pij_idx+4] = x[1]^2
        u_curr[pij_idx+5] = x[2]^2
        u_curr[pij_idx+6] = x[3]
        u_curr[pij_idx+7] = x[4]
    end

    return 0, avg_minor_it / nline
end


text: # https://people.csail.mit.edu/rivest/pubs/RST01.pdf


text:
function smoothed_aggregation(A::TA,
                        ::Type{Val{bs}}=Val{1},
                        symmetry = HermitianSymmetry(),
                        strength = SymmetricStrength(),
                        aggregate = StandardAggregation(),
                        smooth = JacobiProlongation(4.0/3.0),
                        presmoother = GaussSeidel(),
                        postsmoother = GaussSeidel(),
                        improve_candidates = GaussSeidel(iter=4),
                        max_levels = 10,
                        max_coarse = 10,
                        diagonal_dominance = false,
                        keep = false,
                        coarse_solver = Pinv) where {T,V,bs,TA<:SparseMatrixCSC{T,V}}

    n = size(A, 1)
    # B = kron(ones(n, 1), eye(1))
    B = ones(T,n)

    #=max_levels, max_coarse, strength =
        levelize_strength_or_aggregation(max_levels, max_coarse, strength)
    max_levels, max_coarse, aggregate =
        levelize_strength_or_aggregation(max_levels, max_coarse, aggregate)

    improve_candidates =
        levelize_smooth_or_improve_candidates(improve_candidates, max_levels)=#
    # str = [stength for _ in 1:max_levels - 1]
    # agg = [aggregate for _ in 1:max_levels - 1]
    # sm = [smooth for _ in 1:max_levels]

    levels = Vector{Level{TA, TA, Adjoint{T, TA}}}()
    bsr_flag = false
    w = MultiLevelWorkspace(Val{bs}, eltype(A))
    residual!(w, size(A, 1))

    while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
        A, B, bsr_flag = extend_hierarchy!(levels, strength, aggregate, smooth,
                                improve_candidates, diagonal_dominance,
                                keep, A, B, symmetry, bsr_flag)
        coarse_x!(w, size(A, 1))
        coarse_b!(w, size(A, 1))
        #=if size(A, 1) <= max_coarse
            break
        end=#
        residual!(w, size(A, 1))
    end
    #=A, B = extend_hierarchy!(levels, strength, aggregate, smooth,
                            improve_candidates, diagonal_dominance,
                            keep, A, B, symmetry)=#
    MultiLevel(levels, A, coarse_solver(A), presmoother, postsmoother, w)
end

struct HermitianSymmetry
end

function extend_hierarchy!(levels, strength, aggregate, smooth,
                            improve_candidates, diagonal_dominance, keep,
                            A, B,
                            symmetry, bsr_flag)

    # Calculate strength of connection matrix
    if symmetry isa HermitianSymmetry
        S, _T = strength(A, bsr_flag)
    else
        S, _T = strength(adjoint(A), bsr_flag)
    end

    # Aggregation operator
    AggOp = aggregate(S)
    # b = zeros(eltype(A), size(A, 1))

    # Improve candidates
    b = zeros(size(A,1))
    improve_candidates(A, B, b)
    T, B = fit_candidates(AggOp, B)

    P = smooth(A, T, S, B)
    R = construct_R(symmetry, P)
    push!(levels, Level(A, P, R))

    A = R * A * P

    dropzeros!(A)

    bsr_flag = true

    A, B, bsr_flag
end
construct_R(::HermitianSymmetry, P) = P'

function fit_candidates(AggOp, B, tol = 1e-10)

    A = adjoint(AggOp)
    n_fine, n_coarse = size(A)
    n_col = n_coarse

    R = zeros(eltype(B), n_coarse)
    Qx = zeros(eltype(B), nnz(A))
    # copy!(Qx, B)
    for i = 1:size(Qx, 1)
        Qx[i] = B[i]
    end
    # copy!(A.nzval, B)
    for i = 1:n_col
        for j in nzrange(A,i)
            row = A.rowval[j]
            A.nzval[j] = B[row]
        end
    end
    k = 1
    for i = 1:n_col
        norm_i = norm_col(A, Qx, i)
        threshold_i = tol * norm_i
        if norm_i > threshold_i
            scale = 1 / norm_i
            R[i] = norm_i
        else
            scale = 0
            R[i] = 0
        end
        for j in nzrange(A, i)
            row = A.rowval[j]
            # Qx[row] *= scale
            #@show k
            # Qx[k] *= scale
            # k += 1
            A.nzval[j] *= scale
        end
    end

    # SparseMatrixCSC(size(A)..., A.colptr, A.rowval, Qx), R
    A, R
end
function norm_col(A, Qx, i)
    s = zero(eltype(A))
    for j in nzrange(A, i)
        if A.rowval[j] > length(Qx)
            val = 1
        else
            val = Qx[A.rowval[j]]
        end
        # val = A.nzval[A.rowval[j]]
        s += val*val
    end
    sqrt(s)
end


text:
# Full model:
# Data:  f(x,t): pixel intensity at position x, time t
# Model: f(x,t) = a0 + Σi a(i,t) exp(-b(i,t) |x-x(i,t)|^2)
# a(i,t) = N[a(i,t-1), va]
# b(i,t) = N[b(i,t-1), vb]
# x(i,t) = N[x(i,t-1), vx]

include("explore.jl")
include("model.jl")
include("optim.jl")

function track(y)

end

max_stars_repo_name: JuliaTagBot/NeuronTracking

max_issues_repo_name: JuliaTagBot/NeuronTracking

max_forks_repo_name: JuliaTagBot/NeuronTracking


text:
module StreamStatsTest
using OnlineStats, StatsBase, Base.Test

x = randn(500)
x1 = randn(500)
x2 = randn(501)
xs = hcat(x1, x)

@testset "StreamStats" begin
@testset "HyperLogLog" begin
    o = HyperLogLog(5)
    y = rand(Bool, 10000)
    for yi in y
        fit!(o, yi)
    end
    @test_approx_eq_eps value(o) 2 .5
end
@testset "Bootstrap" begin
    @testset "BernoulliBootstrap" begin
        o = Mean()
        o = BernoulliBootstrap(o, mean, 1000)
        fit!(o, rand(10000))
        cached_state(o)
        mean(o)
        std(o)
        var(o)
        confint(o)
        confint(o, .95, :normal)
        @test_throws Exception confint(o, .95, :fakemethod)
        replicates(o)
    end
    @testset "PoissonBootstrap" begin
        o = Mean()
        o = PoissonBootstrap(o, mean, 1000)
        fit!(o, rand(1000))
        cached_state(o)
        mean(o)
        std(o)
        var(o)
        confint(o)
        replicates(o)
        # vector input
        o = Means(2)
        o = PoissonBootstrap(o, mean, 1000)
        for i in 1:10 fit!(o, rand(2)) end
        cached_state(o)
        @test length(mean(o)) == 2
        std(o)
        var(o)
        #confint(o) not supported right now
        replicates(o)
    end
    @testset "FrozenBootstrap" begin
        o = Mean()
        o = BernoulliBootstrap(o, mean, 1000)

        o2 = Mean()
        o2 = BernoulliBootstrap(o2, mean, 1000)
        fit!(o, randn(100))
        fit!(o2, randn(100) + 3)
        d = o - o2
        mean(d)
        var(d)
        std(d)
        confint(d)
        replicates(o)
    end
end # bootstrap
end # streamstats
end # module


text:
# Basic structures

export SampleData,SampleArray,SampleVector

export duration
export stoi,itos
export channel,times

"""
    stoi(freq::Real,s::{Real,Second,Millisecond})

Convert `s` seconds to data index if data is sampled at `freq`"""
stoi(freq::K,s::T) where {K<:Real,T<:Real} = round(Int,s*freq)
stoi(freq::K,s::Second) where {K<:Real} = stoi(freq,Dates.value(s))
stoi(freq::K,s::Millisecond) where {K<:Real} = stoi(freq,Dates.value(s)/1000)
stoi(freq::K,s::AbstractArray{T,1}) where {K<:Real,T<:Union{Real,Second,Millisecond}} = [stoi(freq,ss) for ss in s]

"""
    itos(freq::Real,i::Integer)

Convert data index `i` to `Millisecond` if data is sampled at `freq`"""
itos(freq::K,i::T) where {K<:Real,T<:Integer} = Millisecond(round(Int,1000*i/freq))

"""
    itos(freq::Real,i::Integer,start::DateTime)
    itos(freq::Real,i::Vector{Integer},start::DateTime)

Convert data index `i` to `DateTime` if data is sampled at `freq` and started at `start`"""
itos(freq::K,i::T,start::DateTime) where {K<:Real,T<:Integer} = itos(freq,i) + start
itos(freq::K,i::AbstractArray{T,1},start::DateTime) where {K<:Real,T<:Integer} = [itos(freq,ii,start) for ii in i]

"""
    times(num::Int,freq::Real,start::DateTime=DateTime(0))

Create a DateTime vector the length `num`, assuming the data was collected at `freq` Hz and started at `start`"""
times(num::Int,freq::K,start::DateTime=DateTime(0)) where K<:Real = [itos(freq,i-1,start) for i in 1:num]

"Abstract type containing data sampled at a certain frequency, along with timing and metadata"
abstract type SampleData{T,K} end

"""
    SampleVector(data::AbstractArray{Number,1}, [times::AbstractArray{{Float,DateTime},1},] freq::Float64; name::String="", code::Symbol=:unknown, start::DateTime=DateTime(0))

Struct of a single `Vector` of sampled data.
* `data`: raw data points
* `times`: corresponding times of raw data points, either as a `Number` or `DateTime`.

   If omitted, it will be calculated as `DateTime` using the given `freq` starting at `start`
* `freq`: data frequency
* `name`: arbitrary name as `String` (default = `""`)
* `code`: arbitrary data type as `Symbol` (default = `:unknown`)

Data can be accessed using []-notation
"""
mutable struct SampleVector{T<:Number,K<:Union{AbstractFloat,DateTime}} <: SampleData{T,K}
    data::AbstractArray{T,1}
    times::AbstractArray{K,1}
    freq::Float64
    name::String
    code::Symbol

    function SampleVector(data::AbstractArray{T,1},times::AbstractArray{K,1},freq::F,name::String,code::Symbol) where {T<:Number,K<:Union{AbstractFloat,DateTime},F<:Real}
        if length(data) != length(times)
            error("""Trying to create a SampleVector with data length $(length(data)) and times length $(length(times))""")
        end
        return new{T,K}(data,times,Float64(freq),name,code)
    end
end

SampleVector(data::AbstractArray{T,1},times::AbstractArray{K,1},freq::F;name::String="",code::Symbol=:unknown) where {T<:Number,K<:Union{AbstractFloat,DateTime},F<:Real} = SampleVector(data,times,Float64(freq),name,code)
SampleVector(data::AbstractArray{T,1},freq::F;name::String="",code::Symbol=:unknown,start::DateTime=DateTime(0)) where {T<:Number,K<:Union{AbstractFloat,DateTime},F<:Real} = SampleVector(data,times(length(data),freq,start),Float64(freq),name,code)

"""
    SampleArray(data::AbstractArray{Number,2}, [times::AbstractArray{{Float,DateTime},1},] freq::Float64; names::Vector{String}=[], code::Symbol=:unknown, start::DateTime=DateTime(0))

Struct of an `Array` frequency data. Makes the assumption `data` was collected at same frequency, synced in time.

* `data`: raw data points; arranged as (time, channels), column-major
* `times`: corresponding times of raw data points, either as a `Number` or `DateTime`.

   If omitted, it will be calculated as `DateTime` using the given `freq` starting at `start`
* `freq`: data frequency
* `names`: `Vector{String}` of names (default = ["","" ...])
* `code`: arbitrary data type as `Symbol` (default = `:unknown`)

Data can be accessed by []-notation, or using the helper function `channel` to return a single `SampleVector`
"""
mutable struct SampleArray{T<:Number,K<:Union{AbstractFloat,DateTime}} <: SampleData{T,K}
    data::AbstractArray{T,2}
    times::AbstractArray{K,1}
    freq::Float64
    names::Vector{String}
    code::Symbol

    function SampleArray(data::AbstractArray{T,2},times::AbstractArray{K,1},freq::F,names::Vector{String},code::Symbol) where {T<:Number,K<:Union{AbstractFloat,DateTime},F<:Real}
        if size(data,1) != length(times)
            error("""Trying to create a SampleArray with data length $(size(data,1)) and times length $(length(times))
                    Data must be arranged in columns""")
        end
        if size(data,2) != length(names)
            error("""Trying to create a SampleArray with data size $(size(data)) and names length $(length(names))
                    Data must be arranged in columns""")
        end
        return new{T,K}(data,times,Float64(freq),names,code)
    end
end

function SampleArray(data::AbstractArray{T,2},times::AbstractArray{K,1},freq::F;names::Vector{String}=[],code::Symbol=:unknown) where {T<:Number,K<:Union{AbstractFloat,DateTime},F<:Real}
    if length(names)==0
        names = ["" for i=1:size(data,2)]
    end
    SampleArray(data,times,freq,names,code)
end

function SampleArray(data::AbstractArray{T,2},freq::F;names::Vector{String}=[],code::Symbol=:unknown,start::DateTime=DateTime(0)) where {T<:Number,K<:Union{AbstractFloat,DateTime},F<:Real}
    SampleArray(data,times(size(data,1),freq,start),freq,names=names,code=code)
end

"""
    SampleArray(v::SampleVector...)

Combine multiple `SampleVector`s into a single `SampleArray`. Is only valid for combining `SampleVector`s that have
the same `code` and `freq`, and that each contain the same number of samples.
"""
function SampleArray(v::SampleVector{T,K}...) where {T,K}
    codes = collect(Set([vv.code for vv in v]))
    freqs = collect(Set([vv.freq for vv in v]))
    length(codes)!=1 && error("Trying to combine SampleVectors containing different codes")
    length(freqs)!=1 && error("Trying to combine SampleVectors with different frequencies")

    lengths = collect(Set([length(vv.data) for vv in v]))
    length(lengths)!=1 && error("Trying to combine SampleVectors of different lengths (stacking together into a SampleArray requires the same length SampleVector's)")

    names = [vv.name for vv in v]
    data = Array{T,2}(undef,(lengths[1],length(names)))
    times = Array{K,2}(undef,(lengths[1],length(names)))

    for i=1:length(names)
        data[:,i] = v[i].data
        times[:,i] = v[i].times
    end

    SampleArray(data,times,freqs[1],names,codes[1])
end

import Base.length,Base.getindex,Base.lastindex,Base.==

itype = Union{Vector{Int},UnitRange,Colon}

length(v::SampleVector) = length(v.data)
length(a::SampleArray) = size(a.data,1)

lastindex(v::SampleData) = length(v)
lastindex(a::SampleArray,d::Int) = size(a.data,d)

function ==(a::SampleVector,b::SampleVector)
    ((a.code != b.code) || (a.freq != b.freq) || (a.name != b.name)) && return false
    return a.data == b.data
end

function ==(a::SampleArray,b::SampleArray)
    ((a.code != b.code) || (a.freq != b.freq) || (a.names != b.names)) && return false
    return a.data == b.data
end

"""length of `SampleData` in `Millisecond`s"""
duration(v::SampleData) = v.times[end] - v.times[1] + Millisecond(round(Int,1000/v.freq))

getindex(v::SampleVector,I::itype) = SampleVector(view(v.data,I),view(v.times,I),v.freq,v.name,v.code)
getindex(a::SampleArray,I::itype) = SampleArray(view(a.data,I,Colon()),view(a.times,I),a.freq,a.names,a.code)

getindex(v::SampleVector,I::StepRange) = SampleVector(view(v.data,I),view(v.times,I),v.freq/I.step,v.name,v.code)
getindex(a::SampleArray,I::StepRange) = SampleArray(view(a.data,I,Colon()),view(a.times,I),a.freq/I.step,a.names,a.code)

getindex(v::SampleVector,I::Int) = v.data[I]
getindex(a::SampleArray,I::Int) = SampleArray(view(a.data,I:I,Colon()),view(a.times,I:I),a.freq,a.names,a.code)

# If passing a Float or DateTime as the index (and the times are corresponding), return the first index at or after given time
getindex(a::SampleData{T,K},f::F) where {T,K<:AbstractFloat,F<:AbstractFloat} = getindex(a,searchsortedfirst(a.times,f))
getindex(a::SampleData{T,K},f::F) where {T,K<:DateTime,F<:DateTime} = getindex(a,searchsortedfirst(a.times,f))
getindex(a::SampleData{T,K},f::F) where {T,K<:DateTime,F<:TimePeriod} = getindex(a,searchsortedfirst(a.times,f+a.times[1]))

"""Extract a specific channel from `SampleArray` and returns a single `SampleVector`"""
channel(a::SampleArray,c::Int,I::itype) = SampleVector(view(a.data,I,c),view(a.times,I),a.freq,a.names[c],a.code)
channel(a::SampleArray,c::Int,I::Int) = a.data[I,c]
channel(a::SampleArray,c::Int) = channel(a,c,Colon())
channel(a::SampleArray,c::String,I::itype) = channel(a,findfirst(isequal(c),a.names),I)
channel(a::SampleArray,c::String,I::Int) = a.data[I,findfirst(isequal(c),a.names)]
channel(a::SampleArray,c::String) = channel(a,c,Colon())

"""Extract a set of channels from `SampleArray` and returns another `SampleArray`"""
channel(a::SampleArray,c::itype,I::itype) = SampleArray(view(a.data,I,c),view(a.times,I),a.freq,a.names[c],a.code)
channel(a::SampleArray,c::itype) = channel(a,c,Colon())
channel(a::SampleArray,c::Vector{String}) = channel(a,[findfirst(isequal(cc),a.names) for cc in c])


text:
@noinline function sinthesin(x)
    sin(sin(x))
end

@noinline function inaloop(y)
    for i = 1:y
	sinthesin(i)
    end
end

type averymutabletype
    a::Int
end

@noinline sinthesin(a::averymutabletype) = sinthesin(a.a)

function testforline(branch)
    a = gcd(4, 25)
    if branch
        b = gcd(7, 21)
    else
        b = gcd(6, 24)
    end
    gcd(a, b)
end

function optional(n = sin(1))
    x = asin(n)
    cos(x)
end


text:
#
# This file containst the functions for single-sets, that is for those cases where
# the list of minimum-distances is between the molecules of a single component.
#
function update_list!(
    i,
    j,
    d2,
    mol_index::F,
    list::List,
) where {F<:Function}
    imol = mol_index(i)
    jmol = mol_index(j)
    if imol != jmol
        d = sqrt(d2)
        if d < list[imol].d
            list[imol] = MinimumDistance(true, i, j, d)
        end
        if d < list[jmol].d
            list[jmol] = MinimumDistance(true, j, i, d)
        end
    end
    return list
end

"""

```
function minimum_distances!(
    mol_index::F,
    x_list::AbstractVector{<:MinimumDistance},
    box, cl;
    parallel = true,
    list_threaded = nothing
) where F <: Function
```

Compute the list of minimum distances given the precomputed cell lists, auxiliary vectors, and
lists of indexes of molecules of each atom. Should be preferred for multiple calls of this function,
with the outer update of the cell lists. The `mol_index` function must return, for each atom index,
the corresponding molecular index. To obtain faster an minimally-allocating calls,
preallocate the `list_threaded` variable.

# Examples

```julia-repl
julia> using CellListMap, StaticArrays

julia> x = [ rand(SVector{2,Float32}) for _ in 1:100 ];

julia> box = Box([1,1],0.2f0);

julia> cl = CellList(x,box);

julia> list = init_list(x, i -> mol_index(i,5));

julia> list_threaded = [ copy(list) for _ in 1:nbatches(cl) ]

julia> minimum_distances!(i -> mol_index(i,5), list, box, cl)
20-element Vector{MinimumDistance{Float32}}:
 MinimumDistance{Float32}(true, 5, 56, 0.033878796f0)
 MinimumDistance{Float32}(true, 9, 69, 0.009123626f0)
 MinimumDistance{Float32}(true, 15, 100, 0.054782398f0)
 ⋮
 MinimumDistance{Float32}(true, 93, 7, 0.03434341f0)
 MinimumDistance{Float32}(true, 96, 29, 0.03825064f0)

```

"""
function minimum_distances!(
    mol_index::F,
    x_list::List,
    box::Box,
    cl::CellList;
    parallel = true,
    list_threaded = nothing
) where {F<:Function}
    reset!(x_list)
    reset!(list_threaded)
    map_pairwise!(
        (x, y, i, j, d2, list) -> update_list!(i, j, d2, mol_index, list),
        x_list,
        box,
        cl;
        parallel = parallel,
        reduce = reduce_list!,
        output_threaded = list_threaded
    )
    return x_list
end

"""

```
function minimum_distances(
    x, n_atoms_per_molecule_x, box;
    parallel = true
)
```

Compute the list of minimum distances given the vector of atomic coordinates `x`,
and the number of atoms of the molecules, the box (of type `CellListMap.Box`). This is the simplest
interface, but it allocates a new list. For iterative computations use `minimum_distances!`
instead.

## Example

```
julia> x = [ rand(SVector{3,Float64}) for _ in 1:100 ];

julia> box =  Box([1,1,1],0.2);

julia> x_list = minimum_distances(x, 5, box);

julia> x_list
20-element Vector{MinimumDistance{Float64}}:
 MinimumDistance{Float64}(true, 1, 84, 0.12675823752297105)
 MinimumDistance{Float64}(true, 9, 20, 0.059338463002696136)
 MinimumDistance{Float64}(true, 13, 53, 0.035903943974538326)
 ⋮
 MinimumDistance{Float64}(true, 94, 34, 0.06721955132313923)
 MinimumDistance{Float64}(true, 99, 86, 0.08788836756712907)

```

Each entry of the output contains the indexes of the atoms involved in the contact, and
the corresponding distance, for each molecule associated to vector `x` of coordinates.

A more general interface, if the molecules do not have the same number of atoms, allows providing a function
that returns the index of each molecule instead of the number of atoms of each molecule:

```
function minimum_distances(
    x, mol_index::F, box;
    parallel = true
) where F <: Function
```

See the user guide for further information.

"""
function minimum_distances(
    x::AbstractVector{<:AbstractVector{T}},
    mol_index::F,
    box::Box{UnitCellType,N,T};
    parallel = true
) where {UnitCellType,N,T,F<:Function}
    cl = CellList(x, box, parallel = parallel)
    x_list = init_list(x, mol_index)
    minimum_distances!(mol_index, x_list, box, cl; parallel = parallel)
    return x_list
end

function minimum_distances(
    x::AbstractVector{<:AbstractVector{T}},
    n_atoms_per_molecule_x::Int,
    box::Box{UnitCellType,N,T};
    parallel = true
) where {UnitCellType,N,T}
    return minimum_distances(
        x,
        i -> mol_index(i, n_atoms_per_molecule_x),
        box;
        parallel = parallel,
    )
end

function minimum_distances(
    x::AbstractVector{AbstractVector{T}},
    n_atoms_per_molecule_x::Int,
    cutoff::T;
    parallel = true
) where {T}
    box = Box(limits(x),cutoff)
    return minimum_distances(
        x,
        i -> mol_index(i, n_atoms_per_molecule_x),
        box;
        parallel = parallel,
    )
end


text:
```
The Accounting Trinity or [CAL]!

[]: important (somehow) - it's objective alright
(): optional read , acts as a comment
me: personal words, to add an overall flavor & meaning to the comment itself - as needed
1.Assets k let's start - here goes nothing

2.Capital :  an abstract entity me: i.e. bank account value- me: has a Direct Relation with either an 'Asset' or a Liability A/C "Account"
3.Assets: a Tangible object, of Value, has an/a owner(s)

below is the explination, non-academic - based on the Personal View -- for your eyes only,ok
# 3.1 thinking about it - I'm a comment after all oh wait maybe I'm not -- who Am ... I ??? - stop talking this is your brain police stop --or what ? - we'll see
#3.2 Ownership:

the actual Owner: for whom the Account belongs to - usualy 1 or more Seperate entities (persons, companies - or mix of the 2 )

1. One Owner: 1 is the lonliest number, if that's you , kudos to you!
2. A couple, partnerships or best Friends (usually Spoils in, or before the end (whatever that means)  you know,  as nothing stands still on its own for good- hate to spoil that on you!) - but sb gotta spill the Beans, right ;) -- sb: somebody
3. 3 ot more: fits the bill for  family-related  members  - [Time: for either long-term or short-term]

3.Liabilities: Debt,  shares US ()(intl.) - how to say in the UK???(if you're authorized to issue that ), given by other people willingly (hopefully they know what they're doing, right) i.e.  Accounts payable A/P  , me: the 'pay me back later, okay?' accounts  - [either long-term or short-term]

together, they make up the following Golden (brown) Equation:

                Assets = Liabilites + Capital
```
# the Accounting Equation , :
## Assets = Liabilites + Capital

include
## ---
export calcCapital, capital,liability
module findtheThird
#nothing, yet - chop chop, gotta go, are you there yet, where you at!
end
 #syntax: incomplete: "module" at /Users/ahmadlutfi/Desktop/src/src/accountingTrinity.jl:3 requires end
#TODO: working input ref ...
capital = parsemissing( input("please insert total $capital"))  # insert a value only
# print("capital:   $capital") # TODO: Ignore Alpha words -- REGEX  - maybe later -- I guessed it
#@assert capital >= 0

#=
function calcCapital(totalAssets,totalLiabilities)
asset = abs(totalAssets);  liability = abs(totalLiabilites) #edit: I don't recall writing those , conciously -wow
=#

return max(asset, liability) - min(asset, liability) #that's another thing I've been thinking about - to add or not the add -- that is the Question
return asset - liability #== capital
end

#---
#capital Calculation

#TODO:calculate Asset  = capital - liability  #Done!
#"""user inserts  totalCapital"""
#global capital(totalCapital) =abs(totalCapital);

global asset(totalAssets) = abs(totalLiabilities)
global liability(totalLiabilities) = abs(totalLiabilities)
#TODO:finish  #Done! #2

function calcAssets(totalCapital, totalLiabilities)
    tota
capital(totalCapital) = capital = parsemissing( input("please insert total $capital")) tryparse(:Number, input("please insert total Capital"))  # capital(totalCapital)
liability(totalLiabilities) = tryparse(:Number, input("please insert total Liabilities")) #abs(totalLiabilities)
#return max(capital, liability) - min(capital, liability)
return capital - liability # Asset
end

function calcCapital(totalAssets, totalLiabilities)
capital(totalCapital) = tryparse(:Number, input("please insert total Capital"))  # capital(totalCapital)
liability(totalLiabilities) = tryparse(:Number, input("please insert total Liabilities")) #abs(totalLiabilities)
#return max(capital, liability) - min(capital, liability)
return capital - liability # Asset
end
end

#--- prompt Area

using Parsers
 function baseline(s,delimiter=' ')
           pieces = split(s, delimiter, keepempty=false)
           map(pieces) do piece
               parse(Float64, piece)
           end
       end

read  = readline()

 #p = Parsers.DELIMITED(" ", ignorerepeated=true);
=#

#=
....
function _sth parse(::Type{T}, s::AbstractString) where T
    y = tryparse(T, s)
    y === nothing && throw(ArgumentError("Cannot parse as $T: \"$s\""))
    return y
end
=#

    function prompt(message="How are you feeling?")
    print(message*"\n")
    line = readline() #read
    # splits on spaces and casts to T
    result = produce(map(x -> parse(T, x), split(line)))
    return result # parsebin(line)
    end

#---
# num = parse(Number, line)
# num =map(x -> parse(Real{Integer},line)

#--- map parse a long dump, long way to go

# num = map(x-> parse(Real{Integer}),line)
#num = custom_parse(T, x::AbstractString) = parse(T, x)  #nope
#num = custom_parse(::Type{T}, x::AbstractString) where #

#{T <: AbstractString} = T(x)
#num = parse(Integer,line)
#my_parse(T, x) = ((T <: AbstractString) && (x isa AbstractString)) ? x : parse(T, x)

#num = map(x -> my_parse(x...), zip((Int64,Float64,String),("123","234","234.7")))
#num = map(x -> my_parse(x), zip((Int64),(line)))
#num = map(x->parse(Integer,x),line)#!null

#if isnumeric(num)#assumes numeric = Number
#    return "$num"*" " *typeof(num)
#else
#end

prompt("How are you feeling?") #it did NOT Prompt - I Swear it's over

#---  ap(x -> parse(T, x)
_out =""
i=0
while i<1 #!eof(f)
    print("anything else?")
    line = readline()#readline(f)
    # splits on spaces and casts to T
    _out = id, s = split(s, "")#produce(map(x -> parse(T, x), split(line)))
i=i+1
end
print(_out)

#--- testing Area: parse stuff - input: file

    function parsestuff(file)
    if (file isa nothing)
        result = ""
    else
        _buffer = read(file)
        len = length(_buffer)
        pos = 1
        opts = Parsers.Options(delim = ' ', wh1 = 0x00)
        while pos <= len
            res = Parsers.xparse(Int, buf, pos, len, opts)
            if Parsers.ok(res.code)
                # parsing succeeded, do stuff with res.val
            end
            pos += res.tlen
        end
    end
end # needed

result = parsestuff(input("how are you doing?"))  # for file or directory - knock, knock you're there ?
print(result)

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using Abaco
using Test

interval=10

function onresult(ts, ne, name, value, inputs)
    @debug "timestamp [$ts]: ne: [$ne] $name = $value"
end

function onresult(timestamp, ne, name, value::Abaco.PValue, inputs)
    @debug "progressive timestamp [$ts]: ne: [$ne] $name = $value"
end

abaco = abaco_init(onresult, interval=interval)

@debug "starting demo with time span: $interval, values timestamp: $ts"

setup_settings(abaco, "state", oncomplete=onresult)

formula(abaco, "total", "sum(region.city.footprint)", "state")

italy = node(abaco, "italy", "state")
veneto = node(abaco, italy, "veneto", "region")
trentino = node(abaco, italy, "trentino", "region")

node(abaco, veneto, "belluno", "city")
node(abaco, veneto, "padova", "city")
node(abaco, trentino, "rovereto", "city")
node(abaco, trentino, "pergine", "city")
node(abaco, trentino, "mori", "city")
node(abaco, trentino, "bondone", "mountain")

nodes = split("region.city.footprint", ".")[1:end-1]
count = Abaco.deep_count(abaco, italy.ne, nodes)
@test count == 5


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractJacobianFunc <: Function end

Abstract type of the jacobian calculators.

Represents the jacobian of the function. The instance can be called at the point `θ`,
a LinearMap representing the jacobian is returned.
"""
abstract type AbstractJacobianFunc <: Function end

struct JacVecProdTag{F, T} end

function _dual_along(f::F, x::AbstractVector{T1}, δ::AbstractVector{T2}) where {F, T1, T2}
    T =  promote_type(T1, T2)
    T_Dual = JacVecProdTag{F, T}
    f(ForwardDiff.Dual{T_Dual}.(x, δ))
end

"""
    FullJacobianFunc(f)

Construct of the Jacobian with ForwardDiff.

When called at point `θ`, the Jacobian matrix is fully instantiated
and stored explicitly in memory.

# Examples
```julia
# forward_model: θ -> Distribution
jacobian_func = FullJacobianFunc(forward_model)
jacobian = jacobian_func(θ)  # LinearMap
jacobian * v  # act on vector
```
"""
struct FullJacobianFunc{F<:Function} <: AbstractJacobianFunc
    f::F
end

function (jf::FullJacobianFunc)(θ::AbstractVector{T}) where T
    λs = jf.f
    jac = ForwardDiff.jacobian(λs, θ)
    LinearMap{T}(jac, isposdef=false, issymmetric=false, ishermitian=false)
end

"""
    FwdRevADJacobianFunc(f)

Construct of the Jacobian with ForwardDiff for direct action and Zygote.pullback for the adjoint

The Jacobian action is computed on the fly, no matrix is stored in memory at any time.

# Examples
```julia
# forward_model: θ -> Distribution
jacobian_func = FwdRevADJacobianFunc(forward_model)
jacobian = jacobian_func(θ)  # LinearMap
jacobian * v  # act on vector as if it was a matrix
```
"""
struct FwdRevADJacobianFunc{F<:Function} <: AbstractJacobianFunc
    f::F
end

function (jf::FwdRevADJacobianFunc)(θ::AbstractVector{T}) where T
    λs = jf.f
    num_λs = size(λs(θ), 1)
    grad_along(f, x, δ) = ForwardDiff.partials.(_dual_along(f, x, δ), 1)
    jvd(δ) = grad_along(λs, θ, δ)
    vjd(δ) = first(Zygote.pullback(λs, θ)[2](δ))
    LinearMap{T}(jvd, vjd, num_λs, size(θ, 1), isposdef=false, issymmetric=false, ishermitian=false)
end

"""
    FwdDerJacobianFunc(f)

Construct of the Jacobian with ForwardDiff for the direct action, and twice applied
ForwardDiff for the adjoint action

Adjoint is implemented by introducing a placeholder parametric vector ``\\vec{t}``:

`` \\frac{d}{d\\vec{t}} \\vec{x} \\cdot (A \\vec{t}) = A^{T} \\vec{x} ``

This allows to implement both, the direct and adjoint actions, using
ForwardDiff, without instantiating full Jacobian at any point.
"""
struct FwdDerJacobianFunc{F<:Function} <: AbstractJacobianFunc
    f::F
end

function (jf::FwdDerJacobianFunc)(θ::AbstractVector{T}) where T
    λs = jf.f
    num_θs = size(θ, 1)
    num_λs = size(λs(θ), 1)
    grad_along(f, x, δ) = ForwardDiff.partials.(_dual_along(f, x, δ), 1)
    jvd(δ) = grad_along(λs, θ, δ)
    vjd(δ) = ForwardDiff.gradient(t -> dot(δ, jvd(t)), zeros(num_θs))
    LinearMap{T}(jvd, vjd, num_λs, num_θs, isposdef=false, issymmetric=false, ishermitian=false)
end

export FullJacobianFunc,
       FwdRevADJacobianFunc,
       FwdDerJacobianFunc


text:
using PyPlot
using JLD2
using Statistics

@load "../../../data/T1_libration_laplace.jld2" eloft1 eloft2 eloft3
cp = ["C0","C1","C2","C3","C4","C5","C6","C7"]
pname=["b","c","d","e","f","g","h"]
nskip0=[1,1,8,8,8,8,8,8]
clf()
figure(figsize=(6.3,6))
i1=1; i20=[210,210,8233,8233,8353,8383,8083]
for j=1:7;
i2 = i20[j]
nskip = nskip0[j]
ecos1 = eloft1[i1:nskip:i2,3,j].*cos.(eloft1[i1:nskip:i2,4,j])
esin1 = eloft1[i1:nskip:i2,3,j].*sin.(eloft1[i1:nskip:i2,4,j])
ecos2 = eloft2[i1:nskip:i2,3,j].*cos.(eloft2[i1:nskip:i2,4,j])
esin2 = eloft2[i1:nskip:i2,3,j].*sin.(eloft2[i1:nskip:i2,4,j])
ecos3 = eloft3[i1:nskip:i2,3,j].*cos.(eloft3[i1:nskip:i2,4,j])
esin3 = eloft3[i1:nskip:i2,3,j].*sin.(eloft3[i1:nskip:i2,4,j])
#plot(ecos1 ,esin1,".")
#plot(ecos2 ,esin2,".")
#plot(ecos3 ,esin3,".")
plot(ecos1 .- mean(ecos1),esin1 .- mean(esin1),color=cp[j],label=pname[j],linewidth=2)
#plot(ecos2 .- mean(ecos2),esin2 .- mean(esin2),".",color=cp[j],alpha=0.5)
#plot(ecos3 .- mean(ecos3),esin3 .- mean(esin3),".",color=cp[j],alpha=0.5)
xlabel(L"e \cos{\omega} - \langle e \cos{\omega}\rangle",fontsize=15)
ylabel(L"e \sin{\omega} - \langle e \sin{\omega}\rangle",fontsize=15)
end
legend(fontsize=10)
tight_layout()
axis([-0.01,0.01,-0.01,0.01])
savefig("../T1_evector_forced.pdf",bbox_inches="tight")
#read(stdin,Char)

clf()
figure(figsize=(6.3,6))
for j=7:-1:1;
i2 = i20[j]
nskip = nskip0[j]
ecos1 = eloft1[i1:nskip:i2,3,j].*cos.(eloft1[i1:nskip:i2,4,j])
esin1 = eloft1[i1:nskip:i2,3,j].*sin.(eloft1[i1:nskip:i2,4,j])
ecos2 = eloft2[i1:nskip:i2,3,j].*cos.(eloft2[i1:nskip:i2,4,j])
esin2 = eloft2[i1:nskip:i2,3,j].*sin.(eloft2[i1:nskip:i2,4,j])
ecos3 = eloft3[i1:nskip:i2,3,j].*cos.(eloft3[i1:nskip:i2,4,j])
esin3 = eloft3[i1:nskip:i2,3,j].*sin.(eloft3[i1:nskip:i2,4,j])
plot(ecos1 ,esin1,color=cp[j],label=pname[j],linewidth=2)
plot(ecos2 ,esin2,color=cp[j],linewidth=2)
plot(ecos3 ,esin3,color=cp[j],linewidth=2)
xlabel(L"e \cos{\omega}",fontsize=15)
ylabel(L"e \sin{\omega}",fontsize=15)
end
legend(fontsize=12)
tight_layout()
axis([-0.02,0.01,-0.015,0.015])
savefig("../T1_evector.pdf",bbox_inches="tight")


text:
module TestPaging

  using GBIF
  using Test

  # No query
  set = occurrences()
  occurrences!(set)
  @test length(set) == 40

  # No query, different limit
  set = occurrences("limit" => 40)
  occurrences!(set)
  @test length(set) == 80
  @test length(set) == length(unique([o.key for o in set]))

  # Query and different limit
  setQ = occurrences(taxon("Iris versicolor", rank=:SPECIES), "limit" => 10)
  occurrences!(setQ)
  @test length(setQ) == 20
  occurrences!(setQ)
  @test length(setQ) == length(unique([o.key for o in setQ]))

end


text:
#test for simpleprover(now in dvcreso.jl)
using Test
include("load_cheaplogic.jl")

@testset "basic prover" begin
 r001,c001=simpleprover(IOBuffer(
"[x].[-P(x)]
[x].[+P(x),-Q(x)]
[x].[+Q(x)]",
), 5,1)

 @test length(contradictionsof(c001))==1

 r002,c002=simpleprover(IOBuffer(
"[X].[-P(X)]
[x].[+P(x),-Q(x)]
[].[+Q(a)]",
), 5,1)

 @test length(contradictionsof(c002))==1

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


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

using Markdown
using InteractiveUtils

# ╔═╡ c1777303-1af2-4be8-8b26-254b70f6263c
using PlutoTest

# ╔═╡ aa7f1ad0-b0e8-11eb-3904-bb242d044d08
md"""
# [Day 15: Science for Hungry People](https://adventofcode.com/2015/day/15)

Today, you set out on the task of perfecting your milk-dunking cookie recipe. All you have to do is find the right balance of ingredients.

Your recipe leaves room for exactly `100` teaspoons of ingredients. You make a list of the remaining ingredients you could use to finish the recipe (your puzzle input) and their properties per teaspoon:

-   `capacity` (how well it helps the cookie absorb milk)
-   `durability` (how well it keeps the cookie intact when full of milk)
-   `flavor` (how tasty it makes the cookie)
-   `texture` (how it improves the feel of the cookie)
-   `calories` (how many calories it adds to the cookie)

You can only measure ingredients in whole-teaspoon amounts accurately, and you have to be accurate so you can reproduce your results in the future. The total score of a cookie can be found by adding up each of the properties (negative totals become `0`) and then multiplying together everything except calories.

For instance, suppose you have these two ingredients:

```
Butterscotch: capacity -1, durability -2, flavor 6, texture 3, calories 8
Cinnamon: capacity 2, durability 3, flavor -2, texture -1, calories 3

```

Then, choosing to use `44` teaspoons of butterscotch and `56` teaspoons of cinnamon (because the amounts of each ingredient must add up to `100`) would result in a cookie with the following properties:

-   A `capacity` of `44*-1 + 56*2 = 68`
-   A `durability` of `44*-2 + 56*3 = 80`
-   A `flavor` of `44*6 + 56*-2 = 152`
-   A `texture` of `44*3 + 56*-1 = 76`

Multiplying these together (`68 * 80 * 152 * 76`, ignoring `calories` for now) results in a total score of `62842880`, which happens to be the best score possible given these ingredients. If any properties had produced a negative total, it would have instead become zero, causing the whole score to multiply to zero.

Given the ingredients in your kitchen and their properties, **what is the total score of the highest-scoring cookie you can make?**
"""

# ╔═╡ 861a2d99-e0a6-44a9-9b18-ca1b25d1eabd
puzzle_input = "Sprinkles: capacity 5, durability -1, flavor 0, texture 0, calories 5
PeanutButter: capacity -1, durability 3, flavor 0, texture 0, calories 1
Frosting: capacity 0, durability -1, flavor 4, texture 0, calories 6
Sugar: capacity -1, durability 0, flavor 0, texture 2, calories 8"

# ╔═╡ ad735981-4fdb-4527-9869-a6a748e2e5dd
example_input = """
Butterscotch: capacity -1, durability -2, flavor 6, texture 3, calories 8
Cinnamon: capacity 2, durability 3, flavor -2, texture -1, calories 3
""" |> strip

# ╔═╡ 34edc8a5-edc7-4a8d-8bbe-c9d354d8561d
function parseline(line)
	m = match(r"(\w+): capacity (-?\d+), durability (-?\d+), flavor (-?\d+), texture (-?\d+), calories (-?\d+)", line)
	m
	(name=m[1], capacity=parse(Int, m[2]), durability=parse(Int, m[3]), flavor=parse(Int, m[4]), texture=parse(Int, m[5]), calories=parse(Int, m[6]))
end

# ╔═╡ 4122ea48-c078-4b72-81ca-0f6163d7b474
example_ingredients = map(parseline, split(example_input, "\n"))

# ╔═╡ cf9d9179-cbd0-4f1d-aee7-5eb7eaed4559
function cookiescore(ingredients, ratios)
	#if sum(ratios) != 100
#		return -1
#end
	scores = [
		sum(map(((idx, ingredient),)-> ingredient.capacity * ratios[idx], enumerate(ingredients))),
		sum(map(((idx, ingredient),)-> ingredient.durability * ratios[idx], enumerate(ingredients))),
		sum(map(((idx, ingredient),)-> ingredient.flavor * ratios[idx], enumerate(ingredients))),
		sum(map(((idx, ingredient),)-> ingredient.texture * ratios[idx], enumerate(ingredients))),
	]
	prod(map(score -> max(0, score), scores))
end

# ╔═╡ 782f3c0c-d95f-402f-8e1a-ff640b885c2b
function cookiecalories(ingredients, ratios)
	calories = 0
	for (idx, ratio) in enumerate(ratios)
		calories += ingredients[idx][6] * ratio
	end
	calories
end

# ╔═╡ b9a51427-a73f-4b53-aac0-35c49846ec3e
@assert cookiescore(example_ingredients, (44, 56)) == 62842880

# ╔═╡ 65abcb4d-b961-4ad7-b78a-1d40aa64c867
cookiescore(example_ingredients, (44, 56))

# ╔═╡ 0e81de2d-45d0-4627-8479-9dc92ce4433e
function bestcookiescore1(ingredients)
	best = 0
	for a in 1:100
		for b in 1:100-a
			for c in 1:100-(a+b)
				d = 100 - a - b -c
				score = cookiescore(ingredients, (a,b,c,d))
				if score > best
					best = score
				end
			end
		end
	end
	best
end

# ╔═╡ 6e394cd1-dbdc-4540-ab6d-39bdf0299a9c
@assert bestcookiescore1(example_ingredients) == 62842880

# ╔═╡ c3d471b9-f0b5-4e29-937e-dcb917cf2a9c
bestcookiescore1(example_ingredients)

# ╔═╡ 59938b46-ff7d-441f-95ba-1a49fcaeefe7
puzzle_ingredients = map(parseline, split(puzzle_input, "\n"))

# ╔═╡ 0d4650b4-1035-43f4-b423-8824cb3436d2
part1 = bestcookiescore1(puzzle_ingredients)

# ╔═╡ fb3a890b-211c-4fa5-a85e-472183509fa8
@test part1 == 13882464

# ╔═╡ 117c9ac8-6f1e-4dc3-9b74-bb2cabe9b171
md"Your puzzle answer was `13882464`."

# ╔═╡ 9e397f78-a72b-406b-905c-dc1f79e852d0
md"""
# Part Two

Your cookie recipe becomes wildly popular! Someone asks if you can make another recipe that has exactly `500` calories per cookie (so they can use it as a meal replacement). Keep the rest of your award-winning process the same (100 teaspoons, same ingredients, same scoring system).

For example, given the ingredients above, if you had instead selected `40` teaspoons of butterscotch and `60` teaspoons of cinnamon (which still adds to `100`), the total calorie count would be `40*8 + 60*3 = 500`. The total score would go down, though: only `57600000`, the best you can do in such trying circumstances.

Given the ingredients in your kitchen and their properties, **what is the total score of the highest-scoring cookie you can make with a calorie total of `500`**?

"""

# ╔═╡ 3ae1cb15-f704-4559-9043-186c2874a035
function bestcookiescore2(ingredients)
	best = 0
	for a in 1:100
		for b in 1:100-a
			for c in 1:100-(a+b)
				d = 100 - a - b -c
				score = cookiescore(ingredients, (a,b,c,d))
				if score > best && cookiecalories(ingredients, (a,b,c,d)) == 500
					best = score
				end
			end
		end
	end
	best
end

# ╔═╡ a10ef2a7-c9d9-41f0-9cfc-6921a390306f
part2 = bestcookiescore2(puzzle_ingredients)

# ╔═╡ 9274f947-ede1-446e-9845-9e776b2bf42c
@test part2 == 11171160

# ╔═╡ da313ad1-afdf-45cd-9601-5e6dec1f92bc
md"Your puzzle answer was `11171160`."

# ╔═╡ Cell order:
# ╠═c1777303-1af2-4be8-8b26-254b70f6263c
# ╟─aa7f1ad0-b0e8-11eb-3904-bb242d044d08
# ╟─861a2d99-e0a6-44a9-9b18-ca1b25d1eabd
# ╠═ad735981-4fdb-4527-9869-a6a748e2e5dd
# ╠═34edc8a5-edc7-4a8d-8bbe-c9d354d8561d
# ╠═4122ea48-c078-4b72-81ca-0f6163d7b474
# ╠═cf9d9179-cbd0-4f1d-aee7-5eb7eaed4559
# ╠═782f3c0c-d95f-402f-8e1a-ff640b885c2b
# ╠═b9a51427-a73f-4b53-aac0-35c49846ec3e
# ╠═65abcb4d-b961-4ad7-b78a-1d40aa64c867
# ╠═0e81de2d-45d0-4627-8479-9dc92ce4433e
# ╠═6e394cd1-dbdc-4540-ab6d-39bdf0299a9c
# ╠═c3d471b9-f0b5-4e29-937e-dcb917cf2a9c
# ╠═59938b46-ff7d-441f-95ba-1a49fcaeefe7
# ╠═0d4650b4-1035-43f4-b423-8824cb3436d2
# ╠═fb3a890b-211c-4fa5-a85e-472183509fa8
# ╟─117c9ac8-6f1e-4dc3-9b74-bb2cabe9b171
# ╟─9e397f78-a72b-406b-905c-dc1f79e852d0
# ╠═3ae1cb15-f704-4559-9043-186c2874a035
# ╠═a10ef2a7-c9d9-41f0-9cfc-6921a390306f
# ╠═9274f947-ede1-446e-9845-9e776b2bf42c
# ╟─da313ad1-afdf-45cd-9601-5e6dec1f92bc


text:
using DiffEqOperators, Test

@testset "Auxiliary functions for GhostDerivativeOperator" begin
    # check length and size functions
    L = CenteredDifference(2, 2, 1., 10)
    Q1 = RobinBC((-1.2, .3, .5), (.2, .5, -30.5), 1.)
    Q2 = PeriodicBC(Float64)

    A = L * Q1
    B = L * Q2

    @test ndims(A) == 2
    @test size(A) == (10, 10)
    @test size(A, 1) == 10
    @test length(A) == 100

    @test ndims(B) == 2
    @test size(B) == (10, 10)
    @test size(B, 1) == 10
    @test length(B) == 100
end


text:
using JuMP
import MathProgBase

m = Model()
@variable(m, x1)
@NLconstraint(m, cons1, x1^2 <= 1)
@NLconstraint(m, cons2, x1^2 >= 2)
@NLconstraint(m, cons2, x1^2 == 2)
z = [1.0]

d = JuMPNLPEvaluator(m)
MathProgBase.initialize(d, [:Grad])

g_val= zeros(3)
MathProgBase.eval_g(d, g_val, z)
@show g_val

I,J = MathProgBase.jac_structure(d)
jac_val = zeros(size(J))
MathProgBase.eval_jac_g(d, jac_val, z)
jac = full(sparse(I, J, jac_val))
@show jac

lb_consts, ub_consts = JuMP.getConstraintBounds(m)

println("----")