text:
using Test
using moomin
using JuMP
using DelimitedFiles

solvers = [:CPLEX; :Gurobi]
solverAvailable = falses(length(solvers))
for i in eachindex(solvers)
    try
        eval(:(using $(solvers[i])))
        solverAvailable[i] = true
        @info "$(solvers[i]) available."
    catch
    end
end

@testset "Solver availability" begin
    @test any(solverAvailable)
end

exampleData() = Data(["g1"; "g2"; "g3"; "g4"; "g5"],
                        [1; 0.95; 0; 0.5; 0.9],
                        [1; -1; -2; 3; -1])

exampleReactions() = Reactions(["r1"; "r2"; "r3"; "r4"; "r5"],                  #IDs
                               ["reac1"; "reac2"; "reac3"; "reac4"; "reac5"],   #names
                               zeros(5,5),                                      #S
                               fill(false, 5),                                  #reversible
                               ["g2"; "g6"; "g1"; "g3"; "g4"],                  #grRules
                               ["a"; "a"; "b"; "c"; "a"],                       #subSystems
                               [[[1]]; [[2;1]]; [[1;2;3]]; [[4;5]]; [[]]])      #assGenes

exampleGenes() = Genes(["g2"; "g6"; "g1"; "g3"; "g4"])

@testset "Types" begin
    data = exampleData()
    @test data isa Data
    reactions = exampleReactions()
    @test reactions isa Reactions
    genes = exampleGenes()
    @test genes isa Genes
    model = MoominModel(data, genes, reactions)
    @test model isa MoominModel
end

@testset "Weight function" begin
    @test weight(1, 0.9, 1, 7) == 4.605170
    @test weight(0.95, 0.9, 3, 7) == 1.3862940
    @test weight(0.5, 0.9, 1, 3) == -3.220
    @test weight(0, 0.5, 1, 7) == -1.3862940
end

@testset "Reading data from file" begin
    data = readData("data/toyData.txt", delimiter='\t', geneIDhead="GeneID",
        PPDEhead="PPDE", logFChead="logFC")
    @test data isa Data
    @test data.geneIDs == ["g1"; "g2"; "g3"; "g4"; "g5"]
    @test data.PPDE == [1; 0.95; 0; 0.5; 0.9]
    @test data.logFC == [1; -1; -2; 3; -1]
end

@testset "Mapping data" begin
    data = exampleData()
    genes = exampleGenes()
    mapDataToGenes!(genes, data)
    @test genes.PPDE == [0.95; 0; 1; 0; 0.5]
    @test genes.logFC == [-1; 0; 1; -2; 3]

    colour!(genes, 0.9, 1, 7)
    @test genes.colours == [-1; 0; 1; 0; 0]
    @test genes.weights == [1.386294;
                          -4.60517;
                          4.60517;
                          -4.60517;
                          -3.218876]

    reactions = exampleReactions()
    colour!(reactions, genes, 0.9, 1, 7)
    @test reactions.inputColours == [-1; -1; 0; 0; 0]
    @test reactions.weights == [1.386294;
                                1.386294;
                                -3.218876;
                                -3.218876;
                                -4.60517]
    @test reactions.leadingGenes == [1; 1; 0; 0; 0]
end

@testset "Reading models" begin
    rules = ["x(1)";
             "( x(2) | x(3) )";
             "(( x(2) & x(3) ) & ( x(4) & x(5) & x(6) ) & x(75) )";
             ""]
    @test getAssGenes(rules) == [[[1]]; [[2;3]]; [[2;3;4;5;6;75]]; [[]]]

    model = readModel("data/toyData.txt", "data/toyModel.mat",
    	delimiter='\t', geneIDhead="GeneID", PPDEhead="PPDE", logFChead="logFC",
    	modelStr="model",
    	pThresh=0.9, alpha=1, precision=7)
    @test model isa MoominModel
    @test model.reactions.inputColours == [-1; -1; 0; 0; 0; 0; 0]
    @test model.reactions.weights == [1.386294;
                                      1.386294;
                                      -3.218876;
                                      -3.218876;
                                      -4.60517;
                                      -4.60517;
                                      -4.60517]
end

@testset "MILP" begin
    for i in eachindex(solvers)
        if solverAvailable[i]
            @info "Testing $(solvers[i])..."
            if solvers[i] == :CPLEX
                optimizer = CPLEX.Optimizer
            elseif solvers[i] == :Gurobi
                optimizer = Gurobi.Optimizer
            end

            model = readModel("data/toyData.txt", "data/toyModel.mat")
            MILP = createMILP(model, optimizer, stoichiometry=true)
            @test lower_bound.(MILP[:v]) == fill(-100, 7)
            @test upper_bound.(MILP[:v]) == [100; 0; fill(100, 5)]

            xPlus = [1; 0; 0; 0; 1; 0]
            xMinus = [0; 1; 0; 0; 0; 1]
            inputColours = [0; 0; 0; 0; -1; 0]
            reversible = [false; false; false; false; true; true]
            @test interpretSolution(xPlus, xMinus, inputColours, reversible) ==
                [1; -1; 0; 0; -2; 6]

            outputColours = [1  0  0;
                             0  2  2;
                             -1 -2 0;
                             0  0  0;
                             1  1  1]
            fillOutputs!(model, outputColours)
            @test model.reactions.outputColours == outputColours
            @test model.reactions.outputFrequency == round.([1/3; 2/3; 2/3; 0; 1], digits=3)
            @test model.reactions.combinedOutput == [1; 2; 6; 0; 1]

            model = runMoomin("data/ecoliData.txt", "data/ecoli.mat", optimizer,
                        delimiter='\t', geneIDhead="GeneID", PPDEhead="PPDE", logFChead="logFC",
                        modelStr="model",
                        pThresh=0.9, alpha=1, precision=7,
                        stoichiometry=true, enumerate=10)
            @test size(model.reactions.outputColours, 2) == 1
            sol = readdlm("data/stoich_sol.txt")
            @test sol[:] == (model.reactions.outputColours .!= 0)

            model = runMoomin("data/ecoliData.txt", "data/ecoli.mat", optimizer,
                        stoichiometry=true, enumerate=10)
            @test size(model.reactions.outputColours, 2) == 3
            sol = readdlm("data/stoich_freq.txt")
            @test all(abs.(sol[:] .- model.reactions.outputFrequency) .< 0.01)

            model = runMoomin("data/ecoliData.txt", "data/ecoli.mat", optimizer,
                        stoichiometry=false, enumerate=10)
            @test size(model.reactions.outputColours, 2) == 1
            sol = readdlm("data/topo_sol.txt")
            @test sol[:] == (model.reactions.outputColours .!= 0)

            model = runMoomin("data/ecoliData.txt", "data/ecoli.mat", optimizer,
                        stoichiometry=false, alpha=1, enumerate=10)
            @test size(model.reactions.outputColours, 2) == 4
            sol = readdlm("data/topo_freq.txt")
            @test all(abs.(sol[:] .- model.reactions.outputFrequency) .< 0.01)
        end
    end
end

@testset "Printing output" begin
    data = exampleData()
    genes = exampleGenes()
    reactions = exampleReactions()
    mapDataToGenes!(genes, data)
    colour!(genes, 0.9, 1, 7)
    colour!(reactions, genes, 0.9, 1, 7)
    model = MoominModel(data, genes, reactions)
    outputColours = [ 1  0;
                     -1  2;
                     -2  0;
                      1  1;
                      6  1]
    fillOutputs!(model, outputColours)

    writeOutput(model, "data/outputExample")
    @test isfile("data/outputExample")
    output = readdlm("data/outputExample", header=true)
    @test output[2] == ["ID" "output"]
    @test output[1][:, 1] == ["r1"; "r2"; "r3"; "r4"; "r5"]
    @test output[1][:, 2] == ["red"; "blue"; "r.blue"; "red"; "yellow"]
    rm("data/outputExample")
    writeOutput(model, "data/outputExample", nSol=2, type="combined", asString=false)
    @test isfile("data/outputExample")
    output = readdlm("data/outputExample", header=true)
    @test output[2] == ["ID" "combined"]
    @test output[1][:, 2] == [1; 6; -2; 1; 6]
    rm("data/outputExample")

    writeOutput(model, "data/outputExample", format="full")
    @test isfile("data/outputExample")
    output = readdlm("data/outputExample", header=true)
    @test output[2] == ["ID" "input" "sol#1" "combined" "frequency" "subsystem" "leadingGene" "PPDE" "logFC"]
    @test output[1][:, 6] == ["a"; "a"; "b"; "c"; "a"]
    @test output[1][:, 7] == ["g2"; "g2"; "NA"; "NA"; "NA"]
    @test output[1][:, 8] == [0.95; 0.95; "NA"; "NA"; "NA"]
    @test output[1][:, 9] == [-1; -1; "NA"; "NA"; "NA"]
    rm("data/outputExample")

    writeOutput(model, "data/outputExample", format="json")
    @test isfile("data/outputExample")
    f = open("data/outputExample")
    output = read(f, String)
    close(f)
    @test output == "{\n\"r1\": red,\n\"r2\": blue,\n\"r3\": r.blue,\n\"r4\": red,\n\"r5\": yellow\n}"
    rm("data/outputExample")

    writeOutput(model, "data/outputExample", type="input")
    @test isfile("data/outputExample")
    rm("data/outputExample")
    writeOutput(model, "data/outputExample", type="frequency")
    @test isfile("data/outputExample")
    rm("data/outputExample")
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
# based on BSON.jl
function modpath(x::Module)
    y = parentmodule(x)
    result = x == y ? [nameof(x)] : [modpath(y)..., nameof(x)]

    if result[1] == :Main
        popfirst!(result)
    end

    return result
end

# based on BSON.jl
typepath(x::DataType) = [modpath(x.name.module)..., x.name.name]
typepathref(x::DataType) = string.(typepath(x))

# Returns a vector of tuples (nm, t),
# where nm is a Symbol for field name, and t a DataType.
function nametypetuples(t::DataType)
    @nospecialize t

    _fieldnames = fieldnames(t)
    _fieldtypes = fieldtypes(t)
    @assert length(_fieldnames) == length(_fieldtypes)
    return [ (_fieldnames[i], _fieldtypes[i]) for i in 1:length(_fieldnames) ]
end

function codegen_serialize(datatype::DataType) :: Expr
    @nospecialize datatype

    # "fieldname" => val.val.fieldname
    function field_value_pair_expr(nm::Symbol, @nospecialize(tt::Type{T})) :: Expr where {T}
        # "fieldname" => BSONSerializer.encode(val.val.fieldname, T)
        return :($("$nm") => BSONSerializer.encode(val.val.$nm, $T))
    end

    # "f1" => val.val.1, "f2" => val.val.f2, ...
    field_value_pairs = Expr(:tuple,
        [ field_value_pair_expr(nm, t) for (nm, t) in nametypetuples(datatype) ]...)

    type_encoded = typepathref(datatype)

    quote
        function BSONSerializer.serialize(val::BSONSerializer.Serializable{$datatype})
            return BSONSerializer.BSON("type" => $type_encoded, "args" => BSONSerializer.BSON($field_value_pairs...))
        end
   end
end

function codegen_deserialize(datatype::DataType) :: Expr
    @nospecialize datatype

    function arg_expr(nm::Symbol, @nospecialize(tt::Type{T})) :: Expr where {T}
        return :(BSONSerializer.decode( args[$("$nm")], $T))
    end

    # singleton or method
    if isdefined(datatype, :instance)
        return quote
            function BSONSerializer.deserialize(bson::Union{BSONSerializer.BSON, Dict}, ::Type{BSONSerializer.Serializable{$datatype}})
                ($datatype).instance
            end
        end
    else
        arg_list = Expr(:tuple,
            [ arg_expr(nm, t) for (nm, t) in nametypetuples(datatype) ]...)

        return quote
            function BSONSerializer.deserialize(bson::Union{BSONSerializer.BSON, Dict}, ::Type{BSONSerializer.Serializable{$datatype}})
                args = bson["args"]
                ($datatype)($arg_list...)
            end
        end
    end
end

is_type_reference(@nospecialize(m), s::Symbol) = isa(m.eval(s), DataType)

function is_type_reference(caller_module::Module, expr::Expr)
    @nospecialize caller_module expr

    is_module_name(expr::Symbol) = true
    is_module_name(expr::QuoteNode) = is_module_name(expr.value)
    is_module_name(other) = false
    function is_module_name(expr::Expr)
        if expr.head == :.
            @assert length(expr.args) == 2
            return is_module_name(expr.args[1]) && is_module_name(expr.args[2])
        end
    end

    is_type_name(expr::Symbol) = true
    is_type_name(expr::QuoteNode) = is_type_name(expr.value)
    is_type_name(other) = false

    # let's go slowly, because we're going to eval some expressions...
    if expr.head == :.
        @assert length(expr.args) == 2
        possibly_module_name = expr.args[1]
        possibly_type_name = expr.args[2]

        if is_module_name(possibly_module_name)
            type_owner_module = caller_module.eval(possibly_module_name)
            if isa(type_owner_module, Module)
                if is_type_name(possibly_type_name)
                    return isa(caller_module.eval(expr), DataType)
                end
            end
        end
    end
    return false
