text:
"""
    Limbdark

Fast, analytic, and numerically stable transit light curves for stars with
arbitrary order limb darkening.

"""
module Limbdark

using SpecialFunctions, LinearAlgebra

export transit_init, transit_poly, transit_poly!, transit_poly_g, transit_poly_g!
export integrate_lightcurve!

include("transit_structure.jl")
include("Mn_coeff.jl")
include("Nn_coeff.jl")
include("compute_g_n_struct.jl")
include("define_constants.jl")
include("cel_bulirsch.jl")
include("s2.jl")
include("Mn_compute.jl")
include("Nn_compute.jl")
include("transit_poly_struct.jl")
include("simpson_vec.jl")
include("integrate_transit_simpson_vec.jl")
include("integrate_lightcurve.jl")

end

max_issues_repo_name: rodluger/limbdark


text:
# Examples borrowed from: http://arxiv.org/pdf/1211.3759.pdf

# 2/3 N(x|0,1) + 1/3 N(x|0,1/10)
function kurtotic(theta)
    if theta[3] <= 0 || theta[6] <= 0
        return -Inf
    end
    llk = 0
    for i in 1:length(y)
      llk += log(theta[1] * pdf(Normal(theta[2], theta[3]), y[i]) +
                 theta[4] * pdf(Normal(theta[5], theta[6]), y[i]))
    end
    return llk
end

function bimodal_density(x)
    lik = (1/2) * pdf(Normal(-1,2/3),x) + (1/2) * pdf(Normal(1,2/3),x)
    return log(lik)
end

# 1/2 N(x|-1,2/3) + 1/2 N(x|1,2/3)
function bimodal_loglikelihood(theta)
    if theta[3] < 0 || theta[6] < 0
        return -Inf
    end
    llk = 0
    for i in 1:length(y)
      llk += log(theta[1] * pdf(Normal(theta[2], theta[3]), y[i]) +
                 theta[4] * pdf(Normal(theta[5], theta[6]), y[i]))
    end
    return -llk
end

# requires y and x
function logistic_regression(theta)
    llk = 0
    eta = [1 / (1 + exp(- x[i,:] * theta))[1] for i in 1:length(y)]
    for i = 1:length(y)
        llk += y[i] * log(eta[i]) + (1-y[i]) * log(1-eta[i])
    end
    return llk
end

ext: jl


text:
tests = ["BSplines.jl"];

for test in tests
    include(test)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE

using UMAT, Test, Materials, Libdl, Tensors

pkg_dir = dirname(Base.find_package("UMAT"))
if Sys.iswindows()
    lib_dir = joinpath(pkg_dir,"..","deps","usr","bin")
else
    lib_dir = joinpath(pkg_dir,"..","deps","usr","lib")
end

# E = 160e9; NU = 0.3; SYIELD = 100
PROPS = [160.e9, 0.3, 100.]
material = UmatMaterial(NTENS=6, NSTATV=0, NPROPS=3,
                        parameters=UmatParameterState(NPROPS=3,PROPS=PROPS),
                        lib_path=joinpath(lib_dir,"libmises_umat." * dlext),
                        behaviour=:mises_umat_)

dtime = 0.1
stresses = [copy(tovoigt(material.variables.stress))]
dstrain11 = 1e-3*dtime
dtimes = [dtime, dtime, dtime, dtime, 1.0]
dstrains11 = [dstrain11, dstrain11, dstrain11, -dstrain11, -4*dstrain11]

expected_S11 = [50. 100. 150. 100. -100.]
expected_stress = zeros(6)

expected_strains = [[0.00025, -7.5e-5, -7.5e-5, 0.0, 0.0, 0.0],
                    [0.0005, -0.00015, -0.00015, 0.0, 0.0, 0.0],
                    [0.00075, -0.000225, -0.000225, 0.0, 0.0, 0.0],
                    [0.0005, -0.00015, -0.00015, 0.0, 0.0, 0.0],
                    [-0.0005, 0.00015, 0.00015, 0.0, 0.0, 0.0]]

for i in 1:length(dtimes)
    dstrain11 = dstrains11[i]
    dtime = dtimes[i]
    uniaxial_increment!(material, dstrain11, dtime)
    update_material!(material)
    expected_stress[1] = expected_S11[i]
    @info(tovoigt(mat.variables.stress), stresses_expected[i])
    @info(tovoigt(material.drivers.strain; offdiagscale=2.0))
    #@test isapprox(tovoigt(material.variables.stress), expected_stress)
    #@test isapprox(tovoigt(material.drivers.strain; offdiagscale=2.0), expected_strains[i])
end


text:
function _build_ic_array_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # build the text fragments -
    ic_buffer=Array{String,1}()
    default_dictionary = ir_dictionary["default_dictionary"]
    system_type = ir_dictionary["system_type"]

    # go ...
    +(ic_buffer,"initial_condition_array = [";suffix="\n")

    # iterate through species, build ic records -
    df_species = ir_dictionary["model_species_table"]
    (number_of_species,_) = size(df_species)
    for species_index = 1:number_of_species

        species_symbol = df_species[species_index,:symbol]
        species_type = df_species[species_index, :type]

        # what is the value?
        value = 0.0
        if (species_type == :DNA)
            value = default_dictionary["parameter_dictionary"]["DNA_initial_condition"]
        elseif (species_type == :RNA)
            value = default_dictionary["parameter_dictionary"]["RNA_initial_condition"]
        elseif (species_type == :PROTEIN)
            value = default_dictionary["parameter_dictionary"]["PROTEIN_initial_condition"]
        end

        +(ic_buffer,"$(value)\t;\t#\t$(species_index)\t$(species_symbol)\tunits: nM\n"; prefix="\t\t\t")
    end

    # ok, if the system type is ANY type of CF_* then lets generate a translation
    # capacity term as an extra species -
    if (contains(system_type,"CF_") == true)
        +(ic_buffer,"\n")
        +(ic_buffer,"# translation capacity - \n"; prefix="\t\t\t")
        +(ic_buffer,"100.0\t;\t#\t$(number_of_species+1)\ttranslation capacity\tunits: NA\n"; prefix="\t\t\t")
    end

    +(ic_buffer,"]"; prefix="\t\t")

    # flatten and return -
    flat_ic_buffer = ""
    [flat_ic_buffer *= line for line in ic_buffer]
    return flat_ic_buffer
end

function _build_system_species_concentration_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # initialize -
    system_buffer = Array{String,1}()

    # get the system block -
    system_dictionary_array = ir_dictionary["system_dictionary"]["list_of_system_species"]

    # here we go - open
    +(system_buffer, "system_concentration_array = ["; suffix="\n")
    for (index,system_species_dictionary) in enumerate(system_dictionary_array)

        # ok, block will be a dictionary with a value, and units keys -
        value = system_species_dictionary["concentration"]
        units = system_species_dictionary["units"]
        symbol = system_species_dictionary["symbol"]

        # build the line -
        +(system_buffer, "$(value)\t;\t#\t$(symbol)\tunits: $(units)"; suffix="\n", prefix="\t\t\t")

    end

    # close -
    +(system_buffer,"]"; prefix="\t\t")

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in system_buffer]
    return flat_buffer
end