end

macro BSONSerializable(expr::Union{Expr, Symbol})
    @nospecialize expr

    #println("macro input: $expr, type $(typeof(expr))")

    if is_type_reference(__module__, expr)
        #println("$expr is a type reference.")
        datatype = __module__.eval(expr)
        expr_serialize_method = codegen_serialize(datatype)
        expr_deserialize_method = codegen_deserialize(datatype)

        return quote
            $expr_serialize_method
            BSONSerializer.encode(val::$datatype, ::Type{$datatype}) = BSONSerializer.serialize(BSONSerializer.Serializable(val))
            BSONSerializer.encode_type(::Type{$datatype}) = BSONSerializer.BSON

            $expr_deserialize_method
            BSONSerializer.decode(val::Union{BSONSerializer.BSON, Dict}, ::Type{$datatype}) = BSONSerializer.deserialize(val, BSONSerializer.Serializable{$datatype})
        end

    elseif isa(expr, Expr) && expr.head == :struct
        __module__.eval(expr) # a bit of a hack...
        struct_name = expr.args[2]
        return quote
            @BSONSerializable($struct_name)
        end

    else
        error("Couldn't apply @BSONSerialize to $expr.")
    end
end

ext: jl


text:
using SharedDataStructures
using Semaphores
using Base.Test
using BenchmarkTools

function deletewithlock!(D, key)
    withlock(D.lck) do
        SharedDataStructures.delete!(D, key)
    end
end

function setindexwithlock!(D, val, key)
    withlock(D.lck) do
        SharedDataStructures.setindex!(D, val, key)
    end
end

function getindexwithlock(D, key)
    withlock(D.lck) do
        SharedDataStructures.getindex(D, key)
    end
end

function test_shm_dict(keepkeys::Bool)
    path = pwd()
    master_shmdict = ShmDict(path, 10^3, 64; create=true, keepkeys=keepkeys)
    close(master_shmdict)
    delete!(master_shmdict)
    master_shmdict = ShmDict(path, 10^3, 64; create=true)

    try
        for N in 1:10^3
            for idx in 1:10^3
                str = string(idx)
                setindexwithlock!(master_shmdict, str, str)
                @test haskey(master_shmdict, str)
                @test getindexwithlock(master_shmdict, str) == convert(Vector{UInt8}, str)
            end
        end
        close(master_shmdict)
    finally
        delete!(master_shmdict)
    end
end

println("with no key stored:")
@btime test_shm_dict(false)
println("with key stored:")
@btime test_shm_dict(true)


text:
"""
    LookupArrays

Module for [`LookupArrays`](@ref) and [`Selector`]s used in DimensionalData.jl

`LookupArrays` defines traits and `AbstractArray` wrappers
that give specific behaviours for a lookup index when indexed with [`Selector`](@ref).

For example, these allow tracking over array order so fast indexing works evne when
the array is reversed.

To load LookupArrays types and methods into scope:

```julia
using DimensionalData
using .LookupArrays
```
"""
module LookupArrays

using Dates, IntervalSets
import Adapt, ConstructionBase

using Base: tail, OneTo, @propagate_inbounds

export order, sampling, span, bounds, locus, hasselection, transformdim,
    metadata, units, sort, selectindices, val, index, reducelookup, shiftlocus,
    maybeshiftlocus

export Selector
export At, Between, Contains, Near, Where, All
export ..

export LookupArrayTrait
export Order, Ordered, ForwardOrdered, ReverseOrdered, Unordered, AutoOrder
export Sampling, Points, Intervals, AutoSampling, NoSampling
export Span, Regular, Irregular, Explicit, AutoSpan, NoSpan
export Locus, Center, Start, End, AutoLocus
export Metadata, NoMetadata
export AutoStep, AutoBounds, AutoIndex

export LookupArray
export AutoLookup, NoLookup
export Aligned, AbstractSampled, Sampled, AbstractCategorical, Categorical
export Unaligned, Transformed

const StandardIndices = Union{AbstractArray{<:Integer},Colon,Integer,CartesianIndex}

rebuild(x; kw...) = ConstructionBase.setproperties(x, (; kw...))

include("metadata.jl")
include("lookup_traits.jl")
include("lookup_arrays.jl")
include("selector.jl")
include("indexing.jl")
include("methods.jl")
include("utils.jl")
include("set.jl")
include("show.jl")

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of basic streams
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael de P. Soares
#* $Id$
#*---------------------------------------------------------------------
type simple_source
	simple_source()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Protected=>true
			)),(NComp)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Valid Phases for Flash Calculation",
				:Valid=>["Vapour-Only", "Liquid-Only","Vapour-Liquid"],
				:Default=>"Vapour-Liquid"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.5256,
				:Symbol=>"_{out}",
				:Protected=>true
			)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Stream Molar Composition"
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Sum of Stream Composition",
				:Protected=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Stream Molar Flow Rate"
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Stream Temperature"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Temperature in  C",
				:Lower=>-250,
				:Upper=>5000
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Stream Pressure"
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Fraction",
				:Hidden=>true
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar Fraction",
				:Hidden=>true
			)),(NComp)),
			[
				:(SumOfComposition = sum(MolarComposition)),
				:(Outlet.z = MolarComposition/sum(MolarComposition)),
				:(Outlet.v = 0),
				:(x = Outlet.z),
				:(y = Outlet.z),
				:(Outlet.h = PP.LiquidEnthalpy(Outlet.T, Outlet.P, x)),
				:(Outlet.v = 1),
				:(x = Outlet.z),
				:(y = Outlet.z),
				:(Outlet.h = PP.VapourEnthalpy(Outlet.T, Outlet.P, y)),
				:([Outlet.v, x, y] = PP.Flash(Outlet.T, Outlet.P, Outlet.z)),
				:(Outlet.h = (1-Outlet.v)*PP.LiquidEnthalpy(Outlet.T, Outlet.P, x) + Outlet.v*PP.VapourEnthalpy(Outlet.T, Outlet.P, y)),
				:(T_Cdeg = Outlet.T/"K" - 273.15),
				:(Outlet.F = F),
				:(Outlet.P = P),
				:(Outlet.T = T),
			],
			[
				"Sum of Composition","Stream Molar Composition","Vapour Fraction","Liquid Composition","Vapour Composition","Overall Enthalpy","Vapor Fraction","Liquid Composition","Vapour Composition","Overall Enthalpy","Flash Calculation","Overall Enthalpy","Temperature in  C","Equate Flow","Equate Pressures","Equate Temperatures",
			],
			[:PP,:NComp,:M,:ValidPhases,],
			[:Outlet,:MolarComposition,:SumOfComposition,:F,:T,:T_Cdeg,:P,:x,:y,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	M::Array{molweight}
	ValidPhases::DanaSwitcher
	Outlet::stream
	MolarComposition::Array{positive}
	SumOfComposition::positive
	F::flow_mol
	T::temperature
	T_Cdeg::DanaReal
	P::pressure
	x::Array{fraction}
	y::Array{fraction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export simple_source
function set(in::simple_source)
	M = PP.MolecularWeight()

end
function setEquationFlow(in::simple_source)
	addEquation(1)
	addEquation(2)
	let switch=ValidPhases
		if switch=="Liquid-Only"
			addEquation(3)
			addEquation(4)
			addEquation(5)
			addEquation(6)
		elseif switch=="Vapour-Only"
			addEquation(7)
			addEquation(8)
			addEquation(9)
			addEquation(10)
		elseif switch=="Vapour-Liquid"
			addEquation(11)
			addEquation(12)
		end
	end
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
end
function atributes(in::simple_source,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Source"
	fields[:Brief]="Simple Material stream source"
	fields[:Info]="
	This model should be used for boundary streams.
	Usually these streams are known and come from another process
	units.

	The user should specify:
	 * Total molar flow
	 * Temperature
	 * Pressure
	 * Molar composition
"
	drive!(fields,_)
	return fields
end
simple_source(_::Dict{Symbol,Any})=begin
	newModel=simple_source()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# CUDA forward compatibility package as taken from the datacenter driver
#
# - https://docs.nvidia.com/deploy/cuda-compatibility/index.html#forward-compatibility-title
# - https://docs.nvidia.com/datacenter/tesla/index.html
# - https://www.nvidia.com/Download/index.aspx

using BinaryBuilder, Pkg

include("../../../fancy_toys.jl")

name = "CUDA_compat"
version = v"0.0.1"

dependencies = []

script = raw"""
sh installer.run -x
cd NVIDIA-*
find .

install_license LICENSE

mkdir -p ${libdir}

ln -s libcuda.so.* libcuda.so
mv libcuda.so* ${libdir}

ln -s libnvidia-ptxjitcompiler.so.* libnvidia-ptxjitcompiler.so
mv libnvidia-ptxjitcompiler.so* ${libdir}
"""

products = [
    LibraryProduct("libcuda", :libcuda),
    LibraryProduct("libnvidia-ptxjitcompiler", :libnvidia_ptxjitcompiler),
]

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

cuda_versions = [v"11.4"]
for cuda_version in cuda_versions
    cuda_tag = "$(cuda_version.major).$(cuda_version.minor)"
    include("build_$(cuda_tag).jl")
end

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
using PyPlot

if !isinteractive()
    DIR = @__DIR__
else
    # assume root project
    DIR = "./tutorial"
end

@everywhere begin
@everywhere @eval DIR=$DIR
using WaveletInformation
using Distributions, StatsBase, MultivariateStats, Wavelets

resolution = 0.001 # 1 ms
τ = range(-8.192, # TODO this shouldn't be global
          stop=8.192,
          step=resolution)
n_trials = 200

include("$DIR/decoders.jl")
include("$DIR/neurons.jl")
include("$DIR/tutorial_lib.jl")

neuron_responses = [
                    NeuronResponse("Gaussian around stimulus",
                                   get_spikes_from_distribution(Normal(0,2))),
                    NeuronResponse("Gaussian offset from stimulus",
                                   get_spikes_from_distribution(Normal(2,2))),
                    NeuronResponse("Ramp after stimulus",
                                   N -> wsample(τ,ramp(τ),N)),
                    NeuronResponse("Oscillate after stimulus",
                                   N -> wsample(τ,oscillate(τ),N)),
                    NeuronResponse("Short response and refractory",
                                   N -> wsample(τ,refractory(τ),N))
                   ]

decoders = [
            ResponseDecoder("Spike Count",
                            spike_count_score),
            ResponseDecoder("PCA",
                            pca_score),
            ResponseDecoder("Wavelet Information",
                            (x,y)->wavelet_decoder_score(x,y,WIopts(50, :max,
                                                                    95, 5, 2)
                                                        ))
           ]

num_spikes_array = [10 .^ (range(1,step=1,stop=5)) 5*
                    (10 .^ (range(1,step=1,stop=5)))]'[:]

num_trials_array = [10 .^ (range(1,step=1,stop=3)) 5*
                    (10 .^ (range(1,step=1,stop=3)))]'[1:end-1]
end

modulation_scores = modulation_score(neuron_responses[1:1],
                                     decoders,
                                     1000, #spikes
                                     100, #trials
                                     10, #surrogates
                                    )

@time begin
modulation_scores_by_num_spikes = pmap(num_spikes->
                                       modulation_score(neuron_responses,
                                     decoders,
                                     num_spikes, #spikes
                                     100, #trials
                                     10, #surrogates
                                    ), num_spikes_array)
end
#
@time begin
modulation_scores_by_num_trials = pmap(num_trials->
                                       modulation_score(neuron_responses,
                                     decoders,
                                     1000, #spikes
                                     num_trials, #trials
                                     10, #surrogates
                                    ), num_trials_array)
end
#"/tmp/julia6n8pMj"

function plot_modulation!(ax, x, modulation, decoders)
    for (di,decoder) in enumerate(decoders)
        m = [m[decoder.name] for m in modulation]
        mean_m = [mean(mx) for mx in m]
        ax[:plot](x, mean_m,
                             label=decoder.name,
                             color=decoder_colors[di])
        for (oi,obs) in enumerate(m)
            ax[:scatter](repeat(x[oi:oi], length(obs)),
                         obs,
                         s=0.8,
                         color=decoder_colors[di],
                         alpha=0.5)
        end
        ax[:set_xscale]("log")
    end
    ax[:set_ylim](0,1.1)
end

include("$DIR/plotting.jl")
plt[:style][:use]("/home/data/.config/matplotlib/stylelib/alex.mplstyle")
ioff()
decoder_colors = [[156,0,255]./256,[9,178,85]./256,[255,173,25]./256]
f,axs = subplots(4,length(neuron_responses),
                 figsize = (7,5),
                 gridspec_kw = Dict(:height_ratios=>[2,2,3,3]),
                 sharey="row"
                )
for (ni,neuron_response) in enumerate(neuron_responses)
    raster_ax = axs[1,ni]
    spikes = [neuron_response.get_spikes(100) for _ in 1:10]
    raster_plot!(raster_ax, spikes)
    #
    histogram_ax = axs[2,ni]
    histogram_ax[:hist](reduce(vcat,spikes),bins=50)
    #
    modulation_spikes_ax = axs[3,ni]
    modulation = [m[neuron_response.name]
                  for m in modulation_scores_by_num_spikes]
    plot_modulation!(modulation_spikes_ax, num_spikes_array,
                     modulation, decoders)
    modulation_trials_ax = axs[4,ni]
    modulation = [m[neuron_response.name]
                  for m in modulation_scores_by_num_trials]
    plot_modulation!(modulation_trials_ax, num_trials_array,
                     modulation, decoders)
    #
    if ni == 1
        raster_ax[:set_xlabel]("Stimulus Onset (s)")
        raster_ax[:set_yticks]([])
        histogram_ax[:set_xlabel]("Stimulus Onset (s)")
        histogram_ax[:set_ylabel]("Spike Count")
        modulation_spikes_ax[:set_xlabel]("Number of Spikes")
        modulation_spikes_ax[:set_ylabel]("Decoder\nAccuracy")
        modulation_trials_ax[:set_xlabel]("Number of Trials")
        modulation_trials_ax[:set_ylabel]("Decoder\nAccuracy")
    else
        nothing
    end
end
f[:tight_layout]()
f[:subplots_adjust](hspace=0.5)
if true#false
f[:savefig]("/mnfs/vtad1/data/amorley_figs/Thesis/wavelet_decoder_comparison.svg",
           )
f[:savefig]("/mnfs/vtad1/data/amorley_figs/Thesis/wavelet_decoder_comparison.png",
           )
end
#=
    #
    x1 = spike_count_score(activity_matrix_split, before_after_identifiers)
    #
    opts = WaveletInformation.WIopts(50,:max, 95, 5, 2)
    x2 = wavelet_decoder_score(activity_matrix_split, before_after_identifiers,
                               WaveletInformation.WIopts(50,:max, 95, 5, 2))
    println("""
            Spike Count Score: $x1
            Wavelet Decoder Score: $x2
            """
           )
=#


text:
"""
Implement the `Data.Source` interface in the `DataStreams.jl` package.

"""
mutable struct Source{I} <: Data.Source
    schema::Data.Schema
    options::Options
    io::I
    fullpath::String
    datapos::Int # the position in the IOBuffer where the rows of data begins
    currentline::Vector{Union{Missing,String}}
    lastcol::Int
    malformed::Bool # file is malformed so slow parsing
    eolpad::Int #Number of characters to pad at end of line (typically=0, 1 for CRLF files)
    line_len::Int # first line length successfully parsed by readsplitline!
end

function Base.show(io::IO, f::Source)
    println(io, "FWF.Source: ", f.fullpath)
    println(io, "Currentline   : ")
    show(io, f.currentline)
    println(io)
    show(io, f.options)
    show(io, f.schema)
end

function union_missing(m::Bool, t::Type)
    m ? (Union{t, Missing}) : (t)
end

# Negative values will break these functions

# function row_calc(lines::Int, rows::Int, skip::Int, header::Bool)
#     return row_calc(lines, rows, skip) - (header ? 1 : 0)
# end

# function row_calc(lines::Int, rows::Int, skip::Int, header::T) where {T}
#     return row_calc(lines, rows, skip)
# end

# function row_calc(lines::Int, rows::Int, skip::Int)
#     # rows to process, subtract skip and header if they exist
#     rows = rows <= 0 ?  lines : ( (lines < rows) ? (lines) : (rows))
#     return skip > 1 ? rows - skip : rows
# end

function calculate_ranges(columnwidths::Union{Vector{UnitRange{Int}}, Vector{Int}})
    rangewidths = Vector{UnitRange{Int}}(undef, length(columnwidths))
    if isa(columnwidths, Vector{Int})
        l = 0
        for i in eachindex(columnwidths)
            columnwidths[i] < 1 && (throw(ArgumentError("A column width less than 1")))
            rangewidths[i] = l+1:l+columnwidths[i]
            l=last(rangewidths[i])
        end
    else
        #Validate the passed ranges
        first(first(columnwidths)) <= 0 && (throw(ArgumentError("Columns must start > 0")))
        for i in 1:length(columnwidths)
            first(columnwidths[i]) ≤ last(columnwidths[i]) || throw(ArgumentError(string("Negative range ", columnwidths[i])))
            rangewidths[i] = columnwidths[i]
            i==1 && (continue)
            (last(columnwidths[i-1]) < first(columnwidths[i])) || (throw(ArgumentError("Non-increasing ranges "*string(columnwidths[i-1])*","*string(columnwidths[i]))))
        end
    end
    return rangewidths
end

# Create a source data strcuture.  To do this we need to do the following
# * Ensure file exists and open for reading
# * Determine column names
# * Determine column types
# * Determine number of rows we will read
# * Convert missings to a dictionary, if applicable
# * Convert date formats to a dictionary, if applicable

function Source(
    fullpath::Union{AbstractString, IO},
    columnwidths::Union{Vector{UnitRange{Int}}, Vector{Int}}
    ;
    usemissings::Bool=true,
    trimstrings::Bool=true,
    errorlevel::Symbol=:parse,
    unitbytes::Bool=true,
    use_mmap::Bool=true,
    skip::Int=0,
    rows::Int=0,
    types::Vector=Vector(),
    header::Union{Bool, Vector{String}}=Vector{String}(),
    missings::Vector{String}=Vector{String}()
    )
    # Attempting to re-create all objects here to minimize outside tampering
    datedict = Dict{Int, DateFormat}()
    typelist = Vector{DataType}()
    rangewidths = Vector{UnitRange{Int}}()
    malformed = false

    isa(fullpath, AbstractString) && (isfile(fullpath) || throw(ArgumentError("\"$fullpath\" is not a valid file")))

    # open the file and prepare for procesing
    if isa(fullpath, IOBuffer)
        source = fullpath
        fs = nb_available(fullpath)
        fullpath = "<IOBuffer>"
    elseif isa(fullpath, IO)
        source = IOBuffer(Base.read(fullpath))
        fs = nb_available(source)
        fullpath = isdefined(fullpath, :name) ? fullpath.name : "__IO__"
    else
        source = open(fullpath, "r") do f
            IOBuffer(use_mmap ? Mmap.mmap(f) : Base.read(f))
        end
        fs = filesize(fullpath)
    end

    # Number of columns = # of widths
    isempty(columnwidths) && throw(ArgumentError("No column widths provided"))
    columns = length(columnwidths)

    rangewidths = calculate_ranges(columnwidths)

    # Starting position
    startpos = position(source)
    # rows to process, subtract skip and header if they exist
    lines, eolpad =  row_countlines(source)
    skip < 0 && (skip = 0)
    rows -= skip
    # rows = row_calc(lines, rows,skip, header)
    lines ≤ (isa(header, Bool) && header) + skip && (throw(ArgumentError("More skips than rows available")))
    # Go back to start
    seek(source, startpos)

    # Don't think this is necessary, but just in case...BOM character detection
    if fs > 0 && Base.peek(source) == 0xef
         read(source, UInt8)
         read(source, UInt8) == 0xbb || seek(source, startpos)
         read(source, UInt8) == 0xbf || seek(source, startpos)
    end

    # reposition iobuffer
    tmp = skip
    while (!eof(source)) && (tmp > 0)
        readline(source)
        tmp -= 1
    end
    datapos = position(source)
    line_len = -1

    # Figure out headers
    if isa(header, Bool) && header
        # first row is heders
        header_tmp = Vector{Union{Missing,String}}()
        line_len = FWF.readsplitline!(header_tmp, source, rangewidths, true, unitbytes)
        headerlist = Vector{String}(length(header_tmp))
        datapos = position(source)
        for i in eachindex(headerlist)
            if ismissing(header_tmp[i]) || isempty(header_tmp[i])
                headerlist[i] = "Column$i"
            else
                headerlist[i] = header_tmp[i]
            end
        end
    elseif (isa(header, Bool) && !header) || isempty(header)
        # number columns
        headerlist = ["Column$i" for i = 1:columns]
    elseif !isempty(header)
        length(header) != columns && (throw(ArgumentError("Header count doesn't match column count")))
        headerlist = copy(header)
    else
        throw(ArgumentError("Can not determine headers"))
    end

    # Type is set to String if types are not passed in
    # Otherwise iterate through copying types & creating date dictionary
    if isempty(types)
        typelist = [union_missing(usemissings, String) for i = 1:columns]
    else
        length(types) != columns && throw(ArgumentError("Wrong number of types: "*string(length(types))))
        @static if VERSION < v"0.7.0-DEV"
            typelist = Vector{Type}(columns)
        else
            typelist = Vector{Type}(uninitialized, columns)
        end
        for i in 1:length(types)
            if (isa(types[i], DateFormat))
                typelist[i] = union_missing(usemissings, Date)
                datedict[i] = types[i]
            elseif (isa(types[i], Type))
                !(types[i] in (Int, Float64, String, Missing)) && (throw(ArgumentError("Invalid Type: "*string(types[i]))))
                isa(types[i], Missing) ? Missing : typelist[i] = union_missing(usemissings, types[i])
            else
               throw(ArgumentError("Found type that is not a DateFormat or DataType"))
            end
        end
    end

    # Convert missings to Set for faster lookup later.
    missingset=Set{String}(missings)

    sch = Data.Schema(typelist, headerlist, ifelse(rows ≤ 0, missing, rows))
    opt = Options(usemissings=usemissings, trimstrings=trimstrings,
                    errorlevel=errorlevel, unitbytes=unitbytes, skip=skip, missingvals=missingset,
                    dateformats = datedict,
                    columnrange=rangewidths)
    return Source(sch, opt, source, string(fullpath), datapos, Vector{Union{Missing,String}}(), 0, malformed, eolpad, line_len)
end

# needed? construct a new Source from a Sink
#Source(s::CSV.Sink) = CSV.Source(fullpath=s.fullpath, options=s.options)
Data.reset!(s::FWF.Source) = (seek(s.io, s.datapos); return nothing)
Data.schema(source::FWF.Source) = source.schema
Data.accesspattern(::Type{<:FWF.Source}) = Data.Sequential
@inline Data.isdone(io::FWF.Source, row, col, rows, cols) = eof(io.io) || (!ismissing(rows) && row > rows)
#@inline Data.isdone(io::Source, row, col) = Data.isdone(io, row, col, size(io.schema)...)
Data.streamtype(::Type{<:FWF.Source}, ::Type{Data.Column}) = false
Data.streamtype(::Type{<:FWF.Source}, ::Type{Data.Field}) = true
@inline Data.streamfrom(source::FWF.Source, ::Type{Data.Field}, ::Type{T}, row, col::Int) where {T} = FWF.parsefield(source, T, row, col)
Data.streamfrom(source::FWF.Source, ::Type{Data.Column}, ::Type{T}, col::Int) where {T} = FWF.parsecol(source, T, col)
Data.reference(source::FWF.Source) = source.io.data


text:
module Blog

using JSServe
using Markdown
using Hyperscript
using JSServe: Asset, ES6Module, AssetFolder, Routes
using GitHub, Downloads

include("gh-utils.jl")

site_path(files...) = normpath(joinpath(@__DIR__, "..", "docs", files...))
markdown(files...) = joinpath(@__DIR__, "pages", "blogposts", files...)
asset(files...) = Asset(normpath(joinpath(@__DIR__, "assets", files...)))
const Highlight = ES6Module(joinpath(@__DIR__,  "assets", "libs", "highlight", "highlight.pack.js"))

function list(elements)
    style = "display: flex; flex-direction: column; justify-content: start; align-items: start"
    return DOM.div(JSServe.jsrender.(elements)..., style=style)
end

function make_app(dom)
    return JSServe.App() do
        assets = asset.([
            "css/franklin.css",
            "css/makie.css",
            "css/minimal-mistakes.css",
            "libs/highlight/github.min.css",
            "css/style.css"])

        highlight = DOM.div(
            DOM.script(src=asset("libs", "highlight", "highlight.pack.js")),
        DOM.script("hljs.highlightAll()")
        )
        return DOM.html(
            DOM.head(
                DOM.meta(charset="UTF-8"),
                DOM.meta(name="viewport", content="width=device-width, initial-scale=1"),
                assets...,
                DOM.link(rel="icon", type="image/x-icon", href=asset("images", "favicon.ico")),
            ),
            DOM.body(dom, highlight)
        )
    end
end

function page(file)
    source = read(file, String)
    md = JSServe.string_to_markdown(source, Blog; eval_julia_code=Blog)
    banner = DOM.a(DOM.img(src = asset("images", "bannermesh_gradient.png")), href="/")
    body = DOM.div(DOM.div(md, class="inner-page"), class="outer-page")
    return make_app(DOM.div(banner, body))
end

function make()
    src = readdir(markdown(), join=true)
    routes = Routes()
    for file in src
        name, ext = splitext(basename(file))
        if ext == ".md"
            routes["blogposts/" * name] = page(file)
        end
    end
    routes["/"] = page(markdown("..", "index.md"))
    JSServe.export_static(site_path(), routes)
end

end

max_stars_repo_name: MakieOrg/Blog

max_issues_repo_name: MakieOrg/Blog

max_forks_repo_name: MakieOrg/Blog


text:
"""
```julia
randToeplitz(d, n;  norm, hermitian)
```
- `d` : entry distribution
- `n` : dimension
- `norm` : default `false`;  if `norm` set to `true`, then the matrix will be normlaized with ``n^{-1/2}``.
- `hermitian`: default `true`; if `true` the matrix will be Hermitian

# Examples

Generate a ``4 \\times 4`` random Hermitian Toeplitz matrix with entries Standard Normal.
```julia
randToeplitz(4)

4×4 Matrix{Float64}:
  1.10207   -0.47292   -0.745498   1.06809
 -0.47292    1.10207   -0.47292   -0.745498
 -0.745498  -0.47292    1.10207   -0.47292
  1.06809   -0.745498  -0.47292    1.10207
```
Generate a ``4 \\times 4`` normalized random Toeplitz matrix with entries `Exponential(1)`.
```julia
using Distributions
randToeplitz(Exponential(1),4, norm = true, hermitian = false)

4×4 Matrix{Float64}:
 0.667888  0.260045  1.48812   0.477305
 1.50374   0.667888  0.260045  1.48812
 1.1475    1.50374   0.667888  0.260045
 0.363966  1.1475    1.50374   0.667888
```
"""
function  randToeplitz(d::D, n::Int;  norm = false::Bool, hermitian=true::Bool)  where D<:S
    if hermitian
      t = rand(d,n)
      t = [t[n:-1:2]'..., t...]
    else
      t = rand(d,2n-1)
    end
    M = zeros(eltype(d),n,n)
    for i in 1:n,j in 1:n
        M[i,j]=t[i-j+n]
    end

    if norm  return M/sqrt(n) end
    M

end

function randToeplitz(n::Int; norm=false::Bool, hermitian=true::Bool)
  return randToeplitz(Normal(),n,norm=norm,hermitian=hermitian)
end

"""
```julia
randHankel(d, n;  norm )

randHankel(n;  norm)
```
- `d` : entry distribution
- `n` : dimension
- `norm` : default `false`,  if `norm` set to `true`, then the matrix will be normlaized with ``n^{-1/2}``.

# Examples

Generate a ``5\\times 5`` random Hankel matrix with entries uniformly distributed on ``\\{1, i, \\pi \\}``
```julia
randHankel((1,im,pi),5)

5×5 Matrix{Number}:
  1   1  im  1   1
  1  im   1  1   π
 im   1   1  π   π
  1   1   π  π   π
  1   π   π  π  im
```

"""
function randHankel(d::D, n::Int;  norm = false::Bool)  where D<:S
  h = rand(d,2n-1)
  M = zeros(eltype(d),n,n)
  for i in 1:n
    M[i,:] = h[i:i+n-1]
  end

  if norm  return M/sqrt(n) end
  M
end

function randHankel(n::Int;  norm = false::Bool)
  return randHankel(Normal(),n,norm=norm)
end


text:
# Where node.

mutable struct WhereNode <: TabularNode
    over::Union{SQLNode, Nothing}
    condition::SQLNode

    WhereNode(; over = nothing, condition) =
        new(over, condition)
end

WhereNode(condition; over = nothing) =
    WhereNode(over = over, condition = condition)

"""
    Where(; over = nothing, condition)
    Where(condition; over = nothing)

The `Where` node filters the input rows by the given `condition`.

```sql
SELECT ...
FROM \$over
WHERE \$condition
```

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person) |>
           Where(Fun(">", Get.year_of_birth, 2000));

julia> print(render(q))
SELECT "person_1"."person_id", "person_1"."year_of_birth"
FROM "person" AS "person_1"
WHERE ("person_1"."year_of_birth" > 2000)
```
"""
Where(args...; kws...) =
    WhereNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Where), pats::Vector{Any}) =
    dissect(scr, WhereNode, pats)

function PrettyPrinting.quoteof(n::WhereNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Where), quoteof(n.condition, ctx))
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::WhereNode) =
    label(n.over)