function _build_model_species_alias_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()

    # get the list of model species -
    df_species = ir_dictionary["model_species_table"]
    (number_of_species,_) = size(df_species)
    for species_index = 1:number_of_species

        # grab the symbol -
        species_symbol = df_species[species_index,:symbol]

        # this weird setup is because to the Mustache layout -
        if (species_index == 1)

            # build the record -
            +(buffer, "$(species_symbol) = x[$(species_index)]"; suffix="\n")
        elseif (species_index > 1 && species_index < number_of_species)

            # build the record -
            +(buffer, "$(species_symbol) = x[$(species_index)]"; prefix="\t", suffix="\n")
        else

            # build the record -
            +(buffer, "$(species_symbol) = x[$(species_index)]"; prefix="\t")
        end
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_system_type_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()

    # get the system block -
    system_type = ir_dictionary["system_dictionary"]["system_type"]

    # close -
    +(buffer,":$(system_type)"; prefix="\t")

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_system_species_alias_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()

    # get the system_dictionary -
    system_dictionary = ir_dictionary["system_dictionary"]
    list_of_system_species = system_dictionary["list_of_system_species"]
    for (index,species_dictionary) in enumerate(list_of_system_species)

        # get data -
        symbol = species_dictionary["symbol"]

        # formulate system species line -
        line = "$(symbol) = system_array[$(index)]"

        prefix="\t"
        if (index==1)
            prefix=""
        end

        # add to buffer -
        +(buffer, line; prefix=prefix, suffix="\n")
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_u_variable_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()

    # internal helper functions: handle activators -
    function _process_list_of_activators(target::String, output::String, polymerase_symbol::String,
        list_of_actor_dictionaries)::String

        # initialize -
        local_buffer = Array{String,1}()

        # generate list of W's -
        +(local_buffer,"W = ["; suffix="\n", prefix="\t")
        +(local_buffer,"W_$(target)\t;"; prefix="\t\t\t", suffix="\n")
        for (index, actor_dictionary) in enumerate(list_of_actor_dictionaries)

            # get actor symbol -
            actor_symbol = actor_dictionary["symbol"]

            # build the entry -
            +(local_buffer, "W_$(target)_$(actor_symbol)\t;"; prefix="\t\t\t", suffix="\n")
        end
        +(local_buffer,"]"; suffix="\n",prefix="\t")

        +(local_buffer, "$(target)_activator_array = Array{Float64,1}()"; suffix="\n", prefix="\t")
        +(local_buffer, "push!($(target)_activator_array, 1.0)"; suffix="\n", prefix="\t")
        for (_, actor_dictionary) in enumerate(list_of_actor_dictionaries)

            # get activator symbol -
            actor_symbol = actor_dictionary["symbol"]
            actor_type = actor_dictionary["type"]

            if (actor_type == "positive_sigma_factor")
                +(local_buffer, "$(actor_symbol)_$(polymerase_symbol) = $(polymerase_symbol)*f($(actor_symbol),K_$(target)_$(actor_symbol),n_$(target)_$(actor_symbol))"; suffix="\n", prefix="\t")
                +(local_buffer, "push!($(target)_activator_array, $(actor_symbol)_$(polymerase_symbol))"; suffix="\n", prefix="\t")
            end
        end

        # build the A term -
        +(local_buffer, "A = sum(W.*$(target)_activator_array)"; suffix="\n", prefix="\t")

        # flatten and return -
        flat_buffer = ""
        [flat_buffer *= line for line in local_buffer]
        return flat_buffer
    end

    function _process_list_of_repressors(target::String, output::String, polymerase_symbol::String,
        list_of_actor_dictionaries)::String

        # initialize -
        local_buffer = Array{String,1}()

        # do we have any repressors -
        if (isempty(list_of_actor_dictionaries) == true)
            +(local_buffer,"R = 0.0"; suffix="\n",prefix="\t")
        else
            # generate list of W's -
            +(local_buffer,"W = ["; suffix="\n", prefix="\t")
            for (index, actor_dictionary) in enumerate(list_of_actor_dictionaries)

                # get actor symbol -
                actor_symbol = actor_dictionary["symbol"]

                # build the entry -
                +(local_buffer, "W_$(target)_$(actor_symbol)\t;"; prefix="\t\t\t", suffix="\n")
            end
            +(local_buffer,"]"; suffix="\n",prefix="\t")

            # repressor -
            +(local_buffer, "$(target)_repressor_array = Array{Float64,1}()"; suffix="\n", prefix="\t")
            for (_, actor_dictionary) in enumerate(list_of_actor_dictionaries)

                # get activator symbol -
                actor_symbol = actor_dictionary["symbol"]
                actor_type = actor_dictionary["type"]
                counter_agent = actor_dictionary["counter_agent"]

                if (actor_type == "negative_protein_repressor")

                    +(local_buffer, "$(actor_symbol)_active = $(actor_symbol)*(1.0 - f($(counter_agent), K_$(target)_$(actor_symbol), n_$(target)_$(actor_symbol)))"; suffix="\n", prefix="\t")
                    +(local_buffer, "push!($(target)_repressor_array,$(actor_symbol)_active)"; prefix="\t", suffix="\n")
                else
                    +(local_buffer, "push!($(target)_repressor_array, f($(actor_symbol),K_$(target)_$(actor_symbol),n_$(target)_$(actor_symbol))"; suffix="\n", prefix="\t")
                end
            end

            # build the A term -
            +(local_buffer, "R = sum(W.*$(target)_repressor_array)"; suffix="\n", prefix="\t")
        end

        # flatten and return -
        flat_buffer = ""
        [flat_buffer *= line for line in local_buffer]
        return flat_buffer
    end

    # get the list of transcription models -
    list_of_transcription_models = ir_dictionary["list_of_transcription_models"]
    for (index, model_dictionary) in enumerate(list_of_transcription_models)

        # get data from the model dictionary -
        title = model_dictionary["title"]
        input = model_dictionary["input"]
        output = model_dictionary["output"]
        polymerase_symbol = model_dictionary["polymerase_symbol"]

        # process: list of activators -
        list_of_activators = model_dictionary["list_of_activators"]
        activator_clause = _process_list_of_activators(input, output, polymerase_symbol,list_of_activators)
        +(buffer,"\n";prefix="\t")
        +(buffer,"# - $(input) --------------------------------------------------------------------------------------------- ";suffix="\n", prefix="\t")
        +(buffer, "# $(input) activation - "; suffix="\n", prefix="\t")
        +(buffer, "$(activator_clause)"; suffix="\n")

        # process: list of repressors -
        list_of_repressors = model_dictionary["list_of_repressors"]
        repressor_clause = _process_list_of_repressors(input, output, polymerase_symbol, list_of_repressors)
        +(buffer, "# $(input) repression - "; suffix="\n", prefix="\t")
        +(buffer, "$(repressor_clause)")
        +(buffer, "push!(u_array, u(A,R))"; suffix="\n", prefix="\t")
        +(buffer,"\n")
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_transcription_kinetic_limit_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()

    # get the species table -
    model_species_table = ir_dictionary["model_species_table"]

    # get the list of transcription models -
    list_of_transcription_models = ir_dictionary["list_of_transcription_models"]
    for (index,transcription_model_dictionary) in enumerate(list_of_transcription_models)

        # get the input et al -
        input_string = transcription_model_dictionary["input"]
        polymerase_symbol = transcription_model_dictionary["polymerase_symbol"]

        # compute the length -
        idx_symbol = findfirst(x->x==input_string,model_species_table[!,:symbol])
        sequence = model_species_table[idx_symbol,:sequence]
        L = FASTX.FASTA.seqlen(sequence)

        # tmp_line =
        push_line = "push!(kinetic_limit_array, r(k_cat_characteristic,LX,$(L),$(polymerase_symbol),𝛕_$(input_string),K_$(input_string),$(input_string)))"

        if (index == 1)
            +(buffer,push_line; suffix="\n")
        elseif (index==length(list_of_transcription_models))
            +(buffer,push_line; prefix="\t")
        else
            +(buffer, push_line, prefix="\t",suffix="\n")
        end
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_translation_kinetic_limit_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()

    # get the species table -
    model_species_table = ir_dictionary["model_species_table"]

    # get the list of transcription models -
    list_of_models = ir_dictionary["list_of_translation_models"]
    for (index, model_dictionary) in enumerate(list_of_models)

        # get the input et al -
        output_string = model_dictionary["output"]
        input_string = model_dictionary["input"]
        polymerase_symbol = model_dictionary["ribosome_symbol"]

        # compute the length -
        idx_symbol = findfirst(x->x==output_string,model_species_table[!,:symbol])
        sequence = model_species_table[idx_symbol,:sequence]
        L = FASTX.FASTA.seqlen(sequence)

        # tmp_line =
        push_line = "push!(kinetic_limit_array, r(k_cat_characteristic,LL,$(L),$(polymerase_symbol),𝛕_$(output_string),K_$(output_string),$(input_string)))"

        if (index == 1)
            +(buffer,push_line; suffix="\n")
        elseif (index==length(list_of_models))
            +(buffer,push_line; prefix="\t")
        else
            +(buffer, push_line, prefix="\t",suffix="\n")
        end
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_model_parameter_array_snippet(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::NamedTuple

     # initialize -
    buffer = Array{String,1}()
    list_of_translation_models = ir_dictionary["list_of_translation_models"]
    list_of_transcription_models = ir_dictionary["list_of_transcription_models"]
    model_species_table = ir_dictionary["model_species_table"]
    parameter_symbol_array = Array{String,1}()

    # start -
    +(buffer,"model_parameter_array = ["; suffix="\n\n")

    # process the list of transcription models -
    pcounter = 1
    for (_, model_dictionary) in enumerate(list_of_transcription_models)

        # go through and formulate the W's -
        input_string = model_dictionary["input"]

        # grab the symbol -
        push!(parameter_symbol_array, "W_$(input_string)")

        # put a comment line -
        +(buffer,"# control parameters: $(input_string)\n";prefix="\t\t\t")
        +(buffer,"0.001\t;\t#\t$(pcounter)\tW_$(input_string)\n";prefix="\t\t\t")

        # update -
        pcounter = pcounter + 1

        # process the list of activators -
        list_of_activators = model_dictionary["list_of_activators"]
        for (_, activator_dictionary) in enumerate(list_of_activators)

            activator_symbol = activator_dictionary["symbol"]

            # grab the parameter symbol -
            push!(parameter_symbol_array, "W_$(input_string)_$(activator_symbol)")

            +(buffer,"1.0\t\t;\t#\t$(pcounter)\tW_$(input_string)_$(activator_symbol)\n"; prefix="\t\t\t")
            pcounter = pcounter + 1
        end

        # process the list of repressors -
        list_of_repressors = model_dictionary["list_of_repressors"]
        for (_, repressor_dictionary) in enumerate(list_of_repressors)

            repressor_symbol = repressor_dictionary["symbol"]

            # grab the parameter symbol -
            push!(parameter_symbol_array, "W_$(input_string)_$(repressor_symbol)")

            +(buffer,"1.0\t\t;\t#\t$(pcounter)\tW_$(input_string)_$(repressor_symbol)\n"; prefix="\t\t\t")
            pcounter = pcounter + 1
        end

        # binding constant, and binding order for activators -
        for (_, activator_dictionary) in enumerate(list_of_activators)

            activator_symbol = activator_dictionary["symbol"]

            # grab the parameter symbol -
            push!(parameter_symbol_array, "K_$(input_string)_$(activator_symbol)")
            push!(parameter_symbol_array, "n_$(input_string)_$(activator_symbol)")

            +(buffer,"1.0\t\t;\t#\t$(pcounter)\tK_$(input_string)_$(activator_symbol)\n"; prefix="\t\t\t")
            pcounter = pcounter + 1

            +(buffer,"1.0\t\t;\t#\t$(pcounter)\tn_$(input_string)_$(activator_symbol)\n"; prefix="\t\t\t")
            pcounter = pcounter + 1
        end

        # binding constant, and binding order for repressors -
        for (_, repressor_dictionary) in enumerate(list_of_repressors)

            repressor_symbol = repressor_dictionary["symbol"]
            +(buffer,"1.0\t\t;\t#\t$(pcounter)\tK_$(input_string)_$(repressor_symbol)\n"; prefix="\t\t\t")
            pcounter = pcounter + 1

            +(buffer,"1.0\t\t;\t#\t$(pcounter)\tn_$(input_string)_$(repressor_symbol)\n"; prefix="\t\t\t")
            pcounter = pcounter + 1

            # grab the parameter symbol -
            push!(parameter_symbol_array, "K_$(input_string)_$(repressor_symbol)")
            push!(parameter_symbol_array, "n_$(input_string)_$(repressor_symbol)")
        end

        # add extra row -
        +(buffer,"\n")
    end

    for (tx_index, model_dictionary) in enumerate(list_of_transcription_models)

        # what is the species -
        input_string = model_dictionary["input"]

        if (tx_index>1)
            +(buffer,"\n")
        end

        # K -
        +(buffer,"# transcription kinetic limit parameters: $(input_string)\n";prefix="\t\t\t")
        +(buffer, "1.0\t\t;\t#\t$(pcounter)\tK_$(input_string)\n"; prefix="\t\t\t")
        pcounter = pcounter + 1

        +(buffer, "1.0\t\t;\t#\t$(pcounter)\t𝛕_$(input_string)\n"; prefix="\t\t\t")
        pcounter = pcounter + 1

        # grab the parameter symbol -
        push!(parameter_symbol_array, "K_$(input_string)")
        push!(parameter_symbol_array, "𝛕_$(input_string)")
    end

    +(buffer,"\n")

    # process the list of translation models -
    for (tl_index, model_dictionary) in enumerate(list_of_translation_models)

        # go through and formulate the W's -
        output_string = model_dictionary["output"]

        if (tl_index>1)
            +(buffer,"\n")
        end

        # K -
        +(buffer,"# translation kinetic limit parameters: $(output_string)\n";prefix="\t\t\t")
        +(buffer, "1.0\t\t;\t#\t$(pcounter)\tK_$(output_string)\n"; prefix="\t\t\t")
        pcounter = pcounter + 1

        +(buffer, "1.0\t\t;\t#\t$(pcounter)\t𝛕_$(output_string)\n"; prefix="\t\t\t")
        pcounter = pcounter + 1

        # grab the parameter symbol -
        push!(parameter_symbol_array, "K_$(output_string)")
        push!(parameter_symbol_array, "𝛕_$(output_string)")
    end

    +(buffer,"\n")
    +(buffer,"# species degradation constants - \n"; prefix="\t\t\t")

    # ok, lastly we need degrdation constants for system species -
    (number_of_species, _) = size(model_species_table)
    for species_index = 1:number_of_species
        species_type = model_species_table[species_index, :type]
        species_symbol = model_species_table[species_index, :symbol]
        if (species_type != :DNA)
            +(buffer, "1.0\t\t;\t#\t$(pcounter)\t𝛳_$(species_symbol)\n"; prefix="\t\t\t")
            pcounter = pcounter + 1

            # grab the parameter symbol -
            push!(parameter_symbol_array, "𝛳_$(species_symbol)")
        end
    end

    # end -
    +(buffer, "]"; suffix="\n", prefix="\t\t")

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]

    results_tuple = (flat_buffer=flat_buffer, parameter_symbol_array=parameter_symbol_array)
    return results_tuple
end