rebase(n::WhereNode, n′) =
    WhereNode(over = rebase(n.over, n′), condition = n.condition)


text:
doc"""
    NegativeBinomial(r,p)

A *Negative binomial distribution* describes the number of failures before the `r`th success in a sequence of independent Bernoulli trials. It is parameterized by `r`, the number of successes, and `p`, the probability of success in an individual trial.

$P(X = k) = {k + r - 1 \choose k} p^r (1 - p)^k, \quad \text{for } k = 0,1,2,\ldots.$

The distribution remains well-defined for any positive `r`, in which case

$P(X = k) = \frac{\Gamma(k+r)}{k! \Gamma(r)} p^r (1 - p)^k, \quad \text{for } k = 0,1,2,\ldots.$

```julia
NegativeBinomial()        # Negative binomial distribution with r = 1 and p = 0.5
NegativeBinomial(r, p)    # Negative binomial distribution with r successes and success rate p

params(d)       # Get the parameters, i.e. (r, p)
succprob(d)     # Get the success rate, i.e. p
failprob(d)     # Get the failure rate, i.e. 1 - p
```

External links:

* [Negative binomial distribution on Wikipedia](http://en.wikipedia.org/wiki/Negative_binomial_distribution)

"""

immutable NegativeBinomial{T<:Real} <: DiscreteUnivariateDistribution
    r::T
    p::T

    function (::Type{NegativeBinomial{T}}){T}(r::T, p::T)
        @check_args(NegativeBinomial, r > zero(r))
        @check_args(NegativeBinomial, zero(p) < p <= one(p))
        new{T}(r, p)
    end

end

NegativeBinomial{T<:Real}(r::T, p::T) = NegativeBinomial{T}(r, p)
NegativeBinomial(r::Real, p::Real) = NegativeBinomial(promote(r, p)...)
NegativeBinomial(r::Integer, p::Integer) = NegativeBinomial(Float64(r), Float64(p))
NegativeBinomial(r::Real) = NegativeBinomial(r, 0.5)
NegativeBinomial() = NegativeBinomial(1.0, 0.5)

@distr_support NegativeBinomial 0 Inf

#### Conversions

function convert{T<:Real}(::Type{NegativeBinomial{T}}, r::Real, p::Real)
    NegativeBinomial(T(r), T(p))
end
function convert{T <: Real, S <: Real}(::Type{NegativeBinomial{T}}, d::NegativeBinomial{S})
    NegativeBinomial(T(d.r), T(d.p))
end

#### Parameters

params(d::NegativeBinomial) = (d.r, d.p)
@inline partype{T<:Real}(d::NegativeBinomial{T}) = T

succprob(d::NegativeBinomial) = d.p
failprob(d::NegativeBinomial) = 1 - d.p

#### Statistics

mean(d::NegativeBinomial) = (p = succprob(d); (1 - p) * d.r / p)

var(d::NegativeBinomial) = (p = succprob(d); (1 - p) * d.r / (p * p))

std(d::NegativeBinomial) = (p = succprob(d); sqrt((1 - p) * d.r) / p)

skewness(d::NegativeBinomial) = (p = succprob(d); (2 - p) / sqrt((1 - p) * d.r))

kurtosis(d::NegativeBinomial) = (p = succprob(d); 6 / d.r + (p * p) / ((1 - p) * d.r))

mode(d::NegativeBinomial) = (p = succprob(d); floor(Int,(1 - p) * (d.r - 1) / p))

#### Evaluation & Sampling

@_delegate_statsfuns NegativeBinomial nbinom r p

rand(d::NegativeBinomial) = convert(Int, StatsFuns.RFunctions.nbinomrand(d.r, d.p))

immutable RecursiveNegBinomProbEvaluator <: RecursiveProbabilityEvaluator
    r::Float64
    p0::Float64
end

RecursiveNegBinomProbEvaluator(d::NegativeBinomial) = RecursiveNegBinomProbEvaluator(d.r, failprob(d))
nextpdf(s::RecursiveNegBinomProbEvaluator, p::Float64, x::Integer) = ((x + s.r - 1) / x) * s.p0 * p
_pdf!(r::AbstractArray, d::NegativeBinomial, rgn::UnitRange) = _pdf!(r, d, rgn, RecursiveNegBinomProbEvaluator(d))

function mgf(d::NegativeBinomial, t::Real)
    r, p = params(d)
    return ((1 - p) * exp(t))^r / (1 - p * exp(t))^r
end

function cf(d::NegativeBinomial, t::Real)
    r, p = params(d)
    return (((1 - p) * cis(t)) / (1 - p * cis(t)))^r
end


text:
#################################
# Weather Alarms API            #
#################################
# Get Weather Alarms for a specific location. AccuWeather Weather Alarms are determined using the daily forecasts for a location. An alarm exists for a location if the forecast weather meets or exceeds specific thresholds (https://developer.accuweather.com/weather-alarm-thresholds).

# Contains a centralized function to make calls to the API.
include("request.jl")

function get_weatheralarms(apikey::String, locationKey::String, time::Int64, language::String)
    # Returns 1/5/10/15 day(s) of weather alarms for a specific location.

    url = "http://dataservice.accuweather.com/alarms/v1/"

    # Add required parameters to the URL
    url = url * string(time) * "day/"
    url = url * locationKey
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end

    return make_API_call(url)
end


text:
@enum SETTING_IDENTIFIER SETTINGS_HEADER_TABLE_SIZE=0x1 SETTINGS_ENABLE_PUSH=0x2 SETTINGS_MAX_CONCURRENT_STREAMS=0x3 SETTINGS_INITIAL_WINDOW_SIZE=0x4 SETTINGS_MAX_FRAME_SIZE=0x5 SETTINGS_MAX_HEADER_LIST_SIZE=0x6

struct SettingsFrame
    is_ack::Bool
    parameters::Union{Nothing, Array{Tuple{SETTING_IDENTIFIER, UInt32}, 1}}
end

SettingsFrame() = SettingsFrame(false,  Array{Tuple{SETTING_IDENTIFIER, UInt32}, 1}())

==(a::SettingsFrame, b::SettingsFrame) =
    a.is_ack == b.is_ack &&
    (isnothing(a.parameters) || a.parameters == b.parameters)

mutable struct UnknownIdentifierError <: Exception end

function decode_settings(header, payload)
    @assert header.stream_identifier == 0x0

    is_ack = header.flags & 0x1 == 0x1

    if is_ack
        @assert length(payload) == 0
        return SettingsFrame(is_ack, nothing)
    else
        parameters = Array{Tuple{SETTING_IDENTIFIER, UInt32}, 1}()
        for i = 1:div(length(payload), 6)
            identifier = UInt16(payload[(i-1)*6+1]) << 8 + UInt16(payload[(i-1)*6+2])
            value = UInt32(payload[(i-1)*6+3]) << 24 + UInt32(payload[(i-1)*6+4]) << 16 +
                UInt32(payload[(i-1)*6+5]) << 8 + UInt32(payload[(i-1)*6+6])
            (identifier > 0x6 || identifier < 0x1) && continue
            push!(parameters, (SETTING_IDENTIFIER(identifier), value))
        end
        return SettingsFrame(is_ack, parameters)
    end
end

function encode_settings(frame)
    if frame.is_ack
        return wrap_payload([], SETTINGS, 0x1, 0x0)
    else
        payload = Array{UInt8, 1}()
        for val in frame.parameters
            append!(payload, [ UInt8(UInt16(val[1]) >> 8);
                               UInt8(UInt16(val[1]) & 0x00ff);
                               UInt8(val[2] >> 24);
                               UInt8(val[2] >> 16 & 0x000000ff);
                               UInt8(val[2] >> 8 & 0x000000ff);
                               UInt8(val[2] & 0x000000ff) ])
        end
        return wrap_payload(payload, SETTINGS, 0x0, 0x0)
    end
end


text:
"""
    step!(integ::DEIntegrator [, dt [, stop_at_tdt]])
Perform one (successful) step on the integrator.

Alternative, if a `dt` is given, then `step!` the integrator until
there is a temporal difference `≥ dt` in `integ.t`.  When `true` is
passed to the optional third argument, the integrator advances exactly
`dt`.
"""
function step!(d::DEIntegrator) error("Integrator stepping is not implemented") end
Base.resize!(i::DEIntegrator,ii::Int) = error("resize!: method has not been implemented for the integrator")
Base.deleteat!(i::DEIntegrator,ii) = error("deleteat!: method has not been implemented for the integrator")
addat!(i::DEIntegrator,ii,val=zeros(length(idxs))) = error("addat!: method has not been implemented for the integrator")
get_tmp_cache(i::DEIntegrator) = error("get_tmp_cache!: method has not been implemented for the integrator")
user_cache(i::DEIntegrator) = error("user_cache: method has not been implemented for the integrator")
u_cache(i::DEIntegrator) = error("u_cache: method has not been implemented for the integrator")
du_cache(i::DEIntegrator) = error("du_cache: method has not been implemented for the integrator")
ratenoise_cache(i::DEIntegrator) = ()
rand_cache(i::DEIntegrator) = ()
full_cache(i::DEIntegrator) = error("full_cache: method has not been implemented for the integrator")
resize_non_user_cache!(i::DEIntegrator,ii::Int) = error("resize_non_user_cache!: method has not been implemented for the integrator")
deleteat_non_user_cache!(i::DEIntegrator,idxs) = error("deleteat_non_user_cache!: method has not been implemented for the integrator")
addat_non_user_cache!(i::DEIntegrator,idxs) = error("addat_non_user_cache!: method has not been implemented for the integrator")
terminate!(i::DEIntegrator) = error("terminate!: method has not been implemented for the integrator")
get_du(i::DEIntegrator) = error("get_du: method has not been implemented for the integrator")
get_du!(out,i::DEIntegrator) = error("get_du: method has not been implemented for the integrator")
get_dt(i::DEIntegrator) = error("get_dt: method has not been implemented for the integrator")
get_proposed_dt(i::DEIntegrator) = error("get_proposed_dt: method has not been implemented for the integrator")
set_proposed_dt!(i::DEIntegrator) = error("modify_proposed_dt!: method has not been implemented for the integrator")
"""
    savevalues!(integrator::DEIntegrator,
      force_save=false) -> Tuple{Bool, Bool}

Try to save the state and time variables at the current time point, or the
`saveat` point by using interpolation when appropriate. It returns a tuple that
is `(saved, savedexactly)`. If `savevalues!` saved value, then `saved` is true,
and if `savevalues!` saved at the current time point, then `savedexactly` is
true.

The saving priority/order is as follows:
  - `save_on`
    - `saveat`
    - `force_save`
    - `save_everystep`
"""
u_modified!(i::DEIntegrator,bool) = error("u_modified!: method has not been implemented for the integrator")
savevalues!(i::DEIntegrator) = error("savevalues!: method has not been implemented for the integrator")
add_tstop!(i::DEIntegrator,t) = error("add_tstop!: method has not been implemented for the integrator")
add_saveat!(i::DEIntegrator,t) = error("add_saveat!: method has not been implemented for the integrator")
set_abstol!(i::DEIntegrator,t) = error("set_abstol!: method has not been implemented for the integrator")
set_reltol!(i::DEIntegrator,t) = error("set_reltol!: method has not been implemented for the integrator")
reinit!(integrator::DEIntegrator,args...; kwargs...) =
       error("reinit!: method has not been implemented for the integrator")
auto_dt_reset!(integrator::DEIntegrator) = error("auto_dt_reset!: method has not been implemented for the integrator")
change_t_via_interpolation!(i::DEIntegrator,args...) = error("change_t_via_interpolation!: method has not been implemented for the integrator")
addsteps!(i::DEIntegrator,args...) = nothing
"""
    reeval_internals_due_to_modification!(integrator::DDEIntegrator)

Recalculate interpolation data and update ODE integrator after changes by callbacks.
"""
reeval_internals_due_to_modification!(integrator::DEIntegrator) = nothing

"""
    set_t!(integrator::DEIntegrator, t::Real)

Set current time point of the `integrator` to `t`.
"""
set_t!(integrator::DEIntegrator, t::Real) =
    error("set_t!: method has not been implemented for the integrator")

"""
    set_u!(integrator::DEIntegrator, u)

Set current state of the `integrator` to `u`.
"""
set_u!(integrator::DEIntegrator, u) =
    error("set_u!: method has not been implemented for the integrator")

"""
    set_ut!(integrator::DEIntegrator, u, t)

Set current state of the `integrator` to `u` and `t`
"""
function set_ut!(integrator::DEIntegrator, u, t)
  DiffEqBase.set_u!(integrator, u)
  DiffEqBase.set_t!(integrator, t)
end

### Addat isn't a real thing. Let's make it a real thing Gretchen

function addat!(a::AbstractArray,idxs,val=nothing)
  if val === nothing
    resize!(a,length(a)+length(idxs))
  else
    error("real addat! on arrays isn't supported yet")
    #=
    flip_range = last(idxs):-1:idxs.start
    @show idxs,flip_range
    splice!(a,flip_range,val)
    =#
  end
end

### Integrator traits

has_reinit(i::DEIntegrator) = false

### Display

Base.summary(I::DEIntegrator) = string(
                  TYPE_COLOR, nameof(typeof(I)),
                  NO_COLOR, " with uType ",
                  TYPE_COLOR, typeof(I.u),
                  NO_COLOR, " and tType ",
                  TYPE_COLOR, typeof(I.t), NO_COLOR)

function Base.show(io::IO, A::DEIntegrator)
  println(io,string("t: ",A.t))
  print(io,"u: ")
  show(io, A.u)
end
function Base.show(io::IO, m::MIME"text/plain", A::DEIntegrator)
  println(io,string("t: ",A.t))
  print(io,"u: ")
  show(io,m,A.u)
end
TreeViews.hastreeview(x::DiffEqBase.DEIntegrator) = true
function TreeViews.treelabel(io::IO,x::DiffEqBase.DEIntegrator,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Base.Text(Base.summary(x)))
end

### Error check (retcode)

last_step_failed(integrator::DEIntegrator) = false

"""
    check_error(integrator)

Check state of `integrator` and return one of the
[Return Codes](http://docs.juliadiffeq.org/latest/basics/solution.html#Return-Codes-(RetCodes)-1)
"""
function check_error(integrator::DEIntegrator)
  # This implementation is intended to be used for ODEIntegrator and
  # SDEIntegrator.
  if isnan(integrator.dt)
    if integrator.opts.verbose
      @warn("NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.")
    end
    return :DtNaN
  end
  if integrator.iter > integrator.opts.maxiters
    if integrator.opts.verbose
      @warn("Interrupted. Larger maxiters is needed.")
    end
    return :MaxIters
  end
  if !integrator.opts.force_dtmin && integrator.opts.adaptive &&
     abs(integrator.dt) <= abs(integrator.opts.dtmin) &&
     ((hasproperty(integrator,:opts) && hasproperty(integrator.opts,:tstops)) ?
     integrator.t + integrator.dt < integrator.tdir*top(integrator.opts.tstops) :
     true)
    if integrator.opts.verbose
      @warn("dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.")
    end
    return :DtLessThanMin
  end
  if integrator.opts.unstable_check(integrator.dt,integrator.u,integrator.p,integrator.t)
    if integrator.opts.verbose
      @warn("Instability detected. Aborting")
    end
    return :Unstable
  end
  if last_step_failed(integrator)
    if integrator.opts.verbose
      @warn("Newton steps could not converge and algorithm is not adaptive. Use a lower dt.")
    end
    return :ConvergenceFailure
  end
  return :Success
end

function postamble! end

"""
    check_error!(integrator)

Same as `check_error` but also set solution's return code
(`integrator.sol.retcode`) and run `postamble!`.
"""
function check_error!(integrator::DEIntegrator)
  code = check_error(integrator)
  if code != :Success
    integrator.sol = solution_new_retcode(integrator.sol, code)
    postamble!(integrator)
  end
  return code
end

### Default Iterator Interface
function done(integrator::DEIntegrator)
  if ! (integrator.sol.retcode in (:Default, :Success))
    return true
  elseif isempty(integrator.opts.tstops)
    postamble!(integrator)
    return true
  elseif integrator.just_hit_tstop
    integrator.just_hit_tstop = false
    if integrator.opts.stop_at_next_tstop
      postamble!(integrator)
      return true
    end
  end
  false
end
function Base.iterate(integrator::DEIntegrator,state=0)
  done(integrator) && return nothing
  state += 1
  step!(integrator) # Iter updated in the step! header
  # Next is callbacks -> iterator  -> top
  return integrator,state
end

Base.eltype(integrator::DEIntegrator) = typeof(integrator)

### Other Iterators

struct IntegratorTuples{I}
 integrator::I
end

function Base.iterate(tup::IntegratorTuples, state=0)
  done(tup.integrator) && return nothing
  step!(tup.integrator) # Iter updated in the step! header
  state += 1
  # Next is callbacks -> iterator  -> top
  return (tup.integrator.u,tup.integrator.t),state
end

RecursiveArrayTools.tuples(integrator::DEIntegrator) = IntegratorTuples(integrator)

"""
$(TYPEDEF)
"""
struct IntegratorIntervals{I}
 integrator::I
end

function Base.iterate(tup::IntegratorIntervals,state=0)
  done(tup.integrator) && return nothing
  state += 1
  step!(tup.integrator) # Iter updated in the step! header
  # Next is callbacks -> iterator  -> top
  return (tup.integrator.uprev,tup.integrator.tprev,tup.integrator.u,tup.integrator.t),state
end

intervals(integrator::DEIntegrator) = IntegratorIntervals(integrator)

struct TimeChoiceIterator{T,T2}
  integrator::T
  ts::T2
end

function Base.iterate(iter::TimeChoiceIterator,state=1)
  state > length(iter.ts) && return nothing
  t = iter.ts[state]
  integrator = iter.integrator
  if isinplace(integrator.sol.prob)
    tmp = first(get_tmp_cache(integrator))
    if t == integrator.t
      tmp .= integrator.u
    elseif t < integrator.t
      integrator(tmp,t)
    else
      step!(integrator,t-integrator.t)
      integrator(tmp,t)
    end
    return (tmp,t),state+1
  else
    if t == integrator.t
      tmp = integrator.u
    elseif t < integrator.t
      tmp = integrator(t)
    else
      step!(integrator,t-integrator.t)
      tmp = integrator(t)
    end
    return (tmp,t),state+1
  end
end

@recipe function f(integrator::DEIntegrator;
                    denseplot=(integrator.opts.calck || typeof(integrator) <: AbstractSDEIntegrator)  && integrator.iter>0,
                    plotdensity =10,
                    plot_analytic=false,vars=nothing)

  int_vars = interpret_vars(vars,integrator.sol)

  if denseplot
    # Generate the points from the plot from dense function
    plott = collect(range(integrator.tprev;step=integrator.t,length=plotdensity))
    plot_timeseries = integrator(plott)
    if plot_analytic
      plot_analytic_timeseries = [integrator.sol.prob.f(Val{:analytic},t,integrator.sol.prob.u0) for t in plott]
    end
  end # if not denseplot, we'll just get the values right from the integrator.

  dims = length(int_vars[1])
  for var in int_vars
    @assert length(var) == dims
  end
  # Should check that all have the same dims!

  plot_vecs = []
  for i in 2:dims
    push!(plot_vecs,[])
  end

  labels = String[]# Array{String, 2}(1, length(int_vars)*(1+plot_analytic))
  for x in int_vars
    for j in 2:dims
      if denseplot
        push!(plot_vecs[j-1], u_n(plot_timeseries, x[j],integrator.sol,plott,plot_timeseries))
      else # just get values
        if x[j] == 0
          push!(plot_vecs[j-1], integrator.t)
        elseif x[j]==1 && !(typeof(integrator.u) <: AbstractArray)
          push!(plot_vecs[j-1], integrator.u)
        else
          push!(plot_vecs[j-1], integrator.u[x[j]])
        end
      end
    end
    add_labels!(labels,x,dims,integrator.sol)
  end

  if plot_analytic
    for x in int_vars
      for j in 1:dims
        if denseplot
          push!(plot_vecs[j], u_n(plot_timeseries, x[j],sol,plott,plot_timeseries))
        else # Just get values
          if x[j] == 0
            push!(plot_vecs[j], integrator.t)
          elseif x[j]==1 && !(typeof(integrator.u) <: AbstractArray)
            push!(plot_vecs[j], integrator.sol.prob.f(Val{:analytic},integrator.t,integrator.sol[1]))
          else
            push!(plot_vecs[j], integrator.sol.prob.f(Val{:analytic},integrator.t,integrator.sol[1])[x[j]])
          end
        end
      end
      add_labels!(labels,x,dims,integrator.sol)
    end
  end

  xflip --> integrator.tdir < 0

  if denseplot
    seriestype --> :path
  else
    seriestype --> :scatter
  end

  # Special case labels when vars = (:x,:y,:z) or (:x) or [:x,:y] ...
  if typeof(vars) <: Tuple && (typeof(vars[1]) == Symbol && typeof(vars[2]) == Symbol)
    xlabel --> vars[1]
    ylabel --> vars[2]
    if length(vars) > 2
      zlabel --> vars[3]
    end
  end
  if getindex.(int_vars,1) == zeros(length(int_vars)) || getindex.(int_vars,2) == zeros(length(int_vars))
    xlabel --> "t"
  end

  linewidth --> 3
  #xtickfont --> font(11)
  #ytickfont --> font(11)
  #legendfont --> font(11)
  #guidefont  --> font(11)
  label --> reshape(labels,1,length(labels))
  (plot_vecs...,)