function _build_model_parameter_symbol_index_map(parameter_symbol_array::Array{String,1})::String

    # initialize -
    buffer = Array{String,1}()

    # build the map -
    +(buffer,"model_parameter_symbol_index_map = Dict{Symbol,Int}()";suffix="\n")
    for (index, parameter_symbol) in enumerate(parameter_symbol_array)
        +(buffer,"model_parameter_symbol_index_map[:$(parameter_symbol)] = $(index)";prefix="\t\t",suffix="\n")
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_inverse_model_parameter_symbol_index_map(parameter_symbol_array::Array{String,1})::String

    # initialize -
    buffer = Array{String,1}()

    # build the map -
    +(buffer,"inverse_model_parameter_symbol_index_map = Dict{Int, Symbol}()";suffix="\n")
    for (index, parameter_symbol) in enumerate(parameter_symbol_array)
        +(buffer,"inverse_model_parameter_symbol_index_map[$(index)] = :$(parameter_symbol)";prefix="\t\t",suffix="\n")
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_parameter_alias_snippet(parameter_symbol_array::Array{String,1})::String

    # initialize -
    buffer = Array{String,1}()

    # build the alias -
    for (index, parameter_symbol) in enumerate(parameter_symbol_array)
        if (index == 1)
            +(buffer,"$(parameter_symbol) = model_parameter_array[model_parameter_index_map[:$(parameter_symbol)]]"; suffix="\n")
        else
            +(buffer,"$(parameter_symbol) = model_parameter_array[model_parameter_index_map[:$(parameter_symbol)]]"; prefix="\t", suffix="\n")
        end
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_degradation_dilution_snippet(model::VLJuliaModelObject,
    ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()
    model_species_table = ir_dictionary["model_species_table"]

    # iterate -
    (number_of_species, _) = size(model_species_table)
    for species_index = 1:number_of_species

        # get type and species -
        species_type = model_species_table[species_index, :type]
        species_symbol = model_species_table[species_index, :symbol]

        # if not DNA
        if (species_type != :DNA)

            prefix_string=""
            if (species_index>1)
                prefix_string="\t"
            end

            +(buffer, "push!(degradation_dilution_array, (μ + 𝛳_$(species_symbol))*$(species_symbol))"; prefix=prefix_string, suffix="\n")
        end
    end

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

function _build_system_dimension_block(model::VLJuliaModelObject,
    ir_dictionary::Dict{String,Any})::String

    # initialize -
    buffer = Array{String,1}()

    # dimensions info for the ir -
    number_of_transcription_models = ir_dictionary["number_of_transcription_models"]
    number_of_translation_models = ir_dictionary["number_of_translation_models"]

    # get transcription block -
    +(buffer, "\n"; prefix="\t")
    +(buffer,"number_of_transcription_processes = $(number_of_transcription_models)"; prefix="\t\t", suffix="\n")
    +(buffer,"number_of_translation_processes = $(number_of_translation_models)"; prefix="\t\t", suffix="\n")

    # flatten and return -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]
    return flat_buffer
end

# == MAIN METHODS BELOW HERE ======================================================================= #
function generate_data_dictionary_program_component(model::VLJuliaModelObject,
    ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    filename = "Parameters.jl"
    template_dictionary = Dict{String,Any}()

    try

        # what is the system type?
        system_type = ir_dictionary["system_type"]

        # small snippet -
        translation_half_life_block=""
        if (contains(system_type, "CF_") == true)
            translation_half_life_block*="# translation capacity half-life (default units: h^-1) -\n"
            translation_half_life_block*="\t\thalf_life_translation_capacity=value(biophysical_parameters_dictionary,:half_life_translation_capacity)\t# units: h^-1\n"
        else
            translation_half_life_block*="# translation capacity half-life (default units: h^-1) -\n"
            translation_half_life_block*="\t\thalf_life_translation_capacity=1000.0\t# units: h^-1\n"
        end

        # build the snippets required by the
        template_dictionary["copyright_header_text"] = build_julia_copyright_header_buffer(ir_dictionary)
        template_dictionary["initial_condition_array_block"] = _build_ic_array_snippet(model, ir_dictionary)
        template_dictionary["system_species_array_block"] = _build_system_species_concentration_snippet(model, ir_dictionary)
        template_dictionary["system_type_flag"] = _build_system_type_snippet(model, ir_dictionary)
        template_dictionary["system_dimension_block"] = _build_system_dimension_block(model, ir_dictionary)
        template_dictionary["translation_half_life_snippet"] = translation_half_life_block

        # we get the parameter array back in addition 2 the flat buffer for this method -
        results_tuple = _build_model_parameter_array_snippet(model, ir_dictionary)
        template_dictionary["model_parameter_array_block"] = results_tuple.flat_buffer

        # symbol index map -
        template_dictionary["model_parameter_symbol_index_map_block"] = _build_model_parameter_symbol_index_map(results_tuple.parameter_symbol_array)
        template_dictionary["inverse_model_parameter_symbol_index_map_block"] = _build_inverse_model_parameter_symbol_index_map(results_tuple.parameter_symbol_array)

        # write the template -
        template=mt"""
        {{copyright_header_text}}

        function update_parameter_dictionary(parameter_dictionary::Dict{String,Any}, varargs...)::Union{Nothing, Dict{String,Any}}

            # Default: return the unmodified dictionary -
            @info "Default implementation of update_parameter_dictionary returns unmodified dictionary"

            # return parameter dictionary -
            return parameter_dictionary
        end

        function generate_default_parameter_dictionary()::Dict{String,Any}

            # initialize -
            parameter_dictionary = Dict{String,Any}()
            system_type_flag = {{system_type_flag}}

            try

                # load the stoichiometric_matrix (SM) and degradation_dilution_matrix (DM) -
                SM = readdlm(joinpath(_PATH_TO_NETWORK,"Network.dat"))
                DM = readdlm(joinpath(_PATH_TO_NETWORK,"Degradation.dat"))

                # build the species initial condition array -
                {{initial_condition_array_block}}

                # setup the system dimension -
                {{system_dimension_block}}

                # build the system species concentration array -
                {{system_species_array_block}}

                # get the biophysical parameters for this system type -
                biophysical_parameters_dictionary = get_biophysical_parameters_dictionary_for_system(system_type_flag)

                # setup the model parameter array -
                {{model_parameter_array_block}}

                # setup the parameter symbol index map -
                {{model_parameter_symbol_index_map_block}}

                # setup the inverse parameter symbol index map -
                {{inverse_model_parameter_symbol_index_map_block}}

                # specific growth rate (default units: h^-1)
                μ = 0.0 # default units: h^-1

                {{translation_half_life_snippet}}

                # == DO NOT EDIT BELOW THIS LINE ========================================================== #
                parameter_dictionary["initial_condition_array"] = initial_condition_array
                parameter_dictionary["number_of_states"] = length(initial_condition_array)
                parameter_dictionary["number_of_transcription_processes"] = number_of_transcription_processes
                parameter_dictionary["number_of_translation_processes"] = number_of_translation_processes
                parameter_dictionary["system_concentration_array"] = system_concentration_array
                parameter_dictionary["biophysical_parameters_dictionary"] = biophysical_parameters_dictionary
                parameter_dictionary["model_parameter_array"] = model_parameter_array
                parameter_dictionary["model_parameter_symbol_index_map"] = model_parameter_symbol_index_map
                parameter_dictionary["inverse_model_parameter_symbol_index_map"] = inverse_model_parameter_symbol_index_map
                parameter_dictionary["stoichiometric_matrix"] = SM
                parameter_dictionary["dilution_degradation_matrix"] = DM
                parameter_dictionary["specific_growth_rate"] = μ
                parameter_dictionary["half_life_translation_capacity"] = half_life_translation_capacity

                # return -
                return parameter_dictionary
                # ========================================================================================= #
            catch error
                throw(error)
            end
        end
        """

        # render step -
        flat_buffer = render(template, template_dictionary)

        # package up into a NamedTuple -
        program_component = (buffer=flat_buffer, filename=filename, component_type=:buffer)

        # return -
        return program_component
    catch error
        rethrow(error)
    end
end

function generate_kinetics_program_component(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    filename = "Kinetics.jl"
    template_dictionary = Dict{String,Any}()

    try

        # we get the parameter array back in addition 2 the flat buffer for this method -
        results_tuple = _build_model_parameter_array_snippet(model, ir_dictionary)
        parameter_symbol_array = results_tuple.parameter_symbol_array

        # build snippets -
        template_dictionary["copyright_header_text"] = build_julia_copyright_header_buffer(ir_dictionary)
        template_dictionary["model_species_alias_block"] = _build_model_species_alias_snippet(model, ir_dictionary)
        template_dictionary["system_species_alias_block"] = _build_system_species_alias_snippet(model, ir_dictionary)
        template_dictionary["transcription_kinetic_limit_block"] = _build_transcription_kinetic_limit_snippet(model, ir_dictionary)
        template_dictionary["translation_kinetic_limit_block"] = _build_translation_kinetic_limit_snippet(model, ir_dictionary)
        template_dictionary["system_parameter_alias_block"] = _build_parameter_alias_snippet(parameter_symbol_array)
        template_dictionary["degradation_dilution_terms_block"] = _build_degradation_dilution_snippet(model, ir_dictionary)

        # setup the template -
        template = mt"""
        {{copyright_header_text}}

        function calculate_transcription_kinetic_limit_array(t::Float64, x::Array{Float64,1},
            parameter_dictionary::Dict{String,Any})::Array{Float64,1}

            # initialize -
            kinetic_limit_array = Array{Float64,1}()
            model_parameter_array = parameter_dictionary["model_parameter_array"]
            model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]
            system_array = parameter_dictionary["system_concentration_array"]
            eX = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["transcription_elongation_rate"].parameter_value)       # default units: nt/s
            LX = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["characteristic_transcript_length"].parameter_value)    # default units: nt
            k_cat_characteristic = (eX/LX)

            # helper function -
            r(kcat, L_char, L, polymerase, 𝛕, K, species) = kcat*(L_char/L)*polymerase*(species/(𝛕*K+(1+𝛕)*species))

            # alias the model species -
            {{model_species_alias_block}}

            # alias the system species -
            {{system_species_alias_block}}

            # alias the system parameters -
            {{system_parameter_alias_block}}

            # compute the transcription kinetic limit array -
            {{transcription_kinetic_limit_block}}

            # return -
            return kinetic_limit_array
        end

        function calculate_translation_kinetic_limit_array(t::Float64, x::Array{Float64,1},
            parameter_dictionary::Dict{String,Any})::Array{Float64,1}

            # initialize -
            kinetic_limit_array = Array{Float64,1}()
            system_array = parameter_dictionary["system_concentration_array"]
            model_parameter_array = parameter_dictionary["model_parameter_array"]
            model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]
            eL = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["translation_elongation_rate"].parameter_value)         # default units: aa/s
            LL = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["characteristic_protein_length"].parameter_value)       # default units: aa
            k_cat_characteristic = (eL/LL)

            # helper function -
            r(kcat, L_char, L, ribosome, 𝛕, K, species) = kcat*(L_char/L)*ribosome*(species/(𝛕*K+(1+𝛕)*species))

            # alias the model species -
            {{model_species_alias_block}}

            # alias the system species -
            {{system_species_alias_block}}

            # alias the system parameters -
            {{system_parameter_alias_block}}

            # compute the translation kinetic limit array -
            {{translation_kinetic_limit_block}}

            # return -
            return kinetic_limit_array
        end

        function calculate_dilution_degradation_array(t::Float64, x::Array{Float64,1},
            parameter_dictionary::Dict{String,Any})::Array{Float64,1}

            # initialize -
            μ = parameter_dictionary["specific_growth_rate"]
            degradation_dilution_array = Array{Float64,1}()
            model_parameter_array = parameter_dictionary["model_parameter_array"]
            model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]

            # alias the model species -
            {{model_species_alias_block}}

            # alias the system parameters -
            {{system_parameter_alias_block}}

            # formulate the degradation and dilution terms -
            {{degradation_dilution_terms_block}}

            # return -
            return degradation_dilution_array
        end
        """

        # render step -
        flat_buffer = render(template, template_dictionary)

        # package up into a NamedTuple -
        program_component = (buffer=flat_buffer, filename=filename, component_type=:buffer)

        # return -
        return program_component
    catch error
        rethrow(error)
    end
end

function generate_balances_program_component(model::VLJuliaModelObject,
    ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    filename = "Balances.jl"
    template_dictionary = Dict{String,Any}()

    try

        # do we need an extra state block -
        system_type = ir_dictionary["system_type"]
        extra_states_block = ""
        if (contains(system_type, "CF_") == true)
            number_of_species = ir_dictionary["number_of_species"]
            extra_states_block*="# extra species: global translation capacity -\n"
            extra_states_block*="\tdx[$(number_of_species+1)] = -(log(2)*(half_life_translation^-1))*x[$(number_of_species+1)]"
        end

        # build snippets -
        template_dictionary["copyright_header_text"] = build_julia_copyright_header_buffer(ir_dictionary)
        template_dictionary["extra_state_block"] = extra_states_block

        # setup the template -
        template = mt"""
        {{copyright_header_text}}
        function Balances(dx,x, parameter_dictionary,t)

            # system dimensions and structural matricies -
            number_of_states = parameter_dictionary["number_of_states"]
            DM = parameter_dictionary["dilution_degradation_matrix"]
            SM = parameter_dictionary["stoichiometric_matrix"]
            half_life_translation = parameter_dictionary["half_life_translation_capacity"]

            # calculate the TX and TL kinetic limit array -
            transcription_kinetic_limit_array = calculate_transcription_kinetic_limit_array(t,x,parameter_dictionary)
            translation_kinetic_limit_array = calculate_translation_kinetic_limit_array(t,x,parameter_dictionary)

            # calculate the TX and TL control array -
            u = calculate_transcription_control_array(t,x,parameter_dictionary)
            w = calculate_translation_control_array(t,x,parameter_dictionary)

            # calculate the rate of transcription and translation -
            rX = transcription_kinetic_limit_array.*u
            rL = translation_kinetic_limit_array.*w
            rV = [rX ; rL]

            # calculate the degradation and dilution rates -
            rd = calculate_dilution_degradation_array(t,x,parameter_dictionary)

            # compute the model equations -
            dxdt = SM*rV + DM*rd

            # package -
            for index = 1:number_of_states
                dx[index] = dxdt[index]
            end

            {{extra_state_block}}
        end
        """

        # render step -
        flat_buffer = render(template, template_dictionary)

        # package up into a NamedTuple -
        program_component = (buffer=flat_buffer, filename=filename, component_type=:buffer)

        # return -
        return program_component
    catch error
        rethrow(error)
    end
end

function generate_control_program_component(model::VLJuliaModelObject,
    ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    filename = "Control.jl"
    template_dictionary = Dict{String,Any}()

    try

        # what system type do we have?
        system_type = ir_dictionary["system_type"]
        number_of_species = ir_dictionary["number_of_species"]

        # we get the parameter array back in addition 2 the flat buffer for this method -
        results_tuple = _build_model_parameter_array_snippet(model, ir_dictionary)
        parameter_symbol_array = results_tuple.parameter_symbol_array

        # formulate the translation control block -
        translation_control_block_buffer=""
        if (contains(system_type,"CF_") == true)
            translation_control_block_buffer*="# translation decay: translation capacity decays over time -\n"
            translation_control_block_buffer*="\tcorrection_term = (x[$(number_of_species+1)]*(100.0^-1))\n"
            translation_control_block_buffer*="\tw = ones(number_of_translation_processes).*(correction_term)\n"
        else
            translation_control_block_buffer*="# default: w = 1 for all translation processes -\n"
            translation_control_block_buffer*="\tw = ones(number_of_translation_processes)\n"
        end

        # build snippets -
        template_dictionary["copyright_header_text"] = build_julia_copyright_header_buffer(ir_dictionary)
        template_dictionary["u_variable_snippet"] = _build_u_variable_snippet(model,ir_dictionary)
        template_dictionary["model_species_alias_block"] = _build_model_species_alias_snippet(model, ir_dictionary)
        template_dictionary["system_species_alias_block"] = _build_system_species_alias_snippet(model, ir_dictionary)
        template_dictionary["system_parameter_alias_block"] = _build_parameter_alias_snippet(parameter_symbol_array)
        template_dictionary["translation_control_variable_snippet"] = translation_control_block_buffer

        # build template -
        template=mt"""
        {{copyright_header_text}}

        # calculate the u-variables -
        function calculate_transcription_control_array(t::Float64, x::Array{Float64,1},
            parameter_dictionary::Dict{String,Any})::Array{Float64,1}

            # initialize -
            number_of_transcription_processes = parameter_dictionary["number_of_transcription_processes"]
            model_parameter_array = parameter_dictionary["model_parameter_array"]
            model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]
            u_array = Array{Float64,1}()

            # local helper functions -
            f(x,K,n) = (x^n)/(K^n+x^n)
            u(A,R) = A/(1+A+R)

            # alias the state vector -
            {{model_species_alias_block}}

            # alias the system species -
            system_array = parameter_dictionary["system_concentration_array"]
            {{system_species_alias_block}}

            # alias the system parameters -
            {{system_parameter_alias_block}}

            # == CONTROL LOGIC BELOW ================================================================= #
            {{u_variable_snippet}}
            # == CONTROL LOGIC ABOVE ================================================================= #

            # return -
            return u_array
        end

        # calculate the w-variables -
        function calculate_translation_control_array(t::Float64, x::Array{Float64,1},
            parameter_dictionary::Dict{String,Any})::Array{Float64,1}

            # initialize -
            number_of_translation_processes = parameter_dictionary["number_of_translation_processes"]

            # encode: encode the w-variables -
            {{translation_control_variable_snippet}}

            # return -
            return w
        end
        """

        # render step -
        flat_buffer = render(template, template_dictionary)

        # package up into a NamedTuple -
        program_component = (buffer=flat_buffer, filename=filename, component_type=:buffer)

        # return -
        return program_component
    catch error
        rethrow(error)
    end
end

function generate_init_program_component(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    filename = "Initialize.jl"
    template_dictionary = Dict{String,Any}()

    try

        # build snippets -
        template_dictionary["copyright_header_text"] = build_julia_copyright_header_buffer(ir_dictionary)

        # build template -
        template=mt"""
        {{copyright_header_text}}

        # stuff we need -
        import Pkg

        # initialization function -
        function __init__()

            # setup paths -
            _PATH_TO_ROOT = pwd()

            # get packages -
            Pkg.activate(_PATH_TO_ROOT)
            Pkg.add(name="DifferentialEquations")
            Pkg.add(name="VLModelParametersDB")
            Pkg.add(name="DelimitedFiles")
            Pkg.add(name="Plots")
            Pkg.instantiate()
        end

        # init me -
        __init__()
        """

        # render step -
        flat_buffer = render(template, template_dictionary)

        # package up into a NamedTuple -
        program_component = (buffer=flat_buffer, filename=filename, component_type=:buffer)

        # return -
        return program_component
    catch error
        rethrow(error)
    end
end

function generate_include_program_component(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    filename = "Include.jl"
    template_dictionary = Dict{String,Any}()

    try

        # build snippets -
        template_dictionary["copyright_header_text"] = build_julia_copyright_header_buffer(ir_dictionary)

        # build template -
        template=mt"""
        {{copyright_header_text}}

        # setup paths -
        const _PATH_TO_ROOT = pwd()
        const _PATH_TO_SRC = joinpath(_PATH_TO_ROOT, "src")
        const _PATH_TO_NETWORK = joinpath(_PATH_TO_SRC, "network")

        # use packages -
        using DifferentialEquations
        using VLModelParametersDB
        using DelimitedFiles
        using Plots

        # include my codes -
        include("./src/Parameters.jl")
        include("./src/Balances.jl")
        include("./src/Kinetics.jl")
        include("./src/Control.jl")
        """

        # render step -
        flat_buffer = render(template, template_dictionary)

        # package up into a NamedTuple -
        program_component = (buffer=flat_buffer, filename=filename, component_type=:buffer)

        # return -
        return program_component
    catch error
        rethrow(error)
    end
end

function generate_driver_program_component(model::VLJuliaModelObject, ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    filename = "Solve.jl"
    template_dictionary = Dict{String,Any}()

    try

        # build snippets -
        template_dictionary["copyright_header_text"] = build_julia_copyright_header_buffer(ir_dictionary)

        # build template -
        template=mt"""
        {{copyright_header_text}}
        # include -
        include("Include.jl")

        function solve_dynamic_problem(time_start::Float64, time_stop::Float64, time_step::Float64,
            parameter_dictionary::Dict{String,Any}, varargs...)

            # Get required stuff from the problem struct -
            time_span = (time_start,time_stop)
            initial_condition_array = parameter_dictionary["initial_condition_array"];

            # build problem object -
            problem_object = ODEProblem(Balances, initial_condition_array, time_span, parameter_dictionary)

            # solve -
            solution = solve(problem_object, AutoTsit5(Rosenbrock23(autodiff=false)), reltol=1e-8,abstol=1e-8)

            # pull solution apart -
            T = solution.t

            # initialize the state array -
            number_of_times_steps = length(T)
            number_of_states = length(initial_condition_array)
            X = zeros(number_of_times_steps,number_of_states)
            for step_index=1:number_of_times_steps

                # grab the solution
                soln_array = solution.u[step_index]
                for state_index = 1:number_of_states
                    X[step_index, state_index] = soln_array[state_index]
                end
            end

            # return -
            return (T,X)
        end

        function solve_static_problem(parameter_dictionary::Dict{String,Any}, varargs...)

            # Get the initial conditions -
            initial_condition_array = problem["initial_condition_array"];

            # calculate the steady-state soln -
            steady_state_prob = SteadyStateProblem(Balances, initial_condition_array, parameter_dictionary)
            steady_state_soln  = solve(steady_state_prob, DynamicSS(AutoTsit5(Rosenbrock23(autodiff=false)); abstol=1e-8, reltol=1e-6, tspan=Inf))

            # return -
            return abs.(steady_state_soln.u)
        end
        """

        # render step -
        flat_buffer = render(template, template_dictionary)

        # package up into a NamedTuple -
        program_component = (buffer=flat_buffer, filename=filename, component_type=:buffer)

        # return -
        return program_component
    catch error
        rethrow(error)
    end

end
# ================================================================================================= #"kinetic_limit_array


text:
using UnitRoot
using Base.Test

# write your own tests here
@test simul_arma()


text:
mutable struct LIF{VT<:Real, IT<:Integer} <: AbstractCell
    # required fields
    voltage::VT
    current::VT

    # model specific fields
    lastt::IT
    τm::VT
    vreset::VT
    R::VT
end

mutable struct ADEXP{VT<:Real, IT<:Integer} <: AbstractCell

    # {VFT=Vector{Float32},VBT=Vector{Bool}}  <: AbstractCell
    # required fields
    voltage::VT
    current::VT

    # model specific fields
    lastt::IT
    τm::VT
    vreset::VT
    R::VT
    #    param::ADEXParameter = ADEXParameter()
    #end
    #@with_kw mutable struct AD
    param::ADEXParameter = ADEXParameter(a,
                                        b,
                                        cm,
                                        v_rest,
                                        tau_m,
                                        tau_w,
                                        v_thresh,
                                        delta_T,
                                        v_reset,
                                        spike_height)
    N::SNNInt = 1
    cnt::SNNInt = 2
    v::Vector{SNNFloat} = fill(param.v_rest, N)
    w::Vector{SNNFloat} = zeros(N)
    fire::Vector{Bool} = zeros(Bool, N)
    I::Vector{SNNFloat} = zeros(N)
    spike_raster::Vector{SNNInt} = zeros(N)
    records::Dict = Dict()
end

@snn_kw struct ADEXParameter{FT=Float32}
    a::FT = 4.0
    b::FT = 0.0805
    cm::FT = 0.281
    v_rest::FT = -70.6
    tau_m::FT = 9.3667
    tau_w::FT = 144.0
    v_thresh::FT = -50.4
    delta_T::FT = 2.0
    v_reset::FT = -70.6
    spike_height::FT = 30
end

"""
LIF(τm, vreset, vth, R = 1.0)
Create a ADEXP neuron with zero initial voltage and empty current queue.
"""
LIF(τm::Real, vreset::Real, R::Real = 1.0) = LIF{Float32, Int}(vreset, 0, 0, τm, vreset, R)

isactive(neuron::LIF, t::Integer; dt::Real = 1.0) = (neuron.current > 0)
getvoltage(neuron::LIF) = neuron.voltage

"""
    excite!(neuron::LIF, current)
    excite!(neurons::AbstractArray{<:LIF}, current)
Excite a `neuron` with external `current`.
"""
excite!(neuron::LIF, current) = (neuron.current += current)
excite!(neurons::T, current) where T<:AbstractArray{<:LIF} = (neurons.current .+= current)

"""
    spike!(neuron::LIF, t::Integer; dt::Real = 1.0)
    spike!(neurons::AbstractArray{<:LIF}, spikes; dt::Real = 1.0)
Record a output spike from the threshold function with the `neuron` body.
Sets `neuron.lastt`.
"""
spike!(neuron::LIF, t::Integer; dt::Real = 1.0) = (t > 0) && (neuron.voltage = neuron.vreset)
function spike!(neurons::T, spikes; dt::Real = 1.0) where T<:AbstractArray{<:LIF}
    map!((x, y, s) -> (s > 0) ? y : x, neurons.voltage, neurons.voltage, neurons.vreset, spikes)
end

"""
    lif(t::Real, I, V; R, tau)
    lif!(t::AbstractArray{<:Real}, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; R::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    lif!(t::CuVector{<:Real}, I::CuVector{<:Real}, V::CuVector{<:Real}; R::CuVector{<:Real}, tau::CuVector{<:Real})
Evaluate a leaky integrate-and-fire neuron.
Use `CuVector` instead of `Vector` to evaluate on GPU.
# Fields
- `t`: time since last evaluation in seconds
- `I`: external current
- `V`: current membrane potential
- `vrest`: resting membrane potential
- `R`: resistance constant
- `tau`: time constant
"""
function lif(t::Real, I, V; vrest, R, tau)
    V = V * exp(-t / tau) + vrest
    V += I * (R / tau)

    return V
end
function lif!(t::AbstractArray{<:Real}, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; vrest::AbstractArray{<:Real}, R::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    # account for leakage
    @. V = V * exp(-t / tau) + vrest

    # apply update step
    @. V += I * (R / tau)

    return V
end
function lif!(t::CuVector{<:Real}, I::CuVector{<:Real}, V::CuVector{<:Real}; vrest::CuVector{<:Real}, R::CuVector{<:Real}, tau::CuVector{<:Real})
    # account for leakage
    @. V = V * exp(-t / tau) + vrest

    # apply update step
    @. V += I * (R / tau)

    return V
end

"""
    evaluate!(neuron::LIF, t::Integer; dt::Real = 1.0)
    (neuron::LIF)(t::Integer; dt::Real = 1.0)
    evaluate!(neurons::AbstractArray{<:LIF}, t::Integer; dt::Real = 1.0)
Evaluate the neuron model at time `t`. Return the resulting membrane potential.
"""
function evaluate!(neuron::LIF, t::Integer; dt::Real = 1.0)
    neuron.voltage = lif((t - neuron.lastt) * dt, neuron.current, neuron.voltage; vrest = 0, R = neuron.R, tau = neuron.τm)
    neuron.lastt = t
    neuron.current = 0

    return neuron.voltage
end
(neuron::LIF)(t::Integer; dt::Real = 1.0) = evaluate!(neuron, t; dt = dt)
function evaluate!(neurons::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:LIF}
    lif!((t .- neurons.lastt) .* dt, neurons.current, neurons.voltage;
         vrest = zero(neurons.voltage), R = neurons.R, tau = neurons.τm)
    neurons.lastt .= t
    neurons.current .= 0

    return neurons.voltage
end

"""
    reset!(neuron::LIF)
    reset!(neurons::AbstractArray{<:LIF})
Reset `neuron` by setting the membrane potential to `neuron.vreset`.
"""
function reset!(neuron::LIF)
    neuron.lastt = 0
    neuron.voltage = neuron.vreset
    neuron.current = 0
end
function reset!(neurons::T) where T<:AbstractArray{<:LIF}
    neurons.lastt .= 0
    neurons.voltage .= neurons.vreset
    neurons.current .= 0
end

Base.show(io::IO, ::MIME"text/plain", neuron::LIF) =
    print(io, """LIF:
                     voltage: $(neuron.voltage)
                     current: $(neuron.current)
                     τm:      $(neuron.τm)
                     vreset:  $(neuron.vreset)
                     R:       $(neuron.R)
    """)
Base.show(io::IO, neuron::LIF) =
    print(io, "LIF(voltage = $(neuron.voltage), current = $(neuron.current))")


text:
## test_fem.jl : test GRF defined on FEM

@testset "finite element mesh      " begin

## Load node/element table ##
(p,t) = star()
m = Matern(0.2,2.0)

## Cholesky ##
grf = GaussianRandomField(CovarianceFunction(2,m),Cholesky(),p,t)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,Matern)
@test ndims(grf.cov) == 2
@test isa(grf,GaussianRandomField{C,Cholesky} where {C})
@test length(grf.pts) == 2
@test size(grf.pts[1],1) == 2
@test size(grf.pts[1],2) == size(p,1)
@test size(grf.pts[2],1) == 3
@test size(grf.pts[2],2) == size(t,1)
#plot_trisurf(grf); close()
#tricontourf(grf); close()
#plot(grf); close()

## Spectral ##
grf = GaussianRandomField(CovarianceFunction(2,m),Spectral(),p,t)
@test isa(grf,GaussianRandomField{C,Spectral} where {C})

## CirculantEmbedding ##
@test_throws ArgumentError GaussianRandomField(CovarianceFunction(2,m),CirculantEmbedding(),p,t)

## KarhunenLoeve ##
grf = GaussianRandomField(CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="nodes")
@test isa(grf,GaussianRandomField{C,KarhunenLoeve{100}} where {C})
#plot_eigenvalues(grf); close()
#plot_eigenfunction(grf,1); close()

## Anisotropic random field ##
a = AnisotropicExponential([1000 0; 0 1000])
cov = CovarianceFunction(2,a)
grf = GaussianRandomField(cov,KarhunenLoeve(500),p,t,quad=GaussLegendre())
@test isa(grf.cov.cov,AnisotropicExponential)

## Center mode ##
grf = GaussianRandomField(CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="center")
@test isa(grf,GaussianRandomField{C,KarhunenLoeve{100}} where {C})
@test_throws ArgumentError GaussianRandomField(CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="cente")
@test_throws ArgumentError GaussianRandomField(1,CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="cente")
@test_throws ArgumentError GaussianRandomField(zeros(size(t,1)),CovarianceFunction(2,m),KarhunenLoeve(100),p,t,mode="cente")

## GaussianRandomField with mean value ##
grf = GaussianRandomField(1,CovarianceFunction(2,m),Cholesky(),p,t,mode="center")
@test isa(grf,GaussianRandomField{C,Cholesky} where {C})
@test length(grf.mean) == size(t,1)
@test all(grf.mean.==1)

grf = GaussianRandomField(ones(size(p,1)),CovarianceFunction(2,m),Spectral(),p,t)
@test isa(grf,GaussianRandomField{C,Spectral} where {C})
@test length(grf.mean) == size(p,1)
@test all(grf.mean.==1)

## Non-bounding box method ##
(p,t) = Lshape()
grf = grf = GaussianRandomField(CovarianceFunction(2,Matern(0.2,1.0)),Spectral(),p,t,n=10)
@test size(grf.data.eigenfunc,2) == 10

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

type Coupling <: BoundaryProblem
    reference_node :: Element{Poi1}
end

function Coupling()
    reference_node = Element(Poi1, Int64[])
    return Coupling(reference_node)
end

function add_coupling_nodes!(problem, elements)
    for element in elements
        push!(problem.elements, element)
    end
    return nothing
end

function add_reference_node!(problem, element)
    problem.properties.reference_node = element
    return nothing
end

function FEMBase.assemble_elements!(problem::Problem{Coupling},
                                    assembly::Assembly,
                                    elements::Vector{Element{Poi1}}, time)

    dimensions = problem.dimension

    function to3d(x)
        if length(x) == 2
            return [x; 0.0]
        end
        return x
    end

    ref_node = problem.properties.reference_node
    if length(get_connectivity(ref_node)) == 0
        error("Reference node node defined. Define reference node using add_reference_node!")
    end
    ref_node_id = first(get_connectivity(ref_node))
    X_r = to3d(first(ref_node("geometry", time)))
    fe = zeros(dimensions)
    weights = Dict{Int64,Float64}()
    X_n = Dict{Int64,Vector{Float64}}()
    for coupling_node in elements
        coupling_node_id = first(get_connectivity(coupling_node))
        X_n[coupling_node_id] = to3d(first(coupling_node("geometry", time)))
        weights[coupling_node_id] = 1.0
    end
    total = 0.0
    for n in keys(weights)
        total = total+weights[n]
    end
    for n in keys(weights)
        weights[n] = weights[n]/total
    end
    xbar = zeros(3)
    for n in keys(weights)
        xbar = xbar+weights[n]*X_n[n]
    end
    r = Dict{Int64,Vector{Float64}}()
    for n in keys(weights)
        r[n] = X_n[n]-xbar
    end
    rR = X_r-xbar
    T = zeros(3,3)
    for n in keys(weights)
        T += weights[n]*(dot(r[n],r[n])*eye(3)-(r[n]*r[n]'))
    end
    if det(T) == 0
        indz = indmin(diag(T))
        T[indz,indz] += 1e-9
    end

    invT = inv(T)

    for coupling_node in elements
        fill!(fe, 0.0)
        coupling_node_id = first(get_connectivity(coupling_node))
        X_n = first(coupling_node("geometry", time))
        gdofs = get_gdofs(problem, coupling_node)
        FR = zeros(3)
        MR = zeros(3)
        for i = 1:3
            if haskey(ref_node, "point moment $i")
                MR[i] = ref_node("point moment $i", time)
            end
            if haskey(ref_node, "point force $i")
                FR[i] = ref_node("point force $i", time)
            end
        end

        MRhat = MR + cross(rR,FR)
        n = coupling_node_id
        Fn = weights[n]*(FR+cross(invT*MRhat,r[n]))
        fe = fe+Fn[1:dimensions]
        add!(assembly.f, gdofs, fe)
    end

    return nothing

end


text:
import Test
import PredictMD

a = PredictMD.version()
Test.@test( typeof(a) == VersionNumber )
Test.@test( typeof(a) === VersionNumber )
Test.@test( a != VersionNumber(0) )
Test.@test( a > VersionNumber(0) )
Test.@test( a > VersionNumber("0.1.0") )
Test.@test( a < VersionNumber("123456789.0.0") )


text:
@testset "Orderings" begin

    import KnuthBendix.set_inversion!
    import Base.Order.lt

    A = Alphabet(['a', 'b', 'c', 'd'])
    set_inversion!(A, 'a', 'b')
    set_inversion!(A, 'c', 'd')

    lenlexord = LenLex(A)

    @test lenlexord isa Base.Order.Ordering

    u1 = Word([1,2])
    u3 = Word([1,3])
    u4 = Word([1,2,3])
    u5 = Word([1,4,2])

    @test lt(lenlexord, u1, u3) == true
    @test lt(lenlexord, u3, u1) == false
    @test lt(lenlexord, u3, u4) == true
    @test lt(lenlexord, u4, u5) == true
    @test lt(lenlexord, u5, u4) == false
    @test lt(lenlexord, u1, u1) == false

    wo = WreathOrder(A)
    @test wo isa KnuthBendix.WordOrdering

    w1 = Word([3,1,4,2,3])
    w5 = Word([1,2,3,2,4,1,2,4])
    w4 = Word([1,3,1,2,4,2,1])
    w6 = Word([1,2,3,2,1,4,2,4])
    w2 = Word([1,3,4,3,2])
    w3 = Word([1,3,2,1,4,1,2])
    ε = one(w1)

    a = [w1, w5, w4, w6, w2, w3, ε]

    @test sort(a, order = wo) == [ε, w1, w2, w3, w4, w5, w6]

    rpo = RecursivePathOrder(A)

    w1 = Word([1])
    w2 = Word([2])
    w14 = Word([1,4])
    w214 = Word([2,1,4])
    w41 = Word([4,1])
    w241 = Word([2,4,1])
    w1224 = Word([1,2,3,2,1,4,2,4])
    w32141 = Word([3,2,1,4,1])

    a = [w14, w1, w1224, w214, ε]
    b = [w32141, w214, w1, ε, w241, w2, w41]

    @test sort(a, order = rpo) == [ε, w1, w14, w214, w1224]
    @test sort(b, order = rpo) == [ε, w1, w2, w214, w41, w241, w32141]

end


text:
# Conversion of rotation matrix to quaternion
# ===========================================
# References
# ----------
# [1] http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche52.html
# --------------------------------------------------------------------------
# rotation matrix R is passed as input; quaternion q is returned
#_______________________________________________________________

function rotqua(R)

  q0 = ( R[1,1] + R[2,2] + R[3,3] + 1.0) / 4.0
  q1 = ( R[1,1] - R[2,2] - R[3,3] + 1.0) / 4.0
  q2 = (-R[1,1] + R[2,2] - R[3,3] + 1.0) / 4.0
  q3 = (-R[1,1] - R[2,2] + R[3,3] + 1.0) / 4.0

  if q0 < 0.0 q0 = 0.0 end
  if q1 < 0.0 q1 = 0.0 end
  if q2 < 0.0 q2 = 0.0 end
  if q3 < 0.0 q3 = 0.0 end

  q0 = sqrt(q0)
  q1 = sqrt(q1)
  q2 = sqrt(q2)
  q3 = sqrt(q3)

  if q0 >= q1 && q0 >= q2 && q0 >= q3
      q0 *= 1.0
      q1 *= sign(R[3,2] - R[2,3])
      q2 *= sign(R[1,3] - R[3,1])
      q3 *= sign(R[2,1] - R[1,2])
  elseif q1 >= q0 && q1 >= q2 && q1 >= q3
      q0 *= sign(R[3,2] - R[2,3])
      q1 *= 1.0
      q2 *= sign(R[2,1] + R[1,2])
      q3 *= sign(R[1,3] + R[3,1])
  elseif q2 >= q0 && q2 >= q1 && q2 >= q3
      q0 *= sign(R[1,3] - R[3,1])
      q1 *= sign(R[2,1] + R[1,2])
      q2 *= 1.0
      q3 *= sign(R[3,2] + R[2,3])
  elseif q3 >= q0 && q3 >= q1 && q3 >= q2
      q0 *= sign(R[2,1] - R[1,2])
      q1 *= sign(R[3,1] + R[1,3])
      q2 *= sign(R[3,2] + R[2,3])
      q3 *= 1.0
  else
      println("rotation_to_quaternion coding error")
  end

  r = sqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3);

  q0 /= r
  q1 /= r
  q2 /= r
  q3 /= r

  return [q0; q1; q2; q3]

end

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
"""
nsolheq(x0, FS, FPS, hdata; diff=:fd)
Internal function to run with CI. Nothing to see here, move along.
"""
function nsolheq(x0, FS, FPS, hdata; diff = :fd)
    if diff == :fd
        heqout = nsol(heqf!, x0, FS, FPS; pdata = hdata, sham = 1)
    else
        heqout = nsol(heqf!, x0, FS, FPS, heqJ!; pdata = hdata, sham = 1)
    end
    return heqout
end


text:
include("addnoise.jl")

import ...Embeddings.AbstractEmbedding
import DelayEmbeddings.Dataset
import Simplices.Delaunay.delaunay
import CausalityToolsBase: CustomReconstruction
import StaticArrays:
    SArray, MArray, SVector, MVector

"""
    DelaunayTriangulation

Contains the vertex indices of the simplices furnishing a Delaunay triangulation.

## Indexing

If `dt` is a `DelaunayTriangulation`, then

- ``dt[1]`` returns the indices of the 1st simplex.
- ``dt[1:4]`` returns a vector of indices for the 1st through 4th simplices.
- ``dt[[5, 8]]`` returns a vector of indices for the 5th and 8th simplices.

Calling `dt[:, i]` will return the same as with single-dimension indexing, but
converts the indices into an
array where each column is a set of simplex vertex indices, i.e

``dt[:, 1]`` returns the indices of the 1st simplex.
"""
struct DelaunayTriangulation <: AbstractDelaunayTriangulation
    indices::Vector{Vector{Int32}}
end

dimension(DT::DelaunayTriangulation) = length(DT.indices[1]) - 1

function DelaunayTriangulation(vertices::Vector{Vector{T}}; joggle = 0.001) where T

    # Slightly joggle points to avoid problems with QHull
    if joggle > 0
        addnoise!(vertices; joggle_factor = joggle)
    end

    DelaunayTriangulation(delaunay(vertices))
end

"""
DelaunayTriangulation(vertices; joggle) where {T} -> DelaunayTriangulation

Construct a Delaunay triangulation from a set of points.
"""
function DelaunayTriangulation(vertices::AbstractArray{Float64, 2}; joggle = 0.001)
    # Slightly joggle points to avoid problems with QHull
    if joggle > 0
        addnoise!(vertices; joggle_factor = joggle)
    end
    # Get the indices of the simplices of the triangulation
    if size(vertices, 1) > size(vertices, 2)
        DelaunayTriangulation(delaunay(vertices))
    elseif size(vertices, 1) < size(vertices, 2)
        DelaunayTriangulation(delaunay(transpose(vertices)))
    end
end

function DelaunayTriangulation(vertices::Vector{SVector{D, T}};
        joggle = 0.001) where {D, T}

    # Slightly joggle points to avoid problems with QHull
    if joggle > 0
        addnoise!(vertices; joggle_factor = joggle)
    end

    simplex_indices = delaunay(vertices)
    return DelaunayTriangulation(simplex_indices)
end

function DelaunayTriangulation(vertices::Vector{MVector{D, T}};
    joggle = 0.001) where {D, T}

# Slightly joggle points to avoid problems with QHull
if joggle > 0
    addnoise!(vertices; joggle_factor = joggle)
end

simplex_indices = delaunay(vertices)
return DelaunayTriangulation(simplex_indices)
end

"""
    DelaunayTriangulation(E::Embedding) -> DelaunayTriangulation

Construct a Delaunay triangulation from the points of an embedding.
"""
function DelaunayTriangulation(E::AbstractEmbedding; joggle = 0.001)
    DelaunayTriangulation(Array(transpose(E.points)), joggle = joggle)
end

"""
    DelaunayTriangulation(D::Dataset) -> DelaunayTriangulation

Construct a Delaunay triangulation from a Dataset.
"""
function DelaunayTriangulation(D::Dataset; joggle = 0.001)
    d = hcat(D.data...,)
    DelaunayTriangulation(d .+ joggle*rand(size(d)...,), joggle = joggle)
end

function DelaunayTriangulation(cr::CustomReconstruction; joggle = 0.001)
    DelaunayTriangulation(cr.reconstructed_pts; joggle = joggle)
end

export DelaunayTriangulation, dimension


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 34012b14-d597-4b9d-b23d-66b638e4c282
using JSON: escape_string

# ╔═╡ a8c1e0d2-3604-4e1d-a87c-c8f5b86b79ed
md"""
# MultiCheckBox
"""

# ╔═╡ c8350f43-0d30-45d0-873b-ff56c5801ac1
md"""
## Definition
"""

# ╔═╡ 144bff17-30eb-458a-8e94-33e1f11edbeb
"Convert a Julia array to a JS array in string form."
function jsarray_string(a::AbstractVector{T}) where {T <: AbstractString}
	string("[\"", join(map(escape_string, a), "\",\""), "\"]")
end

# ╔═╡ 91a08b98-52b5-4a2a-8180-7cba6d7232cd
function jsarray_string(a::AbstractVector{T}) where {T}
	string("[", join(a, ","), "]")
end

# ╔═╡ 430e2c1a-832f-11eb-024a-13e3989fd7c2
begin
	export MultiCheckBox

"""A group of checkboxes (`<input type="checkbox">`) - the user can choose enable or disable of the `options`, an array of `Strings`.
The value returned via `@bind` is a list containing the currently checked items.

See also: [`MultiSelect`](@ref).

`options` can also be an array of pairs `key::String => value::Any`. The `key` is returned via `@bind`; the `value` is shown.

`defaults` specifies which options should be checked initally.

`orientation` specifies whether the options should be arranged in `:row`'s `:column`'s.

`select_all` specifies whether or not to include a "Select All" checkbox.

# Examples
`@bind snacks MultiCheckBox(["🥕", "🐟", "🍌"]))`

`@bind snacks MultiCheckBox(["🥕" => "🐰", "🐟" => "🐱", "🍌" => "🐵"]; default=["🥕", "🍌"])`

`@bind animals MultiCheckBox(["🐰", "🐱" , "🐵", "🐘", "🦝", "🐿" , "🐝",  "🐪"]; orientation=:column, select_all=true)`"""
struct MultiCheckBox
    options::Array{Pair{<:AbstractString,<:Any},1}
    default::Union{Missing,AbstractVector{AbstractString}}
    orientation::Symbol
    select_all::Bool
end

MultiCheckBox(options::Array{<:AbstractString,1}; default=String[], orientation=:row, select_all=false) = MultiCheckBox([o => o for o in options], default, orientation, select_all)
MultiCheckBox(options::Array{<:Pair{<:AbstractString,<:Any},1}; default=String[], orientation=:row, select_all=false) = MultiCheckBox(options, default, orientation, select_all)

function Base.show(io::IO, ::MIME"text/html", multicheckbox::MultiCheckBox)
    if multicheckbox.orientation == :column
        flex_direction = "column"
    elseif multicheckbox.orientation == :row
        flex_direction = "row"
    else
        error("Invalid orientation $orientation. Orientation should be :row or :column")
    end

    js = read(joinpath(@__DIR__, "..", "assets", "multicheckbox.js"), String)
    css = read(joinpath(@__DIR__, "..", "assets", "multicheckbox.css"), String)

    labels = String[]
    vals = String[]
    default_checked = Bool[]
    for (k, v) in multicheckbox.options
        push!(labels, v)
        push!(vals, k)
        push!(default_checked, k in multicheckbox.default)
    end

    print(io, """
    <multi-checkbox class="multicheckbox-container" style="flex-direction:$(flex_direction);"></multi-checkbox>
    <script type="text/javascript">
        const labels = $(jsarray_string(labels));
        const values = $(jsarray_string(vals));
        const checked = $(jsarray_string(default_checked));
        const defaults = $(jsarray_string(multicheckbox.default));
        const includeSelectAll = $(multicheckbox.select_all);
        $(js)
    </script>
    <style type="text/css">
        $(css)
    </style>
    """)
end

Base.get(multicheckbox::MultiCheckBox) = multicheckbox.default
end

# ╔═╡ 8bfaf4c8-557d-433e-a228-aac493746efc
@bind animals MultiCheckBox(["🐰", "🐱" , "🐵", "🐘", "🦝", "🐿" , "🐝",  "🐪"]; orientation=:column, select_all=true)

# ╔═╡ 8e9f3962-d86c-4e07-b5d3-f31ee5361ca2
animals

# ╔═╡ 60183ad1-4919-4402-83fb-d53b86dda0a6
MultiCheckBox(["🐰 &&\\a \$\$", "🐱" , "🐵", "🐘", "🦝", "🐿" , "🐝",  "🐪"])

# ╔═╡ Cell order:
# ╟─a8c1e0d2-3604-4e1d-a87c-c8f5b86b79ed
# ╠═8bfaf4c8-557d-433e-a228-aac493746efc
# ╠═8e9f3962-d86c-4e07-b5d3-f31ee5361ca2
# ╠═60183ad1-4919-4402-83fb-d53b86dda0a6
# ╟─c8350f43-0d30-45d0-873b-ff56c5801ac1
# ╠═430e2c1a-832f-11eb-024a-13e3989fd7c2
# ╠═34012b14-d597-4b9d-b23d-66b638e4c282
# ╠═144bff17-30eb-458a-8e94-33e1f11edbeb
# ╠═91a08b98-52b5-4a2a-8180-7cba6d7232cd


text:
# Case I example: Bad diamond II
# Claudia January 2015

tree = "((((8,10))#H1,7),6,(4,#H1));" # Case I Bad diamond II
#f = open("prueba_tree.txt","w")
#write(f,tree)
#close(f)
net = readTopologyLevel1(tree)

#printEdges(net)
#printNodes(net)


text:
# Read the Bruker XML-based EDS spectrum file format.

"""
    readbrukerspx(fn::AbstractString)::Spectrum
    readbrukerspx(io::IO)::Spectrum

Read a Bruker SPX EDS spectrum file into a `Spectrum`.
"""
function readbrukerspx(filename::String)::Spectrum
    return open(filename, read = true) do io
        readbrukerspx(io)
    end
end

function readbrukerspx(io::IO)::Spectrum
    xml = readxml(io)
    # Read the counts data
    chdata = findall("//TRTSpectrum/ClassInstance/Channels", xml)[1].content
    counts = parse.(Int, strip.(split(chdata, ",")))
    # Read the header data
    props = Dict{Symbol,Any}()
    nrgy = missing
    item = findfirst(
        "//TRTSpectrum/ClassInstance/ClassInstance[@Type='TRTSpectrumHeader']",
        xml,
    )
    if !isnothing(item)
        dmy = Date(DateTime(findfirst("Date", item).content, "d.m.Y")) # <Date>10.5.2017</Date>
        hms = Time(DateTime(findfirst("Time", item).content, "H:M:S")) # <Time>11:36:22</Time>
        props[:AcquisitionTime] = DateTime(dmy, hms)
        @assert parse(Int, findfirst("ChannelCount", item).content) == length(counts)
        off = 1000.0 * parse(Float64, findfirst("CalibAbs", item).content) # <CalibAbs>-9.5639563E-1</CalibAbs>
        gain = 1000.0 * parse(Float64, findfirst("CalibLin", item).content) # <CalibLin>1.0001E-2</CalibLin>
        nrgy = LinearEnergyScale(off, gain)
    end
    item = findfirst(
        "//TRTSpectrum/ClassInstance/ClassInstance[@Type='TRTPSEElementList']",
        xml,
    )
    if !isnothing(item)
        data = findall("ChildClassInstances/ClassInstance", item)
        elmv = Element[
            PeriodicTable.elements[parse(Int, findfirst("Element", datum).content)] for
            datum in data
        ]
        if length(elmv) > 0
            props[:Elements] = elmv
        end
    end
    # item = findall("//TRTSpectrum/ClassInstance/ClassInstance[@Type='TRTResult']",xml) # Someday...
    item = findfirst(
        "//TRTSpectrum/ClassInstance/TRTHeaderedClass/ClassInstance[@Type='TRTSpectrumHardwareHeader']",
        xml,
    )
    item = findfirst(
        "//TRTSpectrum/ClassInstance/TRTHeaderedClass/ClassInstance[@Type='TRTSpectrumHardwareHeader']",
        xml,
    )
    if !isnothing(item)
        props[:RealTime] = 0.001 * parse(Int, findfirst("RealTime", item).content) # <RealTime>310610</RealTime>
        props[:LiveTime] = 0.001 * parse(Int, findfirst("LifeTime", item).content) # <LifeTime>300294</LifeTime>
        props[:BrukerThroughput] = parse(Int, findfirst("ShapingTime", item).content) # <ShapingTime>130000</ShapingTime>
    end
    item = findfirst(
        "//TRTSpectrum/ClassInstance/TRTHeaderedClass/ClassInstance[@Type='TRTDetectorHeader']",
        xml,
    )
    if !isnothing(item)
        props[:DetectorSerialNumber] = findfirst("Serial", item).content
        props[:DetectorModel] = "Bruker " * findfirst("Type", item).content
        thickness = 0.1 * parse(Float64, findfirst("DetectorThickness", item).content) # <DetectorThickness>0.45</DetectorThickness>
        props[:DetectorThickness] = thickness
        sideadlayer = parse(Float64, findfirst("SiDeadLayerThickness", item).content) # <SiDeadLayerThickness>0.029</SiDeadLayerThickness>
        props[:DeadLayerThickness] = 0.1 * sideadlayer
        wls = findfirst("WindowLayers", item)
        window = Film[]
        for i = 0:100
            wl = findfirst("Layer$i", wls)
            if isnothing(wl)
                break
            end
            z = parse(Int, wl["Atom"])  #   <Layer0 Atom="4" Thickness="1.3E1"/>
            thk = 1.0e-7 * parse(Float64, wl["Thickness"]) # nm
            relarea = haskey(wl, "RelativeArea") ? parse(Float64, wl["RelativeArea"]) : 1.0
            if relarea == 1.0
                push!(window, Film(pure(PeriodicTable.elements[z]), thk))
            end
        end
        if length(window) > 0
            props[:Window] = window
        end
    end
    item = findfirst(
        "//TRTSpectrum/ClassInstance/TRTHeaderedClass/ClassInstance[@Type='TRTXrfHeader']",
        xml,
    )
    if !isnothing(item)
        # <Voltage>50</Voltage>
        props[:BeamEnergy] = 1000.0 * parse(Float64, findfirst("Voltage", item).content)
        # <Anode>42</Anode>
        props[:XRFTubeAnode] =
            PeriodicTable.elements[parse(Int, findfirst("Anode", item).content)]
        # <Current>99</Current> ???microamps???
        props[:ProbeCurrent] = 1.0e3 * parse(Float64, findfirst("Current", item).content)
        # <TubeIncidentAngle>8.4E1</TubeIncidentAngle>
        props[:XRFTubeIncidentAngle] =
            deg2rad(parse(Float64, findfirst("TubeIncidentAngle", item).content))
        # <TubeTakeOffAngle>6</TubeTakeOffAngle>
        props[:XRFTubeTakeOffAngle] =
            deg2rad(parse(Float64, findfirst("TubeTakeOffAngle", item).content))
        # <ExcitationAngle>5E1</ExcitationAngle>
        props[:XRFExcitationAngle] =
            deg2rad(parse(Float64, findfirst("ExcitationAngle", item).content))
        # <DetectionAngle>5E1</DetectionAngle>
        props[:Elevation] =
            deg2rad(parse(Float64, findfirst("DetectionAngle", item).content))
        # <ExcitationPathLength>1E1</ExcitationPathLength>
        props[:XRFExcitationPathLength] =
            0.1 * parse(Float64, findfirst("ExcitationPathLength", item).content)
        # <DetectionPathLength>2E1</DetectionPathLength>
        props[:XRFDetectionPathLength] =
            0.1 * parse(Float64, findfirst("DetectionPathLength", item).content)
        props[:DetectorSolidAngle] =
            parse(Float64, findfirst("SolidAngleDetection", item).content)
        # <AzimutAngleAbs>0</AzimutAngleAbs>
        # <DetAzimutAngle>0</DetAzimutAngle>
        # <ChamberPressure>2E1</ChamberPressure>
        props[:ChamberPressure] = parse(Float64, findfirst("ChamberPressure", item).content)
        # <Atmosphere>Air</Atmosphere>
        props[:ChamberAtmosphere] = findfirst("Atmosphere", item).content
        # <TiltAngle>0</TiltAngle>
        props[:XRFSampleTilt] =
            deg2rad(parse(Float64, findfirst("TiltAngle", item).content))
        props[:TakeOffAngle] = props[:Elevation] + props[:XRFSampleTilt]
        welm = PeriodicTable.elements[parse(
            Int,
            findfirst("TubeWindow/AtomicNumber", item).content,
        )]
        wthk = 1.0e-4 * parse(Float64, findfirst("TubeWindow/Thickness", item).content)
        props[:XRFTubeWindow] = Film(pure(welm), wthk)
    end
    item = findfirst(
        "//TRTSpectrum/ClassInstance/TRTHeaderedClass/ClassInstance[@Type='TRTESMAHeader']",
        xml,
    )
    if !isnothing(item)
        # <PrimaryEnergy>2E1</PrimaryEnergy>
        props[:BeamEnergy] =
            1000.0 * parse(Float64, findfirst("PrimaryEnergy", item).content)
        props[:Elevation] = (
            props[:TakeOffAngle] =
                deg2rad(parse(Float64, findfirst("ElevationAngle", item).content))
        )
        props[:WorkingDistance] =
            0.1 * parse(Float64, findfirst("WorkingDistance", item).content)
    end
    #item = findfirst("//TRTSpectrum/ClassInstance/TRTHeaderedClass/ClassInstance[@Type='TRTXrfFPModelHeader']",xml)
    #item = findfirst("//TRTSpectrum/ClassInstance/TRTHeaderedClass/ClassInstance[@Type='TRTQuantitatorConfig']",xml)
    return Spectrum(nrgy, counts, props)
end

function detectbrukerspx(io::IO)
    try
        startswith(strip(readline(io)), "<?xml version=\"1.0\" encoding=") &&
            strip(readline(io)) == "<TRTSpectrum>"
    catch
        return false
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
include("../src/Utils/unionfind.jl")

using Test
import Base.==

==(r::Rep, i::Int) = r.i == i

@testset "UnionFind, mini_test_2         " begin
    u = UnionFind(2)
    @test find(u, 1) == 1
    @test find(u, 2) == 2

    unite!(u, find(u, 1), find(u, 1))
    @test find(u, 1) == 1
    @test find(u, 2) == 2

    unite!(u, find(u, 1), find(u, 2))
    @test find(u, 1) == find(u, 2)
end

@testset "UnionFind, mini_test_3         " begin
    u = UnionFind(3)
    @test find(u, 1) == 1
    @test find(u, 2) == 2
    @test find(u, 3) == 3

    unite!(u, find(u, 3), find(u, 1))

    @test find(u, 3) == find(u, 1)
    @test find(u, 2) != find(u, 1)
    @test find(u, 2) != find(u, 3)

    unite!(u, find(u, 2), find(u, 1))
    @test length(Set(find(u, i) for i = 1:3)) == 1
end

@testset "UnionFind, mini_test_4         " begin
    u = UnionFind(4)
    for i = 1:4
        @test find(u, i) == i
    end

    unite!(u, find(u, 4), find(u, 1))

    @test find(u, 4) == find(u, 1)
    @test find(u, 3) != find(u, 4)
    @test find(u, 3) != find(u, 1)

    unite!(u, find(u, 4), find(u, 3))

    @test find(u, 4) == find(u, 1)
    @test find(u, 3) == find(u, 4)
    @test find(u, 3) == find(u, 1)
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
using Distributed
using OpenTelemetryInstrumentationDistributed

addprocs(1)

@spawnat 1 identity("Hello World!")
@spawnat 2 identity("Hello World!")


text:
abstract type AbstractRisk{T <: DiseaseStateSequence} end

function Base.copy(x::AbstractRisk{SEIR})
  return AbstractRisk{SEIR}(copy(x.sparks),
                            copy(x.susceptibility),
                            copy(x.infectivity),
                            copy(x.transmissibility),
                            copy(x.latency),
                            copy(x.removal))
end

function Base.copy(x::AbstractRisk{SEI})
  return AbstractRisk{SEI}(copy(x.sparks),
                           copy(x.susceptibility),
                           copy(x.infectivity),
                           copy(x.transmissibility),
                           copy(x.latency))
end

function Base.copy(x::AbstractRisk{SIR})
  return AbstractRisk{SIR}(copy(x.sparks),
                           copy(x.susceptibility),
                           copy(x.infectivity),
                           copy(x.transmissibility),
                           copy(x.removal))
end

function Base.copy(x::AbstractRisk{SI})
  return AbstractRisk{SI}(copy(x.sparks),
                          copy(x.susceptibility),
                          copy(x.infectivity),
                          copy(x.transmissibility))
end

function _indices(x::AbstractRisk{T}; zeros::Bool=true, cumulative::Bool=true) where T <: DiseaseStateSequence
  indices = [length(x.sparks)
             length(x.susceptibility)
             length(x.infectivity)
             length(x.transmissibility)]
  if State_E ∈ T
    push!(indices, length(x.latency))
  elseif zeros
    push!(indices, 0)
  end
  if State_R ∈ T
    push!(indices, length(x.removal))
  elseif zeros
    push!(indices, 0)
  end
  if cumulative
    return cumsum(indices)
  else
    return indices
  end
end

function Base.length(x::AbstractRisk{T}) where T <: DiseaseStateSequence
  return _indices(x)[end]
end

function Base.getindex(x::AbstractRisk{T},
                       i::Int64) where T <: DiseaseStateSequence
  indices = _indices(x, zeros = true)
  riskfunc = findfirst(i .<= indices)
  return getfield(x, riskfunc)[end - (indices[riskfunc] - i)]
end


text:
module Bifurcation

using PyPlot

include("model/model.jl")
using .Model
if !isfile("./differential_equation.jl")
    run(`sh create_diffeq.sh`)
end
include("./differential_equation.jl")
include("../../continuation.jl")

function analysis()
    global fp, br

    # BP: name(index) of bifurcation parameter
    global BP = C.InhCDK

    fp = []
    ev = []
    br = []

    for i in 1:4
        global p = param_values()
        # i==1 -> Control
        if i == 2 # Wee1 inhibition
            p[C.kweeS] = 0.0
            p[C.kweeF] = 0.0
        elseif i == 3 # Gwl siRNA
            p[C.Gwtot] = 0.0
        elseif i == 4 # Wee1 inhibition & Gwl siRNA
            p[C.kweeS] = 0.0
            p[C.kweeF] = 0.0
            p[C.Gwtot] = 0.0
        end

        new_curve!(p, direction=false)
        push!(fp,readdlm("./Data/fp.dat",'\t',Float64,'\n'))
        push!(ev,readdlm("./Data/ev.dat",'\t',Float64,'\n'))
        push!(br,get_bistable_regime(ev[i]))
    end
end

function diagram()

    rc("figure",figsize = (20,3))
    rc("font",family = "Arial")
    rc("font",size = 14)
    rc("axes",linewidth = 1.2)
    rc("xtick.major",width = 1.2)
    rc("ytick.major",width = 1.2)
    rc("lines",linewidth = 2)

    for (i,(fixed_point,unstable_ss)) in enumerate(zip(fp,br))
        if length(unstable_ss) > 0
            intermediate_ss = []
            for j=2:length(unstable_ss)
                if unstable_ss[j] - unstable_ss[j-1] != 1
                    intermediate_ss = unstable_ss[j-1]:unstable_ss[j]
                end
            end
        end

        subplot(1,4,i)
        if i == 1
            plot(
                fixed_point[1:unstable_ss[1]-1,VN+1],
                fixed_point[1:unstable_ss[1]-1,V.Subp+1],
                color="royalblue"
            )
            plot(
                fixed_point[unstable_ss[1]:intermediate_ss[1]-1,VN+1],
                fixed_point[unstable_ss[1]:intermediate_ss[1]-1,V.Subp+1],
                color="darkgray","--"
            )
            plot(
                fixed_point[intermediate_ss,VN+1],
                fixed_point[intermediate_ss,V.Subp+1],
                color="darkorange"
            )
            plot(
                fixed_point[intermediate_ss[end]+1:unstable_ss[end],VN+1],
                fixed_point[intermediate_ss[end]+1:unstable_ss[end],V.Subp+1],
                color="darkgray","--"
            )
            plot(
                fixed_point[unstable_ss[end]+1:end,VN+1],
                fixed_point[unstable_ss[end]+1:end,V.Subp+1],
                color="crimson"
            )
            xlim(0,0.81)
            xlabel("1NMPP1 (μM)")
            ylim(-0.05,1.05)
            yticks([0,0.5,1],[0,50,100])
            ylabel("Sub-p (%)")
            title("Control",fontsize=18)
        elseif i == 2
            plot(
                fixed_point[1:unstable_ss[1]-1,VN+1],
                fixed_point[1:unstable_ss[1]-1,V.Subp+1],
                color="royalblue"
            )
            plot(
                fixed_point[unstable_ss,VN+1],
                fixed_point[unstable_ss,V.Subp+1],
                color="darkgray","--"
            )
            plot(
                fixed_point[unstable_ss[end]+1:end,VN+1],
                fixed_point[unstable_ss[end]+1:end,V.Subp+1],
                color="crimson"
            )
            xlim(0,0.81)
            xlabel("1NMPP1 (μM)")
            ylim(-0.05,1.05)
            yticks([0,0.5,1],[0,50,100])
            title("Wee1 inhibition",fontsize=18)
        elseif i == 3
            plot(
                fixed_point[1:unstable_ss[1]-1,VN+1],
                fixed_point[1:unstable_ss[1]-1,V.Subp+1],
                color="royalblue"
            )
            plot(
                fixed_point[unstable_ss,VN+1],
                fixed_point[unstable_ss,V.Subp+1],
                color="darkgray","--"
            )
            plot(
                fixed_point[unstable_ss[end]+1:end,VN+1],
                fixed_point[unstable_ss[end]+1:end,V.Subp+1],
                color="crimson"
            )
            xlim(0,0.81)
            xlabel("1NMPP1 (μM)")
            ylim(-0.05,1.05)
            yticks([0,0.5,1],[0,50,100])
            title("Gwl siRNA",fontsize=18)
        elseif i == 4
            plot(
                fixed_point[:,VN+1],
                fixed_point[:,V.Subp+1],
                color="darkorange"
            )
            xlim(0,0.81)
            xlabel("1NMPP1 (μM)")
            ylim(-0.05,1.05)
            yticks([0,0.5,1],[0,50,100])
            title("Wee1 inhibition & Gwl siRNA",fontsize=18)
        end
    end
    show()
end
end  # module

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function posang(units::Integer, ra1::T, dec1::T, ra2::T, dec2::T) where {T<:AbstractFloat}
    # Convert all quantities to radians.
    if units == 0
        # All radians
        ra1_rad  = ra1
        ra2_rad  = ra2
        dec1_rad = dec1
        dec2_rad = dec2
    elseif units == 1
        # Right ascensions are in hours, declinations in degrees.
        ra1_rad  = ra1*pi/12
        ra2_rad  = ra2*pi/12
        dec1_rad = deg2rad(dec1)
        dec2_rad = deg2rad(dec2)
    elseif units == 2
        # Right ascensions and declinations are in degrees.
        ra1_rad  = deg2rad(ra1)
        ra2_rad  = deg2rad(ra2)
        dec1_rad = deg2rad(dec1)
        dec2_rad = deg2rad(dec2)
    else
        # In any other case throw an error.
        error("units must be 0 (radians), 1 (hours, degrees) or 2 (degrees)")
    end
    sin_radif, cos_radif = sincos(ra2_rad - ra1_rad)
    sin_dec1, cos_dec1 = sincos(dec1_rad)
    angle = atan(sin_radif, cos_dec1 * tan(dec2_rad) -
                  sin_dec1 * cos_radif)
    if units == 0
        return angle
    else
        return rad2deg(angle)
    end
end

"""
    posang(units, ra1, dec1, ra2, dec2) -> angular_distance

### Purpose ###

Compute rigorous position angle of point 2 relative to point 1.

### Explanation ###

Computes the rigorous position angle of point 2 (with given right ascension and
declination) using point 1 (with given right ascension and declination) as the
center.

### Arguments ###

* `units`: integer, can be either 0, or 1, or 2.  Describes units of inputs and
 output:
    * 0: everything (input right ascensions and declinations, and output
      distance) is radians
    * 1: right ascensions are in decimal hours, declinations in decimal degrees,
      output distance in degrees
    * 2: right ascensions and declinations are in degrees, output distance in
      degrees
* `ra1`:  right ascension or longitude of point 1
* `dec1`: declination or latitude of point 1
* `ra2`: right ascension or longitude of point 2
* `dec2`: declination or latitude of point 2

Both `ra1` and `dec1`, and `ra2` and `dec2` can be given as 2-tuples `(ra1,
dec1)` and `(ra2, dec2)`.

### Output ###

Angle of the great circle containing `[ra2, dec2]` from the meridian containing
`[ra1, dec1]`, in the sense north through east rotating about `[ra1, dec1]`.
See `units` argument above for units.

### Method ###

The "four-parts formula" from spherical trigonometry (p. 12 of Smart's Spherical
Astronomy or p. 12 of Green' Spherical Astronomy).

### Example ###

Mizar has coordinates (ra, dec) = (13h 23m 55.5s, +54° 55' 31'').  Its
companion, Alcor, has coordinates (ra, dec) = (13h 25m 13.5s, +54° 59' 17'').
Find the position angle of Alcor with respect to Mizar.

```jldoctest
julia> using AstroLib

julia> posang(1, ten(13, 25, 13.5), ten(54, 59, 17), ten(13, 23, 55.5), ten(54, 55, 31))
-108.46011246802047
```

### Notes ###

* The function `sphdist` provides an alternate method of computing a spherical
 distance.
* Note that `posang` is not commutative: the position angle between A and B is
  \$\\theta\$, then the position angle between B and A is \$180 + \\theta\$.

Code of this function is based on IDL Astronomy User's Library.
"""
posang(units::Integer, ra1::Real, dec1::Real, ra2::Real, dec2::Real) =
    posang(units, promote(float(ra1), float(dec1), float(ra2), float(dec2))...)

### Tuples input
posang(units::Integer, radec1::Tuple{Real, Real}, ra2::Real, dec2::Real) =
    posang(units, radec1..., ra2, dec2)

posang(units::Integer, ra1::Real, dec1::Real, radec2::Tuple{Real, Real}) =
    posang(units, ra1, dec1, radec2...)

posang(units::Integer, radec1::Tuple{Real, Real}, radec2::Tuple{Real, Real}) =
    posang(units, radec1..., radec2...)


text:
using Resparse
using Test
using  ForwardDiff,GenericLinearAlgebra,LinearAlgebra

@testset "Resparse.jl" begin
    # Write your tests here.

    A = randn(2, 2)
    f(x) = svd(x).S |> sum
    g(x) = Resparse.svd_lowrank(x, 2, 2).S |> sum
    @time ForwardDiff.gradient(g, A)

end

max_stars_repo_name: SamanthaWangdl/ReSparse

max_issues_repo_name: SamanthaWangdl/ReSparse

max_forks_repo_name: SamanthaWangdl/ReSparse


text:
export std_dev, variance, dist, rms, angle, gram_schmidt

function std_dev(x)
  return vecnorm(x - mean(x))/sqrt(length(x))
end

function variance(x)
  return sum((x - mean(x)).^2)/length(x)
end

function dist(x, y)
  return vecnorm(x - y)
end

function rms(x)
  return vecnorm(x)/sqrt(length(x))
end

function angle(x, y)
  return acos(dot(x, y)/(vecnorm(x) * vecnorm(y)))
end

# Runs the Gram-Schmidt algorithm for the given set
# of vectors, returning either an orthonormal basis or
# the index of the first vector which can be written
# as a linear combination of the previous vectors
#
# parameters:
#   a   - An array of column vectors
#
# returns:
#   If 'a' is a set of n linearly independent vectors, then
#   gram_schmidt returns an array n of orthonormal vectors.
#
#   If 'a' is not a set of linearly independent vectors, then
#   gram_schmidt returns an array of m < n vectors.
#
# usage:
#   julia> gram_schmidt({[1;0],[0;1][1;1]})
#   1-element Array{Any,1}:
#   [1.0,0.0]
#
#   julia> gram_schmidt({[1;0;0],[0;1;0],[0;0;1]})
#   3-element Array{Any,1}:
#    [1.0,0.0,0.0]
#    [0.0,1.0,0.0]
#    [0.0,0.0,1.0]
#
function gram_schmidt(A)
  epsilon = 1e-6
  k = length(A)
  n = length(A[1])
  qs = []
  for i = 1:k
    q_i = A[i]
    # orthogonalization with previous vectors
    for j = 1:length(qs)
      q_i -= (qs[j]'A[i]).*qs[j]
    end
    # check for dependence, quit if q_i ~= 0
    if norm(q_i) <= epsilon
      break
    end
    # normalize and store the result
    q_i = q_i/norm(q_i)
    push!(qs,q_i)
  end
  return qs
end