end

function step!(integ::DEIntegrator, dt::Real, stop_at_tdt = false)
    (dt * integ.tdir) < 0 && error("Cannot step backward.")
    t = integ.t
    next_t = t+dt
    stop_at_tdt && add_tstop!(integ,next_t)
    while integ.t*integ.tdir < next_t*integ.tdir
        step!(integ)
        integ.sol.retcode in (:Default, :Success) || break
    end
end

has_destats(i::DEIntegrator) = false


text:
using Logging
using JuMP, Gurobi
using DecisionProgramming

const N = 4

@info("Creating the influence diagram.")
diagram = InfluenceDiagram()

add_node!(diagram, ChanceNode("H1", [], ["ill", "healthy"]))
for i in 1:N-1
    # Testing result
    add_node!(diagram, ChanceNode("T$i", ["H$i"], ["positive", "negative"]))
    # Decision to treat
    add_node!(diagram, DecisionNode("D$i", ["T$i"], ["treat", "pass"]))
    # Cost of treatment
    add_node!(diagram, ValueNode("C$i", ["D$i"]))
    # Health of next period
    add_node!(diagram, ChanceNode("H$(i+1)", ["H$(i)", "D$(i)"], ["ill", "healthy"]))
end
add_node!(diagram, ValueNode("MP", ["H$N"]))

generate_arcs!(diagram)

# Add probabilities for node H1
add_probabilities!(diagram, "H1", [0.1, 0.9])

# Declare proability matrix for health nodes H_2, ... H_N-1, which have identical information sets and states
X_H = ProbabilityMatrix(diagram, "H2")
X_H["healthy", "pass", :] = [0.2, 0.8]
X_H["healthy", "treat", :] = [0.1, 0.9]
X_H["ill", "pass", :] = [0.9, 0.1]
X_H["ill", "treat", :] = [0.5, 0.5]

# Declare proability matrix for test result nodes T_1...T_N
X_T = ProbabilityMatrix(diagram, "T1")
X_T["ill", "positive"] = 0.8
X_T["ill", "negative"] = 0.2
X_T["healthy", "negative"] = 0.9
X_T["healthy", "positive"] = 0.1

for i in 1:N-1
    add_probabilities!(diagram, "T$i", X_T)
    add_probabilities!(diagram, "H$(i+1)", X_H)
end

for i in 1:N-1
    add_utilities!(diagram, "C$i", [-100.0, 0.0])
end

add_utilities!(diagram, "MP", [300.0, 1000.0])

generate_diagram!(diagram, positive_path_utility = true)

@info("Creating the decision model.")
model = Model()
z = DecisionVariables(model, diagram)
x_s = PathCompatibilityVariables(model, diagram, z, probability_cut = false)
EV = expected_value(model, diagram, x_s)
@objective(model, Max, EV)

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "IntFeasTol"      => 1e-9,
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)
S_probabilities = StateProbabilities(diagram, Z)
U_distribution = UtilityDistribution(diagram, Z)

@info("Printing decision strategy:")
print_decision_strategy(diagram, Z, S_probabilities)

@info("Printing utility distribution.")
print_utility_distribution(U_distribution)

@info("Printing statistics")
print_statistics(U_distribution)

@info("State probabilities:")
print_state_probabilities(diagram, S_probabilities, [["H$i" for i in 1:N]...])
print_state_probabilities(diagram, S_probabilities, [["T$i" for i in 1:N-1]...])
print_state_probabilities(diagram, S_probabilities, [["D$i" for i in 1:N-1]...])

@info("Conditional state probabilities")
for state in ["ill", "healthy"]
    S_probabilities2 = StateProbabilities(diagram, Z, "H1", state, S_probabilities)
    print_state_probabilities(diagram, S_probabilities2, [["H$i" for i in 1:N]...])
    print_state_probabilities(diagram, S_probabilities2, [["T$i" for i in 1:N-1]...])
    print_state_probabilities(diagram, S_probabilities2, [["D$i" for i in 1:N-1]...])
end


text:
using CSV, StateSpaceModels, Plots, Statistics, Dates

# Load the AirPassengers dataset
AP = CSV.read("AirPassengers.csv")

# Take the log of the series
logAP = log.(Vector{Float64}(AP.Passengers))

# Plot the data
p1 = plot(AP.Date, logAP, label = "Log-airline passengers", legend = :topleft, color = :black)

# Specify the state-space model
model = structural(logAP, 12)

# Estimate the state-space model
ss = statespace(model)

# Analyze its decomposition in trend, slope and seasonal components
p2 = plot(AP.Date, [ss.smoother.alpha[:, 1] ss.smoother.alpha[:, 2] ss.smoother.alpha[:, 3]], layout = (3, 1),
            label = ["Trend" "Slope" "Seasonal"], legend = :topleft)

# Forecast 24 months ahead
N = 24
pred, dist = forecast(ss, N)

# Define forecasting dates
firstdate = AP.Date[end] + Month(1)
newdates = collect(firstdate:Month(1):firstdate + Month(N - 1))

p3 = plot!(p1, newdates, pred, label = "Forecast")

# 10 scenarios for 24 months ahead
sim = simulate(ss, N, 10)
p4 = plot!(p1, newdates, sim[:, 1, :], label = "", color = "Grey", width = 0.2)


text:
# Run `julia --project=@. --color=yes make.jl` in this directory to
# build the documentation. The `Documenter` must be available (i.e.
# installed globally) for this to work.

using Documenter

using SimplexPoly

makedocs(; sitename="SimplexPoly", format=Documenter.HTML(; prettyurls=false))
deploydocs(; repo="github.com/eschnett/SimplexPoly.jl")


text:
using Random
"""
    toy_data_linear(N=100)

# Examples

```julia-repl
toy_data_linear()
```

"""
function toy_data_linear(N=100)
    # Number of points to generate.
    M = round(Int, N / 2)
    Random.seed!(1234)

    # Generate artificial data.
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i = 1:M])
    xt0s = Array([[x1s[i] - 5; x2s[i] - 5] for i = 1:M])

    # Store all the data for later.
    xs = [xt1s; xt0s]
    ts = [ones(M); zeros(M)];
    return xs, ts
end

using Random
"""
    toy_data_non_linear(N=100)

# Examples

```julia-repl
toy_data_non_linear()
```

"""
function toy_data_non_linear(N=100)
    # Number of points to generate.
    M = round(Int, N / 4)
    Random.seed!(1234)

    # Generate artificial data.
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt1s, Array([[x1s[i] - 5; x2s[i] - 5] for i = 1:M]))

    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt0s = Array([[x1s[i] + 0.5; x2s[i] - 5] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt0s, Array([[x1s[i] - 5; x2s[i] + 0.5] for i = 1:M]))

    # Store all the data for later.
    xs = [xt1s; xt0s]
    ts = [ones(2*M); zeros(2*M)];
    return xs, ts
end

using Random
"""
    toy_data_multi(N=100)

# Examples

```julia-repl
toy_data_multi()
```

"""
function toy_data_multi(N=100)
    # Number of points to generate.
    M = round(Int, N / 4)
    Random.seed!(1234)

    # Generate artificial data.
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt1s = Array([[x1s[i] + 1; x2s[i] + 1] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt1s, Array([[x1s[i] - 7; x2s[i] - 7] for i = 1:M]))

    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    xt0s = Array([[x1s[i] + 1; x2s[i] - 7] for i = 1:M])
    x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
    append!(xt0s, Array([[x1s[i] - 7; x2s[i] + 1] for i = 1:M]))

    # Store all the data for later.
    xs = [xt1s; xt0s]
    ts = [ones(M); ones(M).*2; ones(M).*3; ones(M).*4];
    return xs, ts
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ParseGen"
version = v"2.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/sandialabs/parsegen-cpp/archive/refs/tags/v$(version).tar.gz", "c6c7c4958d1c6ab77bf0970b5aacae4b63603f702666492638ee8c0bdf3125c8"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/10.15/MacOSX10.15.sdk.tar.xz",
    "2408d07df7f324d3beea818585a6d990ba99587c218a3969f924dfcc4de93b62")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/parsegen-cpp*/

install_license LICENSE

mkdir build && cd build

if [[ "${target}" == x86_64-apple-darwin* ]]; then
    export CXXFLAGS="-mmacosx-version-min=10.15"
    #install a newer SDK which supports `std::filesystem`
    pushd $WORKSPACE/srcdir/MacOSX10.*.sdk
    rm -rf /opt/${target}/${target}/sys-root/System
    cp -ra usr/* "/opt/${target}/${target}/sys-root/usr/."
    cp -ra System "/opt/${target}/${target}/sys-root/."
    popd
elif [[ "${target}" == aarch64-apple-darwin* ]]; then
    # TODO: we need to fix this in the compiler wrappers
    export CXXFLAGS="-mmacosx-version-min=11.0"
fi

cmake .. \
-DCMAKE_INSTALL_PREFIX=${prefix} \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DBUILD_SHARED_LIBS=ON

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental = true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libparsegen", :libparsegen)
    ExecutableProduct("parsegen-calc", :parsegen_calc)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
#this uses std filesystem, so we need gcc 8 at least
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"9")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
ZeroKernel([tr=IdentityTransform()])

Create a kernel always returning zero
"""
struct ZeroKernel{T,Tr} <: Kernel{T,Tr}
    transform::Tr
    metric::Delta
    function ZeroKernel{T,Tr}(t::Tr) where {T,Tr<:Transform}
        new{T,Tr}(t,Delta())
    end
end

function ZeroKernel(t::Tr=IdentityTransform()) where {Tr<:Transform}
    ZeroKernel{eltype(Tr),Tr}(t)
end

@inline kappa(κ::ZeroKernel,d::T) where {T<:Real} = zero(T)

"""
`WhiteKernel([tr=IdentityTransform()])`

```
    κ(x,y) = δ(x,y)
```
Kernel function working as an equivalent to add white noise.
"""
struct WhiteKernel{T,Tr} <: Kernel{T,Tr}
    transform::Tr
    metric::Delta
    function WhiteKernel{T,Tr}(t::Tr) where {T,Tr<:Transform}
        new{T,Tr}(t,Delta())
    end
end

function WhiteKernel(t::Tr=IdentityTransform()) where {Tr<:Transform}
    WhiteKernel{eltype(Tr),Tr}(t)
end

@inline kappa(κ::WhiteKernel,δₓₓ::Real) = δₓₓ

"""
`ConstantKernel([tr=IdentityTransform(),[c=1.0]])`
```
    κ(x,y) = c
```
Kernel function always returning a constant value `c`
"""
struct ConstantKernel{T,Tr,Tc<:Real} <: Kernel{T,Tr}
    transform::Tr
    metric::Delta
    c::Tc
    function ConstantKernel{T,Tr,Tc}(t::Tr,c::Tc) where {T,Tr<:Transform,Tc<:Real}
        new{T,Tr,Tc}(t,Delta(),c)
    end
end

params(k::ConstantKernel) = (params(k.transform),k.c)
opt_params(k::ConstantKernel) = (opt_params(k.transform),k.c)

function ConstantKernel(c::Tc=1.0) where {Tc<:Real}
    ConstantKernel{Float64,IdentityTransform,Tc}(IdentityTransform(),c)
end

function ConstantKernel(t::Tr,c::Tc=1.0) where {Tr<:Transform,Tc<:Real}
    ConstantKernel{eltype(Tr),Tr,Tc}(t,c)
end

@inline kappa(κ::ConstantKernel,x::Real) = κ.c


text:
using Flora.MLP
using MLDatasets # FashionMNIST

# load full test set
test_x,  test_y  = FashionMNIST.testdata() # 10_000

picks = [
    [0, 3, 9], # 0 T-Shirt, 3 Dress, 9 Ankle boot
    [1, 2, 5],
]

for pick in picks
    p = load(:fashion, pick)
    tX, tY = prepare(test_x, test_y, pick)
    show_accuracy(p.m, p.pick, tX, tY)
end

#=
p.pick = [0, 3, 9]
mean(argmax(p.m(x)) .== argmax(y)) = 0.956
p.pick = [1, 2, 5]
mean(argmax(p.m(x)) .== argmax(y)) = 0.989
=#

#=
julia> using UnicodePlots

julia> spy(reshape(tX[:,11], 28, 28)')
      Sparsity Pattern
      ┌──────────────┐
    1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ > 0
      │⠀⠀⠀⠀⠀⠀⢀⠇⣾⣷⣦⣤⣴⡄│ < 0
      │⠀⠀⠀⠀⢀⠀⠘⢰⣿⣿⣿⣿⣿⠟│
      │⢀⣀⣀⠠⠔⣪⣴⣿⣿⣿⣿⣿⣿⡇│
      │⢠⣶⣶⣾⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿│
      │⡛⠿⢿⣿⣿⣿⣿⣿⠿⠿⣿⣿⣿⣿│
   28 │⠈⠉⠀⠀⠀⠀⠀⠀⠈⠁⠀⠀⠀⠀│
      └──────────────┘
      1             28
         nz = 363
=#


text:
using StructuralCausalModels
using Test

include("test_dag_formulations.jl")
include("test_set_dag_df.jl")

include("test_dagitty_conversion.jl")
include("test_ggm_conversion.jl")

include("test_graphviz_conversions.jl")

include("test_basis_set_01.jl")

include("test_d_separation.jl")

include("test_descendents_01.jl")
include("test_descendents_02.jl")
include("test_descendents_03.jl")

include("test_open_paths_01.jl")
include("test_open_paths_02.jl")
include("test_open_paths_03.jl")
include("test_open_paths_04.jl")

include("test_ancestral_graph.jl")
include("test_ribbon_graph.jl")
include("test_maximize.jl")
#include("test_maximize_02.jl")
include("test_m_separation.jl")
include("test_m_separation_02.jl")
include("test_m_separation_03.jl")

include("test_sr6_4_2.jl")

include("test_sr6_4_3.jl")

include("test_sr6_4_3b.jl")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MPIR_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MPIR")
JLLWrappers.@generate_main_file("MPIR", UUID("426c3dee-e38f-5a7f-9125-e2686f44d6f1"))
end  # module MPIR_jll


text:
const GLOBAL_BUCKET = string(gensym(:global_bucket))
const RATE_LIMIT_RETRY = gensym(:retry)
const RATE_LIMIT_SENTINEL = gensym(:rate_limited)

@enum RateLimitHandler SENTINEL THROW WAIT

mutable struct Bucket
    limit::Int
    remaining::Int
    reset::DateTime

    Bucket() = new(typemax(Int), typemax(Int), typemax(DateTime))
end

struct RateLimiter
    buckets::Dict{String, Bucket}
    routes::Dict{String, String}
    handler::RateLimitHandler

    RateLimiter(handler=WAIT) = new(Dict(), Dict(), handler)
end

struct RateLimitedError <: Exception
    reset::DateTime
    response::Union{Response, Nothing}
end

RateLimitedError(reset) = RateLimitedError(reset, nothing)

on_rate_limit(rl::RateLimiter, reset, resp=nothing) =
    on_rate_limit(Val(rl.handler), reset, resp)
on_rate_limit(::Val{SENTINEL}, reset, resp) = RATE_LIMIT_SENTINEL
on_rate_limit(::Val{THROW}, reset, resp) = throw(RateLimitedError(reset, resp))
function on_rate_limit(::Val{WAIT}, reset, resp)
    interval = max(Millisecond(0), reset - now(UTC))
    @info "Rate limited, waiting for $interval"
    sleep(interval)
    return RATE_LIMIT_RETRY
end

function getbucket(rl::RateLimiter, path)
    return if haskey(rl.routes, path)
        rl.buckets[rl.routes[path]]
    else
        nothing
    end
end

getbucket!(rl::RateLimiter, key) = get!(Bucket, rl.buckets, key)
setpathkey!(rl::RateLimiter, path, key) = rl.routes[path] = key

function check_rate_limits(rl, path)
    bucket = getbucket(rl, path)
    if bucket !== nothing && bucket.remaining == 0 && now(UTC) < bucket.reset
        return on_rate_limit(rl, bucket.reset)
    end
end

function apply_rate_limits!(rl, resp)
    # TODO: I'm almost 100% sure that this isn't entirely correct.
    # The documentation is incredibly unclear.

    # Remove the /api/vN prefix.
    path = resp.request.target[8:end]
    is_global = header(resp, "X-RateLimit-Global") == "true"
    key = if is_global
        GLOBAL_BUCKET
    else
        bucket_key = header(resp, "X-RateLimit-Bucket")
        isempty(bucket_key) && return
        m = match(r"^/api/v\d/(?:channels|guilds|webhooks)/(\d+)", path)
        m === nothing ? bucket_key : "$bucket_key-$(m[1])"
    end

    bucket = getbucket!(rl, key)
    seconds = parse(Float64, header(resp, "X-RateLimit-Reset"))
    bucket.reset = unix2datetime(seconds)
    bucket.limit = parse(Int, header(resp, "X-RateLimit-Limit"))
    bucket.remaining = parse(Int, header(resp, "X-RateLimit-Remaining"))

    if key != GLOBAL_BUCKET
        setpathkey!(rl, path, key)
    end

    if resp.status == 429
        return on_rate_limit(rl, bucket.reset, resp)
    end
end


text:
println(5 + π)

println(typeof(π)) # Echo the type of pi
println(typeof(100)) # Echo the type of 100
println(typeof(6.0000)) # Echo the type of 6.0000

println(Sys.WORD_SIZE) # This variable indicated whether the target system is 32-bit arch or 64-bit arch

# %% typeof
hexX = 0x1
println("typeof(hexX) ~~> ", typeof(hexX))

hexX = 0x123
println("typeof(hexX) ~~> ", typeof(hexX))

hexX = 0x123456789
println("typeof(hexX) ~~> ", typeof(hexX))

x = 0b10
println("typeof(x) ~~> ", typeof(x))

y = 0o010
println("typeof(y) ~~> ", typeof(y))

for T in [Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128]
    println("* $(lpad(T,7)): [$(typemin(T)),$(typemax(T))]")
end
# %% End

# %% Overflow behavior
x = typemax(Int64)
println("x ~~> ", x)
x = x + 1
println("x ~~> ", x)
# x + 1 == typemax(Int64)
println(x + 1 == typemax(Int64))
# %% End

# %% Overflow resolution
println(10^19)
println(big(10)^19)
# %% End

# %% Float
myFloat = 1.777
println("typeof(myFLoat) ~~> ", typeof(myFloat))
println("typeof(Float32(myFloat)) ~~> ", typeof(Float32(myFloat)))
# %% End

# %% Positive and Negetive zero
println("0.0 == -0.0 ~~> ", 0.0 == -0.0)
println("bitstring(0.0) ~~> ", bitstring(0.0))
println("bitstring(-0.0) ~~> ", bitstring(-0.0))
# %% End

# %% Arbitrary Precision Arithmetic
println("BigInt(typemax(Int64)) + 1 ~~> ", BigInt(typemax(Int64)) + 1)
println((big"123456789012345678901234567890" + 1))
println(typeof(big"123456789012345678901234567890" + 1))
parse(BigInt, "123456789012345678901234567890")
println(string(big"2"^200, base = 16))
# %% End

# %% Numeric Literal Coefficients
x = 3
println("x ~~> ", x)
y = 2(x - 1)^2 - 3(x - 1) + 1
println("y ~~> ", y)
println("√16 ~~> ", √16)
# %% End


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

list of predefined matrix polynomials
=#

import DynamicPolynomials
const DP = DynamicPolynomials

function get_psdpoly_data(matpoly::Symbol)
    if matpoly == :matpoly1
        DP.@polyvar x
        M = [
            (x + 2x^3)  1;
            (-x^2 + 2)  (3x^2 - x + 1);
            ]
        MM = M' * M
        return ([x], MM, true)
    elseif matpoly == :matpoly2
        DP.@polyvar x
        poly = x^4 + 2x^2
        return ([x], poly, true)
    elseif matpoly == :matpoly3
        DP.@polyvar x y
        poly = (x + y)^4 + (x + y)^2
        return ([x, y], poly, true)
    elseif matpoly == :matpoly4
        n = 3
        m = 3
        d = 1
        DP.@polyvar x[1:n]
        Z = DP.monomials(x, 0:d)
        M = [sum(rand() * Z[l] for l in 1:length(Z)) for i in 1:m, j in 1:m]
        MM = M' * M
        MM = 0.5 * (MM + MM')
        return (x, MM, true)
    elseif matpoly == :matpoly5
        # example modified from
        # https://github.com/JuliaOpt/SumOfSquares.jl/blob/master/test/simplematrixsos.jl
        # Example 3.77 and 3.79 of Blekherman, G., Parrilo, P. A., & Thomas, R. R. (Eds.),
        # Semidefinite optimization and convex algebraic geometry SIAM 2013
        DP.@polyvar x
        P = [
            (x^2 - 2x + 2)  x;
            x               x^2;
            ]
        return ([x], P, true)
    elseif matpoly == :matpoly6
        # example modified from
        # https://github.com/JuliaOpt/SumOfSquares.jl/blob/master/test/choi.jl
        # verifies that a given polynomial matrix is not a Sum-of-Squares matrix
        # see Choi, M. D., "Positive semidefinite biquadratic forms",
        # Linear Algebra and its Applications, 1975, 12(2), 95-100
        DP.@polyvar x y z
        P = [
            (x^2 + 2y^2)    (-x * y)        (-x * z);
            (-x * y)        (y^2 + 2z^2)    (-y * z);
            (-x * z)        (-y * z)        (z^2 + 2x^2);
            ] .* (x * y * z)^0
        # TODO the (x * y * z)^0 can be removed when
        # https://github.com/JuliaOpt/SumOfSquares.jl/issues/106 is fixed
        return ([x, y, z], P, false)
    elseif matpoly == :matpoly7
        # example modified from
        # https://github.com/JuliaOpt/SumOfSquares.jl/blob/master/test/sosdemo9.jl
        # Section 3.9 of SOSTOOLS User's Manual, see
        # https://www.cds.caltech.edu/sostools/
        DP.@polyvar x y z
        P = hvcat((2, 2),
            (x^4 + x^2 * y^2 + x^2 * z^2),
            (x * y * z^2 - x^3 * y - x * y * (y^2 + 2 * z^2)),
            (x * y * z^2 - x^3 * y - x * y * (y^2 + 2 * z^2)),
            (x^2 * y^2 + y^2 * z^2 + (y^2 + 2 * z^2)^2))
        return ([x, y, z], P, true)
    end
end


text:
"""
    FM90(;c1=0.10, c2=0.70, c3=3.23, c4=0.41, x0=4.60, gamma=0.9)
    FM90(coeffs, x0=4.60, gamma=0.9)

Fitzpatrick & Massa (1990) generative model for ultraviolet dust extinction. The default values are the published values for the Milky Way average.

## Parameters
* `c1` - y-intercept of linear term
* `c2` - slope of liner term
* `c3` - amplitude of 2175 Å bump
* `c4` - amplitude of FUV rise
* `x0` - centroid of 2175 Å bump
* `gamma` - width of 2175 Å bump

If `λ` is a `Unitful.Quantity` it will be automatically converted to Å and the
returned value will be `UnitfulAstro.mag`.

## Examples
```jldoctest
julia> model = FM90(c1=0.2, c2=0.7, c3=3.23, c4=0.41, x0=4.6, gamma=0.99);

julia> model(1500)
5.2521258452800135

julia> FM90()(1500)
5.152125845280013

julia> FM90(c1=0.2, c2=0.7, c3=3.23, c4=0.41, x0=4.6, gamma=0.99).([1000, 1200, 1800])
3-element Array{Float64,1}:
 12.562237969522851
  7.769215017329513
  4.890128210972148

```
# Extended Help

The model has form ``c_1 + c_2x + c_3D(x; \\gamma, x_0) + c_4 F(x)`` where ``x`` is the wavenumber in inverse microns, ``D(x)``
is a Drude profile (modified Lorentzian) used to model the 2175 Å bump with the scale-free parameters ``x_0`` (central wavenumber)
and ``\\gamma`` (damping coefficient), and ``F(x)``, a piecewise function for the far-UV. Note that the coefficients will change
the overall normalization, possibly changing the expected behavior of reddening via the parameter ``A_V``.

## References
[Fitzpatrick & Massa (1990)](https://ui.adsabs.harvard.edu/abs/1990ApJS...72..163F)
"""
@with_kw struct FM90{T<:Number}  <: ExtinctionLaw @deftype T
    c1 = 0.10
    c2 = 0.70
    c3 = 3.23
    c4 = 0.41
    x0 = 4.60
    gamma = 0.99
    @assert x0 ≥ 0 "`x0` must be ≥ 0, got $x0"
    @assert gamma ≥ 0 "`gamma` must be ≥ 0, got $gamma"
end

FM90(c1, c2, c3, c4, x0, gamma) = FM90(promote(c1, c2, c3, c4, x0, gamma)...)
FM90(coeffs, x0=4.60, gamma=0.99) = FM90(coeffs..., x0, gamma)

bounds(::Type{<:FM90}) = (912, 3200)

function (law::FM90)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))

    x = aa_to_invum(wave)
    exvebv = law.c1 + law.c2 * x

    x2 = x^2
    exvebv += law.c3 * (x2 / ((x2 - law.x0^2)^2 + x2 * (law.gamma^2)))

    if x >= 5.9
        y = x - 5.9
        exvebv += law.c4 * @evalpoly y 0 0 0.5392 0.05644
    end

    return exvebv
end


text:
module Networks

using Printf, ..Lif, ..LifConfig

include("./Synapses.jl")
using .Synapses

export LIFNetwork, simulate, reset!, connect!

# by exporting the Synapses module, where 'using' is used to import
# this module, Synapses will be imported as if 'import' had been used
export Synapses
# ============================================================================ #
abstract type NeuralNetwork end

# ============================================================================ #
# LIFNetwork
# ============================================================================ #
mutable struct LIFNetwork{S<:BaseSynapse, N<:BaseNeuron} <: NeuralNetwork
    neurons::Array{N,1}
    synapses::Array{S,2}
end
# ---------------------------------------------------------------------------- #
function LIFNetwork(x::Array{T,1}) where {T<:BaseNeuron}
    n = length(x)
    LIFNetwork(x, Array{BaseSynapse,2}(undef, n, n))
end
# ---------------------------------------------------------------------------- #
function LIFNetwork(n::Integer, xi::Float64=0.2, typ::Type{T}=LIFNeuron) where {T<:BaseNeuron}
    x = Array{T,1}(undef, n)
    for k in eachindex(x)
        x[k] = typ()
        Lif.set!(x[k], :xi, xi)
    end
    LIFNetwork(x, Array{BaseSynapse,2}(undef, n, n))
end
# ---------------------------------------------------------------------------- #
function connect!(net::LIFNetwork)
    for k in 1:length(net)^2
        net.synapses[k] = Synapse()
    end
end
# # ---------------------------------------------------------------------------- #
# function connect!(net::LIFNetwork, args...)
#     d = Dict(args)
#     idx = keys(d)
#     n = length(net)
#
#     for r = 1:n
#         for c = 1:n
#             if (r,c) in idx
#                 net.synapses[r,c] = Synapse(d[(r,c)]...)
#             elseif !isassigned(net.synapses, r, c)
#                 net.synapses[r,c] = Synapse()
#             else
#                 @show(net.synapses[r,c])
#             end
#         end
#     end
# end
# ---------------------------------------------------------------------------- #
function connect!(net::LIFNetwork, k::Tuple{Integer, Integer}, p::Tuple)
    net.synapses[k...] = Synapse(p...)
end
# ---------------------------------------------------------------------------- #
Base.iterate(net::LIFNetwork) = (net[1], 1)
function Base.iterate(net::LIFNetwork, state::Integer)
    if state <= length(net)
        return (net[state], state+1)
    else
        return nothing
    end
end
Base.IteratorSize(net::LIFNetwork) = HasLength()
Base.IteratorEltype(net::LIFNetwork) = HasEltype()
Base.eltype(net::LIFNetwork) = eltype(net.neurons)
Base.length(net::LIFNetwork) = length(net.neurons)

Base.getindex(net::LIFNetwork, k::Integer) = net.neurons[k]
Base.getindex(net::LIFNetwork, k::UnitRange) = net.neurons[k]
Base.getindex(net::LIFNetwork, k::AbstractArray) = net.neurons[k]
# ---------------------------------------------------------------------------- #
function Base.show(io::IO, net::LIFNetwork)
    str = @sprintf("%dx1 Array{LIFNeuron,1}:\n", length(net.neurons))
    print(io, str)
    print(io, net.neurons)
    print(io, "\n\n")
    print(io, net.synapses)
end
# ---------------------------------------------------------------------------- #
function reset!(net::LIFNetwork)
    for x in net.neurons
        #NOTE: not sure why we have to specify the module here...
        Lif.reset!(x)
    end
    for s in net.synapses
        #NOTE: not sure why we have to specify the module here...
        Synapses.reset!(s)
    end
end
# ---------------------------------------------------------------------------- #
function collect_input(net::LIFNetwork, k::Integer, dt::Time, tnow::Time)
    totl = 0.0
    #input is received column-wise
    @inbounds for s in net.synapses[:,k]
        if isconnected(s)
            totl += receive!(net[k], s, dt, tnow)
        end
    end
    return totl
end
# ---------------------------------------------------------------------------- #
function distribute_output(net::LIFNetwork, k::Integer, spk::Bool, tnow::Time)
    #output is distributed row-wise
    @inbounds for s in net.synapses[k,:]
        send!(net[k], s, spk, tnow)
    end
end
# ---------------------------------------------------------------------------- #
# function update!(net::LIFNetwork, stim::Vector{Float64}, dt::Time, tnow::Time)
function update!(net::LIFNetwork, fstim::Function, dt::Time, tnow::Time)
    vm = zeros(length(net),1)
    spk = falses(length(net),1)

    @inbounds for k = 1:length(net)

        #get synaptic and stimulus input
        inp = collect_input(net, k, dt, tnow) + fstim(k, tnow)

        # update Vm
        vm[k], spk[k] = step!(net[k], inp, dt, tnow)

        #propagate output to synapses
        distribute_output(net, k, spk[k], tnow)

        if spk[k]
            vm[k] = net[k].p.vspike
        end
    end

    return vm, spk
end
# ============================================================================ #
# helper functions
# ============================================================================ #
# include("NetworksExt.jl")
# ============================================================================ #
function default_stim(id::Integer, t::Time)
    return 0.0
end
# ============================================================================ #
function simulate(net::T, fstim::Function, duration::Real=500.0,
    record::Vector{Int64}=Int64[]) where {T<:NeuralNetwork}

    tnow = 0.0
    dt = 0.043
    record_vm = false

    reset!(net)

    d = Dict{Symbol, AbstractArray}()

    times = 0.0:dt:duration
    if !isa(record, Array)
        record = [record]
    end
    if !isempty(record) && eltype(record) <: Integer
        filter!(x -> (x > 0 && x <= length(net)), record)
        d[:vm] = zeros(length(times),length(record))
        record_vm = true
    end
    d[:ts] = [Vector{Time}() for k in 1:length(net)]

    for k in eachindex(times)

        vm, spk = update!(net, fstim, dt, times[k])
        if record_vm
            d[:vm][k,:] = vm[record]
        end

        for j in 1:length(net)
            if spk[j]
                #NOTE: timestamps are stored in *SECONDS*
                push!(d[:ts][j], times[k]*1e-3)
            end
        end
    end

    d[:time] = times

    return d
end
# ============================================================================ #
# function test(n=2, xi=0.5)
#
#     net = LIFNetwork(n,xi)
#     for x in net
#         Lif.set!(x, :xi, 0.0)
#     end
#     connect!(net, (1,2)=>(Synapses.Static, .54, 3.0, 2.0))
#     stim(id, t) = t > 100.0 && id == 1 ? 0.2 : 0.0
#
#     d = simulate(net, fstim=stim, duration=500, record=[1])
#     return d
# end
# # ============================================================================ #

end #END MODULE


text:
export texture, textureMode, textureWrap, delTextures

function texture(img::Image)
	imagesep = separate(img)
	imagedata = reinterpret(Float64, map(Float64, data(imagesep)))

	# the following takes the image and turns into a 1D array that opengl
	# expects (i.e., the first 3 elements are the RGB values of the first
	# pixel, the next 3 elements are the RGB values of the second pixel, and
	# so on).
	# there must be an easier and faster way to do this with some sort of
	# opengl trick, though.
	if imagesep.properties["colordim"] == 3
		w = size(imagedata,2)
		h = size(imagedata,1)
		imageGL = Array(Float64, w*h*3)
		@inbounds for y=1:h
			row = 3 * w * (y - 1)
			for x=1:w
				col = 3 * (x - 1)
		        @simd for l=1:3
		            imageGL[row+col+l] = imagedata[y,x,l]
		        end
		    end
		end
	elseif imagesep.properties["colordim"] == 1
		w = size(imagedata,2)
		h = size(imagedata,3)
		imageGL = Array(Float64, w*h*3)
		@inbounds for y=1:h
			row = 3 * w * (y - 1)
			for x=1:w
				col = 3 * (x - 1)
		        @simd for l=1:3
		            imageGL[row+col+l] = imagedata[l,x,y]
		        end
		    end
		end
	end

	tex = GLuint[0]
	glGenTextures(1, tex)
	glActiveTexture(GL_TEXTURE2)
	glBindTexture(GL_TEXTURE_2D, tex[1])
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F, w, h, 0, GL_RGB, GL_FLOAT, map(Float32, imageGL))
	glGenerateMipmap(GL_TEXTURE_2D)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
	max_aniso = [Float32(0.0)]
	glGetFloatv(GL_TEXTURE_MAX_ANISOTROPY_EXT, max_aniso)
	# glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAX_ANISOTROPY_EXT, max_aniso[1])
	# this will require some performance testing. it might be too heavy for
	# psychophysics
	glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAX_ANISOTROPY_EXT, 16)
	glBindTexture(GL_TEXTURE_2D, C_NULL)

	return tex[1], w/state.width, h/state.height
end

function texture(data)
	w = size(data,2)
	h = size(data,1)

	tex = GLuint[0]
	glGenTextures(1, tex)
	glActiveTexture(GL_TEXTURE2)
	glBindTexture(GL_TEXTURE_2D, tex[1])
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F, w, h, 0, GL_RGB, GL_FLOAT, map(Float32, data[:]))
	glGenerateMipmap(GL_TEXTURE_2D)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
	max_aniso = [Float32(0.0)]
	glGetFloatv(GL_TEXTURE_MAX_ANISOTROPY_EXT, max_aniso)
	# glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAX_ANISOTROPY_EXT, max_aniso[1])
	# this will require some performance testing. it might be too heavy for
	# psychophysics
	glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAX_ANISOTROPY_EXT, 16)
	glBindTexture(GL_TEXTURE_2D, C_NULL)

	return tex[1], w/state.width, h/state.height
end

function textureMode(mode)

end

function textureWrap(wrap)
	if wrap == "CLAMP"
		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
	elseif wrap == "REPEAT"
		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT)
		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT)
	end
	return nothing
end

function delTextures(texs)
	glDeleteTextures(length(texs), collect(texs))
end


text:
for p in (:Sparams, :Yparams, :Zparams, :ABCDparams)
    @eval ^(param::($p){T}, N::Int) where {T<:Real} = ($p)(^(param.data,N))
    for f in (:+, :-, :*, :/)
        # Operation between NetworkParams
        @eval ($f)(p1::($p){T}, p2::($p){T}) where {T<:Real} =
            (p1.nPort == p2.nPort)? ($p)(($f)(p1.data, p2.data)) :
            error("The number of ports must be identical in order to perform binary operations")
            @eval ($f)(p1::($p){T}, p2::($p){S}) where {T<:Real, S<:Real} =
            (p1.nPort == p2.nPort)? ($p)(($f)(promote(p1.data, p2.data)...)):
            error("The number of ports must be identical in order to perform binary operations")
    end
end

"""
    permute_ports!(D::NetworkData{S, T}, I_before::Vector{Int},
        I_after::Vector{Int}) where {S<:Real, T<:NetworkParams}
Permute port indices
"""
function permute_ports!(D::NetworkData{S, T}, I_before::Vector{Int},
    I_after::Vector{Int}) where {S<:Real, T<:NetworkParams}
    if length(unique(I_before)) != length(I_before)
        error("Error: The indices contained in `I_before` must be unique")
    end
    if length(unique(I_after)) != length(I_after)
        error("Error: The indices contained in `I_after` must be unique")
    end
    if sort(I_before) != sort(I_after)
        error("Error: The arrays `I_before` and `I_after` must contain the same set of indices")
    end
    D.ports[I_after] = D.ports[I_before]
    for n in 1:D.nPoint
        # permute rows
        D.params[n].data[I_after, :] = D.params[n].data[I_before, :]
        # permute columns
        D.params[n].data[:, I_after] = D.params[n].data[:, I_before]
    end
    return D
end

permute_ports(D::NetworkData{S, T}, I_before::Vector{Int},
    I_after::Vector{Int}) where {S<:Real, T<:NetworkParams} =
    permute_ports!(deepcopy(D), I_before, I_after)

swap_ports!(D::NetworkData{S, T}, i1::Int, i2::Int) where {S<:Real,
    T<:NetworkParams} = permutePorts!(D, [i1, i2], [i2, i1])
swap_ports(D::NetworkData{S, T}, i1::Int, i2::Int) where {S<:Real,
    T<:NetworkParams} = swap_ports!(deepcopy(D), i1, i2)