text:
function get_variable_cost(
    ts::IS.TimeSeriesData,
    component::Component,
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    if start_time === nothing
        start_time = IS.get_initial_timestamp(ts)
    end
    data = IS.get_time_series_array(component, ts, start_time, len = len)
    time_stamps = TimeSeries.timestamp(data)
    return TimeSeries.TimeArray(time_stamps, map(VariableCost, TimeSeries.values(data)))
end

function get_variable_cost(
    device::StaticInjection,
    cost::OperationalCost;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    time_series_key = get_variable(cost)
    if isnothing(time_series_key)
        error("Cost component has a `nothing` stored in field `variable`, Please use `set_variable_cost!` to add variable cost forecast.")
    end
    raw_data = IS.get_time_series_by_key(
        time_series_key,
        device;
        start_time = start_time,
        len = len,
        count = 1,
    )
    cost = get_variable_cost(raw_data, device, start_time, len)
    return cost
end

function get_variable_cost(
    service::ReserveDemandCurve;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    time_series_key = get_variable(service)
    if isnothing(time_series_key)
        error("Cost component has a `nothing` stored in field `variable`, Please use `set_variable_cost!` to add variable cost forecast.")
    end
    raw_data = IS.get_time_series_by_key(
        time_series_key,
        service;
        start_time = start_time,
        len = len,
        count = 1,
    )
    cost = get_variable_cost(raw_data, service, start_time, len)
    return cost
end

function get_services_bid(
    device::StaticInjection,
    cost::MarketBidCost,
    service::Service;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    variable_ts_key = get_variable(cost)
    service_ts_key = IS.TimeSeriesKey(variable_ts_key.time_series_type, get_name(service))
    raw_data = IS.get_time_series_by_key(
        service_ts_key,
        device;
        start_time = start_time,
        len = len,
        count = 1,
    )
    cost = get_variable_cost(raw_data, device, start_time, len)
    return cost
end

function set_variable_cost!(
    sys::System,
    component::StaticInjection,
    time_series_data::IS.TimeSeriesData,
)
    add_time_series!(sys, component, time_series_data)
    key = IS.TimeSeriesKey(time_series_data)
    market_bid_cost = get_operation_cost(component)
    set_variable!(market_bid_cost, key)
    return
end

function set_variable_cost!(
    sys::System,
    component::ReserveDemandCurve,
    time_series_data::IS.TimeSeriesData,
)
    add_time_series!(sys, component, time_series_data)
    key = IS.TimeSeriesKey(time_series_data)
    set_variable!(component, key)
    return
end

function set_service_bid!(
    sys::System,
    component::StaticInjection,
    service::Service,
    time_series_data::IS.TimeSeriesData,
)
    if get_name(time_series_data) != get_name(service)
        error("Name provided in the TimeSeries Data $(get_name(time_series_data)), doesn't match the Service $(get_name(service)).")
    end
    verify_device_eligibility(sys, component, service)
    add_time_series!(sys, component, time_series_data)
    ancillary_services = get_ancillary_services(get_operation_cost(component))
    push!(ancillary_services, service)
    return
end

function verify_device_eligibility(
    sys::System,
    component::StaticInjection,
    service::Service,
)
    contributing_devices = get_contributing_devices(sys, service)
    if !in(component, contributing_devices)
        error("Device $(get_name(component)) isn't eligible to contribute to service $(get_name(service)).")
    end
    return
end


text:
# Tue 13 Nov 2012 04:13:36 PM EST
#
# NeHe Tut 8 - Implement lights and rotate an alpha-blended, textured cube
#
# Q - quit
# B - turn texture alpha-blending on/off
# L - turn lights on/off
# F - change texture filter (linear, nearest, mipmap)
# PageUp/Down - move camera closer/further away from cube
# Up/Down - increase/decrease x-rotation speed
# Left/Right - increase/decrease y-rotation speed

# load necessary GL/SDL routines

using OpenGL
@OpenGL.version "1.0"
@OpenGL.load
using SDL

# initialize variables

bpp                = 16
wintitle           = "NeHe Tut 8"
icontitle          = "NeHe Tut 8"
width              = 640
height             = 480

light              = true
blend              = true
filter             = 3

xrot               = 0.0
yrot               = 0.0
xspeed             = 0.0
yspeed             = 0.0

z                  = -5.0

cube_size          = 1.0

LightAmbient       = [0.5f0, 0.5f0, 0.5f0, 1.0f0]
LightDiffuse       = [1.0f0, 1.0f0, 1.0f0, 1.0f0]
LightPosition      = [0.0f0, 0.0f0, 2.0f0, 1.0f0]

keystate_checked   = false
lastkeycheckTime   = 0
key_repeatinterval = 75 #ms

# open SDL window with an OpenGL context

SDL_Init(SDL_INIT_VIDEO)
#videoInfo = SDL_GetVideoInfo()
videoFlags = (SDL_OPENGL | SDL_GL_DOUBLEBUFFER | SDL_HWPALETTE | SDL_RESIZABLE)
#if videoInfo.hw_available
    videoFlags = (videoFlags | SDL_HWSURFACE)
#else
    #videoFlags = (videoFlags | SDL_SWSURFACE)
#end
#if videoInfo.blit_hw
    videoFlags = (videoFlags | SDL_HWACCEL)
#end
SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1)
SDL_SetVideoMode(width, height, bpp, videoFlags)
SDL_WM_SetCaption(wintitle, icontitle)

glViewport(0, 0, width, height)
glClearColor(0.0, 0.0, 0.0, 0.0)
glClearDepth(1.0)
glDepthFunc(GL_LESS)
glEnable(GL_DEPTH_TEST)
glShadeModel(GL_SMOOTH)

glMatrixMode(GL_PROJECTION)
glLoadIdentity()

gluPerspective(45.0,width/height,0.1,100.0)

glMatrixMode(GL_MODELVIEW)

### auxiliary functions

function cube(size)  # the cube function now includes surface normal specification for proper lighting
  glBegin(GL_QUADS)
    # Front Face
    glNormal(0.0,0.0,1.0)
    glTexCoord(0.0, 0.0)
    glVertex(-size, -size, size)
    glTexCoord(1.0, 0.0)
    glVertex(size, -size, size)
    glTexCoord(1.0, 1.0)
    glVertex(size, size, size)
    glTexCoord(0.0, 1.0)
    glVertex(-size, size, size)

    # Back Face
    glNormal(0.0,0.0,-1.0)
    glTexCoord(1.0, 0.0)
    glVertex(-size, -size, -size)
    glTexCoord(1.0, 1.0)
    glVertex(-size, size, -size)
    glTexCoord(0.0, 1.0)
    glVertex(size, size, -size)
    glTexCoord(0.0, 0.0)
    glVertex(size, -size, -size)

    # Top Face
    glNormal(0.0,1.0,0.0)
    glTexCoord(0.0, 1.0)
    glVertex(-size, size, -size)
    glTexCoord(0.0, 0.0)
    glVertex(-size, size, size)
    glTexCoord(1.0, 0.0)
    glVertex(size, size, size)
    glTexCoord(1.0, 1.0)
    glVertex(size, size, -size)

    # Bottom Face
    glNormal(0.0,-1.0,0.0)
    glTexCoord(1.0, 1.0)
    glVertex(-size, -size, -size)
    glTexCoord(0.0, 1.0)
    glVertex(size, -size, -size)
    glTexCoord(0.0, 0.0)
    glVertex(size, -size, size)
    glTexCoord(1.0, 0.0)
    glVertex(-size, -size, size)

    # Right Face
    glNormal(1.0,0.0,0.0)
    glTexCoord(1.0, 0.0)
    glVertex(size, -size, -size)
    glTexCoord(1.0, 1.0)
    glVertex(size, size, -size)
    glTexCoord(0.0, 1.0)
    glVertex(size, size, size)
    glTexCoord(0.0, 0.0)
    glVertex(size, -size, size)

    # Left Face
    glNormal(-1.0,0.0,0.0)
    glTexCoord(0.0, 0.0)
    glVertex(-size, -size, -size)
    glTexCoord(1.0, 0.0)
    glVertex(-size, -size, size)
    glTexCoord(1.0, 1.0)
    glVertex(-size, size, size)
    glTexCoord(0.0, 1.0)
    glVertex(-size, size, -size)
  glEnd()
end

### end of auxiliary functions

# load textures from images

tex   = Array(Uint32,3) # generating 3 textures

img, w, h = glimread(expanduser("~/.julia/SDL/Examples/NeHe/tut8/glass.bmp"))

glGenTextures(3,tex)
glBindTexture(GL_TEXTURE_2D,tex[1])
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)

glBindTexture(GL_TEXTURE_2D,tex[2])
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)

glBindTexture(GL_TEXTURE_2D,tex[3])
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST)
glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)

gluBuild2DMipmaps(GL_TEXTURE_2D, 3, w, h, GL_RGB, GL_UNSIGNED_BYTE, img)

# initialize lights

glLightfv(GL_LIGHT1, GL_AMBIENT, LightAmbient)
glLightfv(GL_LIGHT1, GL_DIFFUSE, LightDiffuse)
glLightfv(GL_LIGHT1, GL_POSITION, LightPosition)

glEnable(GL_LIGHT1)

# enable texture mapping

glEnable(GL_TEXTURE_2D)
glEnable(GL_LIGHTING)
glEnable(GL_BLEND)

# enable alpha blending for textures

glBlendFunc(GL_SRC_ALPHA, GL_ONE)
glColor(1.0, 1.0, 1.0, 0.5)

# main drawing loop

while true
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    glTranslate(0.0, 0.0, z)

    glRotate(xrot,1.0,0.0,0.0)
    glRotate(yrot,0.0,1.0,0.0)

    glBindTexture(GL_TEXTURE_2D,tex[filter])
    cube(cube_size)

    xrot +=xspeed
    yrot +=yspeed

    SDL_GL_SwapBuffers()

    SDL_PumpEvents()
    if SDL_GetTicks() - lastkeycheckTime >= key_repeatinterval
        keystate         = SDL_GetKeyState()
        keystate_checked = true
        lastkeycheckTime = SDL_GetTicks()
    end

    # Sampling rates for event processing are so fast that a single key press
    # lasts through several iterations of this loop.  This means that one press
    # can be seen as 50 or more presses by the SDL event system.  To correct
    # for this, we only check keypresses every 75ms.

    if keystate_checked == true
        if keystate[SDLK_q] == true
            break
        elseif keystate[SDLK_b] == true
            println("Blend was: $blend")
            blend = (blend ? false : true)
            if blend
                glEnable(GL_BLEND)
                glDisable(GL_DEPTH_TEST)
            else
                glDisable(GL_BLEND)
                glEnable(GL_DEPTH_TEST)
            end
            println("Blend is now: $blend")
        elseif keystate[SDLK_l] == true
            println("Light was: $light")
            light = (light ? false : true)
            println("Light is now: $light")
            if light
                glEnable(GL_LIGHTING)
            else
                glDisable(GL_LIGHTING)
            end
        elseif keystate[SDLK_f] == true
            println("Filter was: $filter")
            filter += 1
            if filter > 3
                filter = 1
            end
            println("Filter is now: $filter")
        elseif keystate[SDLK_PAGEUP] == true
            z -= 0.02
        elseif keystate[SDLK_PAGEDOWN] == true
            z += 0.02
        elseif keystate[SDLK_UP] == true
            xspeed -= 0.01
        elseif keystate[SDLK_DOWN] == true
            xspeed += 0.01
        elseif keystate[SDLK_LEFT] == true
            yspeed -= 0.01
        elseif keystate[SDLK_RIGHT] == true
            yspeed += 0.01
        end
        keystate_checked = false
    end
end

ext: jl


text:
# Dict (lexicalItem, count) to keep track of previous expected counts for vb_update
#type PrevCounts
#    prev_counts :: Dict{ Vector{T} , Dict{LexicalItem{T}, Float64}}
#end
#function PrevCounts()
#    prev_counts
#    PrevCounts(dict())
#end

#function getcatfeat(lex :: LexicalItem{T})
#    return lex.features[length(features)]
#end

function vb_update(counts, prev_counts, grammar)
    # for each lexical item with new expected count
    for (s, dict) in counts
        if !(s in keys(prev_counts))
            prev_counts[s] = Dict()
        end
        for (id, p) in dict
            #get previous count if subtracting it to new count
            prev=0.0
            if id in keys(prev_counts[s])
                prev = prev_counts[s][id]
            end
            #update previous counts with new counts for next iterations
            prev_counts[s][id] = p
            #get category of idical item (will need to change this line when movement is introduced)
            cat = catfeature(features(LexicalItem(id)))
            # update count in dirichlet categorical which will automatically update idical items score
            grammar.lexitem_cond.dists[cat].counts[id] += (p - prev)
        end
    end
    return (prev_counts, grammar)
end

function variational_inference(subset_size, corpus_size, iterations, prev_counts, grammar, corpus, epsilon=true)
    # iterate as many times as desired for updates
    for i in (1:iterations)

        # random sample of subset_size (eg. 20) sentences from corpus called subset
        subset_ids = rand(1:corpus_size, subset_size)
        subset = corpus[subset_ids]
        # Get expected counts for all lexical items in parses of subset
        counts = expected_lexitem_counts_by_sentence(subset, grammar, epsilon)
        # variational Bayesian update given counts
        prev_counts, grammar = vb_update(counts, prev_counts, grammar)

        # write a safety function that saves parameters every so many iterations in case computer crashes
        n = string(i)
        name = string(n,"_MG.txt")
        open(name, "w") do f
            for cat in keys(grammar.lexitem_cond.dists)
                for id in keys(grammar.lexitem_cond.dists[cat].counts)
                    item = grammar.lexitem_cond.dists[cat].counts[id]
                    lex = LexicalItem(id)
                    write(f, "$(lex.phon) ::")
                    map(x -> write(f, " $x"), lex.features)
                    prob = string(grammar.lexitem_cond.dists[cat].counts[id])
                    write(f, " ; $prob")
                    write(f, "\n")
                end
            end
        end
    end
    return grammar
end

##################################
###### FULL CORPUS 1 CORE   ######
##################################

function variational_inference_full_corpus(iterations, prev_counts, grammar, corpus, result_dir, epsilon=true)
    # iterate as many times as desired for updates
    for i in (1:iterations)

        # Get expected counts for all lexical items in parses of corpus
        counts = expected_lexitem_counts_by_sentence(corpus, grammar, epsilon)
        # variational Bayesian update given counts
        prev_counts, grammar = vb_update(counts, prev_counts, grammar)

        # write a safety function that saves parameters every so many iterations in case computer crashes
        n = string(i)
        name = string(result_dir,"/", n,"_MG.txt")
        open(name, "w") do f
            for cat in keys(grammar.lexitem_cond.dists)
                for id in keys(grammar.lexitem_cond.dists[cat].counts)
                    item = grammar.lexitem_cond.dists[cat].counts[id]
                    write(f, id)
                    prob = string(grammar.lexitem_cond.dists[cat].counts[id])
                    write(f, " ; $prob")
                    write(f, "\n")
                end
            end
        end
    end
    return grammar
end

##################################
###### PARALLEL COMPUTATION ######
##################################

"give ranges to partition an array of lengh nb_data into subarrays of length nb_data_per_chunk https://stackoverflow.com/a/43767567"
function partition_array_indices(nb_data::Int, nb_data_per_chunk::Int)
    nb_chunks = ceil(Int, nb_data / nb_data_per_chunk)
    ids = UnitRange{Int}[]
    for which_chunk = 1:nb_chunks
        id_start::Int = 1 + nb_data_per_chunk * (which_chunk - 1)
        id_end::Int = id_start - 1 + nb_data_per_chunk
        if id_end > nb_data
            id_end = nb_data
        end
        push!(ids, id_start:id_end)
    end
    return ids
end

## non-parallel version
#function parse_and_calculate(grammar, corpus :: Vector{Vector{String}}, epsilon :: Bool, exact_counts :: Bool, samples_multiplier=0.5 :: Float64)
#    counts_by_sent = Dict{Vector{String}, Dict{String, Float64}}() # Dict sent => Dict lexid => float exp. count
#    cl = LogProb(1.0) # corpus likelihood
#
#    function p_c(s::Vector{String})
#        f = epsilon ? run_chartparser(s, grammar, epsilon="") : run_chartparser(s, grammar)
#        c = exact_counts ? expected_lexitem_counts(f, grammar, epsilon) : sample_counts(f, grammar, epsilon, true, samples_multiplier) :: Dict{String, Float64}
#        l = length(heads(f)) == 0 ? LogProb(0.0) : score(f)
#        l, c
#    end
#
#    for (j, (l, c)) in enumerate(map(p_c, corpus))
#        cl = cl * l
#        counts_by_sent[corpus[j]] = c
#    end
#    (cl, counts_by_sent)
#end

function parse_and_calculate(grammar, corpus :: Vector{Vector{String}}, epsilon :: Bool, exact_counts :: Bool, cores :: Int, samples_multiplier=0.5 :: Float64)
    counts_by_sent = Dict{Vector{String}, Dict{String, Float64}}() # Dict sent => Dict lexid => float exp. count
    cl = LogProb(1.0) # corpus likelihood

    function p_c(s::Vector{String})
        f = epsilon ? run_chartparser(s, grammar, epsilon="") : run_chartparser(s, grammar)
        c = exact_counts ? expected_lexitem_counts(f, grammar, epsilon) : sample_counts(f, length(s), grammar, epsilon, true, samples_multiplier) :: Dict{String, Float64}
        l = length(heads(f)) == 0 ? LogProb(0.0) : score(f)
        l, c
    end

    println("parse and calculate for each sentence")
    @time likelihoods_counts = pmap(p_c, corpus)
    flush(STDOUT); println("assemble calculations")
    @time cl = prod(c_l[1] for c_l in likelihoods_counts)
    #@time counts_by_sent = merge((c_l[2] for c_l in likelihoods_counts)...)
    @time counts_by_sent = Dict(corpus[i] => likelihoods_counts[i][2] for i in 1:length(corpus))
    flush(STDOUT); println("finished parse and calculate step"); flush(STDOUT)
    (cl, counts_by_sent)
end

## parallel version, spawns non parallel version on each core (simultaneously)
#function parse_and_calculate(grammar, corpus :: Vector{Vector{String}}, epsilon :: Bool, exact_counts :: Bool, cores :: Int, samples_multiplier=0.5 :: Float64)
#    counts_by_sent = Dict() # Dict sent => Dict lexid => float exp. count
#    cl = LogProb(1.0) # corpus likelihood
#    batch_ranges = partition_array_indices(length(corpus), Int(ceil(length(corpus)/cores)))
#    batches = [corpus[r] for r in batch_ranges]
#    p_c(b) = parse_and_calculate(grammar, b, epsilon, exact_counts, samples_multiplier)
#    counts_likelihoods = pmap(p_c, batches)
#    cl = prod(c_l[1] for c_l in counts_likelihoods)
#    counts_by_sent = merge((c_l[2] for c_l in counts_likelihoods)...)
#    println("finished parse and calculate step")
#    (cl, counts_by_sent)
#end

function variational_inference_full_corpus_parallel(iterations,
                                                    prev_counts,
                                                    grammar,
                                                    corpus,
                                                    cores,
                                                    result_dir,
                                                    epsilon=true,
                                                    transitive=true,
                                                    directed=true,
                                                    saveprogress=true,
                                                    exact_counts=true,
                                                    samples_multiplier=0.5)
    # stats to save
    kls = []
    likelihoods = []

    # iterate as many times as desired for updates
    #corpus = [[word.form for word in csent.words] for csent in traincorpus.sents]
    #subset = Int(floor(length(corpus)/cores))
    #println("Parse with grammar 0")
    #@time forests = epsilon ? pmap(s -> run_chartparser(s, grammar, epsilon=""), corpus) : pmap(s -> run_chartparser(s, grammar), corpus)
    #println("Calculate corpus likelihood")
    #@time cl = corpus_likelihood(forests, epsilon, true)
    #push!(likelihoods, cl)
    #@time println("likelihood of the training corpus: $cl\n")
    for i in 1:iterations
        println("BEGIN ITERATION $i")
        println("get counts and corpus likelihood for grammar $(i-1)")
        #@time counts = exact_counts ? expected_lexitem_counts_by_sentence(forests, grammar, epsilon) : sample_counts_by_sentence(forests, grammar, epsilon)
        @time cl, counts = parse_and_calculate(grammar, corpus, epsilon, exact_counts, cores, samples_multiplier)
        println("update counts")
    	@time prev_counts, new_grammar = vb_update(counts, prev_counts, grammar)

        # write a safety function that saves parameters every so many iterations in case computer crashes
        if saveprogress
            n = string(i)
            name = string(result_dir,"/", n,"_MG.txt")
            open(name, "w") do f
                for cat in keys(grammar.lexitem_cond.dists)
                    for id in keys(grammar.lexitem_cond.dists[cat].counts)
                        item = grammar.lexitem_cond.dists[cat].counts[id]
                        write(f, id)
                        prob = string(grammar.lexitem_cond.dists[cat].counts[id])
                        write(f, " ; $prob")
                        write(f, "\n")
                    end
                end
            end
        end

        #println("Parse with grammar $i")
        #@time forests = epsilon ? pmap(s -> run_chartparser(s, new_grammar, epsilon=""), corpus) : pmap(s -> run_chartparser(s, new_grammar), corpus)

        # calculate accuracy at the current iteration
        println("END ITERATAION $i: ")

        println("Calculate KL")
        @time kl = kullback_leibler(grammar, new_grammar)
        push!(kls, kl)
        println("KL divergence: $kl")

        #println("Calculate corpus likelihood")
        #@time cl = corpus_likelihood(forests, epsilon, true)
        push!(likelihoods, cl)
        @time println("likelihood of the training corpus: $cl")

        println("")
        flush(STDOUT)
        grammar = new_grammar
    end

    # print summary of stats
    println("")
    println("END VB: SUMMARY")
    println("===============")
    print("KL:")
    for x in kls
        print("\t$x")
    end
    println("")
    print("likelihoods:")
    for x in likelihoods
        print("\t$x")
    end
    println("")
    print("likelihoods:")
    exp_likelihoods = map(exp, likelihoods)
    for x in exp_likelihoods
        print("\t$x")
    end
    println("")
    flush(STDOUT)

    return grammar
end

##################################
###### LOAD SAVED LEXICON   ######
##################################
function load_lexicon(filename)
    lexicon = ""
    f = open(filename)
        lines = readlines(f)
        for l in lines
            lexicon = string(lexicon, l, "\n")
        end

    close(f)
    return lexicon
end


text:
using Test
using PyBaMM
using SparseArrays, LinearAlgebra
using Sundials
using PyCall

pybamm = pyimport("pybamm")

@testset "Compare with PyBaMM: SPM" begin
    # load model
    model = pybamm.lithium_ion.SPM(name="SPM")
    sim = pybamm.Simulation(model)

    prob = get_ode_problem(sim)

    sol = solve(prob, CVODE_BDF(), saveat=prob.tspan[2] / 100);

    # Calculate voltage in Julia
    V = get_variable(sim, sol, "Terminal voltage [V]")
    t = get_variable(sim, sol, "Time [s]")

    # Solve in python
    sol_pybamm = sim.solve(t)
    V_pybamm = get(sol_pybamm, "Terminal voltage [V]").data
    @test all(isapprox.(V_pybamm, V, atol=1e-3))
end

@testset "Compare with PyBaMM: SPMe" begin
    # load model
    model = pybamm.lithium_ion.SPMe(name="SPMe")
    sim = pybamm.Simulation(model)

    prob = get_ode_problem(sim)

    sol = solve(prob, CVODE_BDF(), saveat=prob.tspan[2] / 100);

    # Calculate voltage in Julia
    V = get_variable(sim, sol, "Terminal voltage [V]")
    t = get_variable(sim, sol, "Time [s]")

    # Solve in python
    sol_pybamm = sim.solve(t)
    V_pybamm = get(sol_pybamm, "Terminal voltage [V]").data
    @test all(isapprox.(V_pybamm, V, atol=1e-4))
end

@testset "Compare with PyBaMM: DFN" begin
    # load model
    model = pybamm.lithium_ion.DFN(name="DFN")
    sim = pybamm.Simulation(model)

    prob = get_dae_problem(sim)

    sol = solve(prob, IDA(), saveat=prob.tspan[2] / 100);

    # Calculate voltage in Julia
    V = get_variable(sim, sol, "Terminal voltage [V]")
    t = get_variable(sim, sol, "Time [s]")

    # Solve in python
    sol_pybamm = sim.solve(t)
    V_pybamm = get(sol_pybamm, "Terminal voltage [V]").data
    @test all(isapprox.(V_pybamm, V, atol=1e-3))
end


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: December 2020
# Copyright: Georgia Institute of Technology, 2020

export _wsave

_wsave(s, fig::Figure; dpi::Int=300) = fig.savefig(s, bbox_inches="tight", dpi=dpi)


text:
using Compat, Calculus, QuadGK
using Test: @test

function psm(nip)
  @test nip >= 1
  # Integrate to compute coefficient in rhs
  rhs(n, nip) = QuadGK.quadgk((x)-> x^nip, -n, n)[1]
  local hn
  if nip == 1
    w = [1.0]
  elseif nip == 2
    w = [1/2, 1/2]
  else
    if isodd(nip)
      n = floor(Int, nip/2)
      hn = range(-n, stop=n, length=nip)
    else
      n = floor(Int, 2*(nip/2)-1)
      hn = range(-n, stop=n, length=nip)
    end
    m = ones(nip, nip)
    for i in 0:nip-1
      for j in 1:nip
        m[i+1, j] = (hn[j])^i
      end
    end
    v = [rhs(n, ni) for ni in 0:nip-1]
    w = m \ v
    iseven(nip) && (w = w / 2.0)
  end
  w'
end

function nc(func, lb, ub, nip)
  @test ub > lb
  @test nip >= 1
  if nip == 1
    h = ub - lb
    x = [lb]
  elseif nip ==2
    h = ub - lb
    x = [lb, ub]
  else
    h = (ub - lb)/(nip-1)
    x = collect(range(lb, stop=ub, length=nip))
  end
  w = psm(nip)
  res = h * w *func.(x)
  (w, res)
end

lb = pi/4
ub = pi/2
func = sin
nip = 7
w, res = nc(func, lb, ub, nip)

@test res[1] ≈ QuadGK.quadgk(sin, pi/4, pi/2)[1]


text:
# # Flock model
# ![](flock.gif)
#
# The flock model illustrates how flocking behavior can emerge when each bird follows three simple rules:
#
# * maintain a minimum distance from other birds to avoid collision
# * fly towards the average position of neighbors
# * fly in the average direction of neighbors

# ## Defining the core structures

# We begin by calling the required packages and defining an agent type representing a bird.

using Agents, Random, LinearAlgebra

mutable struct Bird <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    speed::Float64
    cohere_factor::Float64
    separation::Float64
    seperate_factor::Float64
    match_factor::Float64
    visual_distance::Float64
end

# The fields `id` and `pos` are required for every agent.
# The field `vel` is required for using [`move_agent!`](@ref) in `ContinuousSpace`.
# `speed` defines how far the bird travels in the direction defined by `vel` per `step`.
# `seperation` defines the minimum distance a bird must maintain from its neighbors.
# `visual_distance` refers to the distance a bird can see and defines a radius of neighboring birds.
# The contribution of each rule defined above recieves an importance weight: `cohere_factor`
# is the importance of maintaining the average position of neighbors,
# `match_factor` is the importance of matching the average trajectory of neighboring birds,
# and `seperate_factor` is the importance of maining the minimum
# distance from neighboring birds.

# The function `initialize_model` generates birds and returns a model object using default values.
function initialize_model(;n_birds=100, speed=1.0, cohere_factor=.25, separation=4.0,
    seperate_factor=.25, match_factor=.01, visual_distance=5.0, dims=(100,100))
    space2d = ContinuousSpace(2; periodic=true, extend=dims)
    model = ABM(Bird, space2d, scheduler=random_activation)
    for _ in 1:n_birds
        vel = Tuple(rand(2)*2 .- 1)
        add_agent!(
            model, vel, speed, cohere_factor,separation, seperate_factor,
            match_factor,visual_distance
        )
    end
    index!(model)
    return model
end

# ## Defining the agent_step!
# `agent_step!` is the primary function called for each step and computes velocity
# according to the three rules defined above.
function agent_step!(bird, model)
    ## Obtain the ids of neibhors within the bird's visual distance
    ids = space_neighbors(bird, model, bird.visual_distance)
    ## Compute velocity based on rules defined above
    bird.vel = (bird.vel .+
        cohere(bird, model, ids) .+
        seperate(bird, model, ids) .+
        match(bird, model, ids))./2
    bird.vel = bird.vel ./ norm(bird.vel)
    ## Move bird according to new velocity and speed
    move_agent!(bird, model, bird.speed)
end

distance(a1, a2) = sqrt(sum((a1.pos .- a2.pos).^2))

get_heading(a1, a2) = a1.pos .- a2.pos

# `cohere` computes the average position of neighboring birds, weighted by importance
function cohere(bird, model, ids)
    N = max(length(ids), 1)
    birds = model.agents
    coherence = (0.0,0.0)
    for id in ids
        coherence = coherence .+ get_heading(birds[id], bird)
    end
    return coherence./N.*bird.cohere_factor
end

# `seperate` repells the bird away from neighboring birds
function seperate(bird, model, ids)
    seperation_vec = (0.0,0.0)
    N = max(length(ids), 1)
    birds = model.agents
    for id in ids
        neighbor = birds[id]
        if distance(bird, neighbor) < bird.separation
            seperation_vec = seperation_vec .- get_heading(neighbor, bird)
        end
    end
    return seperation_vec./N.*bird.seperate_factor
end

# `match` computes the average trajectory of neighboring birds, weighted by importance
function match(bird, model, ids)
    match_vector = (0.0, 0.0)
    N = max(length(ids), 1)
    birds = model.agents
    for id in ids
        match_vector = match_vector .+ birds[id].vel
    end
    return match_vector./N.*bird.match_factor
end

# ## Running the model
Random.seed!(23182)
n_steps = 500
model = initialize_model()
step!(model, agent_step!, n_steps)

# ## Plotting the birds
# The great thing about [`plotabm`](@ref) is its flexibility. We can incorporate the
# direction of the birds when plotting them, by making the "marker" function `am`
# create a `Shape`: a triangle with same orientation as the bird's velocity.
# It is as simple as defining the following function:
function bird_triangle(b::Bird)
    φ = atan(b.vel[2], b.vel[1])
    xs = [(i ∈ (0, 3) ? 2 : 1)*cos(i*2π/3 + φ) for i in 0:3]
    ys = [(i ∈ (0, 3) ? 2 : 1)*sin(i*2π/3 + φ) for i in 0:3]
    Shape(xs, ys)
end

# And here is the animation
# %% #src
using AgentsPlots, Plots
Random.seed!(23182)
cd(@__DIR__) #src
model = initialize_model()
anim = @animate for i in 0:100
    i>0 && step!(model, agent_step!, 1)
    p1 = plotabm(model; am = bird_triangle, as = 10)
    title!(p1, "step $(i)")
end
gif(anim, "flock.gif", fps=30)


text:
export PointMass, getpointmass

import Distributions: mean, var, cov, std, insupport, pdf, logpdf, entropy
import Base: ndims, precision, getindex, convert, isapprox, eltype
import SpecialFunctions: loggamma, logbeta

struct PointMass{P}
    point :: P
end

variate_form(::PointMass{T})  where { T <: Real }　                = Univariate
variate_form(::PointMass{V})  where { T, V <: AbstractVector{T} }　= Multivariate
variate_form(::PointMass{M})  where { T, M <: AbstractMatrix{T} }　= Matrixvariate

##

getpointmass(distribution::PointMass) = distribution.point

##

Base.getindex(distribution::PointMass, index...) = Base.getindex(getpointmass(distribution), index...)

Distributions.entropy(distribution::PointMass) = InfCountingReal(eltype(distribution), -1)

# Real-based univariate point mass

Distributions.insupport(distribution::PointMass{T}, x::Real) where { T <: Real } = x == getpointmass(distribution)
Distributions.pdf(distribution::PointMass{T}, x::Real)       where { T <: Real } = Distributions.insupport(distribution, x) ? one(T) : zero(T)
Distributions.logpdf(distribution::PointMass{T}, x::Real)    where { T <: Real } = Distributions.insupport(distribution, x) ? zero(T) : convert(T, -Inf)

Distributions.mean(distribution::PointMass{T}) where { T <: Real } = getpointmass(distribution)
Distributions.var(distribution::PointMass{T})  where { T <: Real } = zero(T)
Distributions.std(distribution::PointMass{T})  where { T <: Real } = zero(T)
Distributions.cov(distribution::PointMass{T})  where { T <: Real } = zero(T)

probvec(distribution::PointMass{T}) where { T <: Real } = error("probvec(::PointMass{ <: Real }) is not defined")

mean(fn::F, distribution::PointMass{T}) where { T <: Real, F <: Function } = fn(mean(distribution))

Base.precision(::PointMass{T}) where { T <: Real } = convert(T, Inf)
Base.ndims(::PointMass{T})     where { T <: Real } = 1

convert_eltype(::Type{ PointMass }, ::Type{T}, distribution::PointMass{R}) where { T <: Real, R <: Real } = PointMass(convert(T, getpointmass(distribution)))

Base.eltype(::PointMass{T}) where { T <: Real } = T

# AbstractVector-based multivariate point mass

Distributions.insupport(distribution::PointMass{V}, x::AbstractVector) where { T <: Real, V <: AbstractVector{T} } = x == getpointmass(distribution)
Distributions.pdf(distribution::PointMass{V}, x::AbstractVector)       where { T <: Real, V <: AbstractVector{T} } = Distributions.insupport(distribution, x) ? one(T) : zero(T)
Distributions.logpdf(distribution::PointMass{V}, x::AbstractVector)    where { T <: Real, V <: AbstractVector{T} } = Distributions.insupport(distribution, x) ? zero(T) : convert(T, -Inf)

Distributions.mean(distribution::PointMass{V}) where { T <: Real, V <: AbstractVector{T} } = getpointmass(distribution)
Distributions.var(distribution::PointMass{V})  where { T <: Real, V <: AbstractVector{T} } = zeros(T, (ndims(distribution), ))
Distributions.std(distribution::PointMass{V})  where { T <: Real, V <: AbstractVector{T} } = zeros(T, (ndims(distribution), ))
Distributions.cov(distribution::PointMass{V})  where { T <: Real, V <: AbstractVector{T} } = zeros(T, (ndims(distribution), ndims(distribution)))

probvec(distribution::PointMass{V}) where { T <: Real, V <: AbstractVector{T} } = mean(distribution)

mean(::typeof(inv), distribution::PointMass{V})       where { T <: Real, V <: AbstractVector{T} } = error("mean of inverse of `::PointMass{ <: AbstractVector }` is not defined")
mean(::typeof(log), distribution::PointMass{V})       where { T <: Real, V <: AbstractVector{T} } = log.(mean(distribution))
mean(::typeof(mirrorlog), distribution::PointMass{V}) where { T <: Real, V <: AbstractVector{T} } = error("mean of mirrorlog of `::PointMass{ <: AbstractVector }` is not defined")
mean(::typeof(loggamma), distribution::PointMass{V})  where { T <: Real, V <: AbstractVector{T} } = loggamma.(mean(distribution))
mean(::typeof(logdet), distribution::PointMass{V})    where { T <: Real, V <: AbstractVector{T} } = error("mean of logdet of `::PointMass{ <: AbstractVector }` is not defined")

Base.precision(distribution::PointMass{V}) where { T <: Real, V <: AbstractVector{T} } = one(T) ./ cov(distribution)
Base.ndims(distribution::PointMass{V})     where { T <: Real, V <: AbstractVector{T} } = length(mean(distribution))

convert_eltype(::Type{ PointMass }, ::Type{T}, distribution::PointMass{R}) where { T <: Real, R <: AbstractVector }           = PointMass(convert(AbstractVector{T}, getpointmass(distribution)))
convert_eltype(::Type{ PointMass }, ::Type{T}, distribution::PointMass{R}) where { T <: AbstractVector, R <: AbstractVector } = PointMass(convert(T, getpointmass(distribution)))

Base.eltype(::PointMass{V}) where { T <: Real, V <: AbstractVector{T} } = T

# AbstractMatrix-based matrixvariate point mass

Distributions.insupport(distribution::PointMass{M}, x::AbstractMatrix) where { T <: Real, M <: AbstractMatrix{T} } = x == getpointmass(distribution)
Distributions.pdf(distribution::PointMass{M}, x::AbstractMatrix)       where { T <: Real, M <: AbstractMatrix{T} } = Distributions.insupport(distribution, x) ? one(T) : zero(T)
Distributions.logpdf(distribution::PointMass{M}, x::AbstractMatrix)    where { T <: Real, M <: AbstractMatrix{T} } = Distributions.insupport(distribution, x) ? zero(T) : convert(T, -Inf)

Distributions.mean(distribution::PointMass{M}) where { T <: Real, M <: AbstractMatrix{T} } = getpointmass(distribution)
Distributions.var(distribution::PointMass{M})  where { T <: Real, M <: AbstractMatrix{T} } = zeros(T, ndims(distribution))
Distributions.std(distribution::PointMass{M})  where { T <: Real, M <: AbstractMatrix{T} } = zeros(T, ndims(distribution))
Distributions.cov(distribution::PointMass{M})  where { T <: Real, M <: AbstractMatrix{T} } = error("Distributions.cov(::PointMass{ <: AbstractMatrix }) is not defined")

probvec(distribution::PointMass{M}) where { T <: Real, M <: AbstractMatrix{T} } = error("probvec(::PointMass{ <: AbstractMatrix }) is not defined")

mean(::typeof(inv), distribution::PointMass{M})       where { T <: Real, M <: AbstractMatrix{T} } = cholinv(mean(distribution))
mean(::typeof(log), distribution::PointMass{M})       where { T <: Real, M <: AbstractMatrix{T} } = log.(mean(distribution))
mean(::typeof(mirrorlog), distribution::PointMass{M}) where { T <: Real, M <: AbstractMatrix{T} } = error("mean of mirrorlog of `::PointMass{ <: AbstractMatrix }` is not defined")
mean(::typeof(loggamma), distribution::PointMass{M})  where { T <: Real, M <: AbstractMatrix{T} } = loggamma.(mean(distribution))
mean(::typeof(logdet), distribution::PointMass{M})    where { T <: Real, M <: AbstractMatrix{T} } = logdet(mean(distribution))

Base.precision(distribution::PointMass{M}) where { T <: Real, M <: AbstractMatrix{T} } = one(T) ./ cov(distribution)
Base.ndims(distribution::PointMass{M})     where { T <: Real, M <: AbstractMatrix{T} } = size(mean(distribution))

convert_eltype(::Type{ PointMass }, ::Type{T}, distribution::PointMass{R}) where { T <: Real, R <: AbstractMatrix }           = PointMass(convert(AbstractMatrix{T}, getpointmass(distribution)))
convert_eltype(::Type{ PointMass }, ::Type{T}, distribution::PointMass{R}) where { T <: AbstractMatrix, R <: AbstractMatrix } = PointMass(convert(T, getpointmass(distribution)))

Base.eltype(::PointMass{M}) where { T <: Real, M <: AbstractMatrix{T} } = T

Base.isapprox(left::PointMass, right::PointMass; kwargs...) = Base.isapprox(getpointmass(left), getpointmass(right); kwargs...)
Base.isapprox(left::PointMass, right; kwargs...) = false
Base.isapprox(left, right::PointMass; kwargs...) = false


text:
"""
    DAMMviz()

Interactive plot of the DAMM model

# To open
```julia-repl
julia> DAMMviz()
```
"""
function DAMMviz(;width = 2200, height = 1600, fontsize = 30)
  fontsize_theme = Theme(fontsize = fontsize, font = "JuliaMono")
  set_theme!(fontsize_theme)

  fig = Figure(resolution = (width, height), figure_padding = 30)
  ax3D = Axis3(fig[1, 2])

  texts = Array{Label}(undef, 9);
  sliderranges = [
    1e8:1e8:1e9, # α or p[1]
    58.0:1.0:70.0, # Ea p[2]
    1e-9:1e-6:1e-5, # kMsx p[3]
    1e-4:1e-3:1e-2, # kMo2 p[4]
    0.3:0.05:0.8, # porosity p[5]
    0.01:0.01:0.1, # sx, p[6]
    0:2:40, # Ts
    0:0.02:0.8, # θ
    0.7:0.1:3.0 # Q10Km
   ]; #
  sliders = [Slider(fig, range = sr) for sr in sliderranges];
  texts[1] = Label(fig, text= lift(X->string("Parameters\n", to_latex("\\alpha_{sx}"), " = ", X, to_latex(" (mgC cm^{-3} h^{-1})")), sliders[1].value), justification = :left, halign = :left, textsize = fontsize, width = Auto(false));
  texts[2] = Label(fig, text= lift(X->string(to_latex("E_a"), " = ", X, to_latex(" (kJ mol^{-1})")), sliders[2].value), justification = :left, halign = :left, textsize = fontsize, width = Auto(false));
  texts[3] = Label(fig, text= lift(X->string(to_latex("kM_{sx}"), " = ", round(X, sigdigits = 2), to_latex(" (gC cm^{-3})")), sliders[3].value), justification = :left, halign = :left, textsize = fontsize, width = Auto(false));
  texts[4] = Label(fig, text= lift(X->string(to_latex("kM_{o2}"), " = ", X, to_latex(" (L L^{-1})")), sliders[4].value), justification = :left, halign = :left, textsize = fontsize, width = Auto(false));
  texts[5] = Label(fig, text= lift(X->string(to_latex("Porosity"), " = ", X, to_latex(" (m^3 m^{-3})")), sliders[5].value), justification = :left, halign = :left, textsize = fontsize, width = Auto(false));
  texts[6] = Label(fig, text= lift(X->string(to_latex("S_x"), " = ", X, to_latex(" (gC cm^{-3})")), sliders[6].value), justification = :left, halign = :left, textsize = fontsize, width = Auto(false));
  texts[7] = Label(fig, text= lift(X->string("Drivers\n", to_latex("T_s"), " = ", X, to_latex(" (°C)")), sliders[7].value), justification = :left, halign = :left, textsize = fontsize, color = :green, width = Auto(false));
  texts[8] = Label(fig, text= lift(X->string(to_latex("θ"), " = ", round(X, sigdigits = 3), to_latex(" (m^3 m^{-3})")), sliders[8].value), justification = :left, halign = :left, textsize = fontsize, color = :green, width = Auto(false));
  texts[9] = Label(fig, text= lift(X->string(to_latex("Q_{10} of kM_{sx}"), " = ", X, to_latex(" (-)")), sliders[9].value), justification = :left, halign = :left, textsize = fontsize, width = Auto(false));

  αsx = sliders[1].value
  set_close_to!(sliders[1], 7e8)
  Ea = sliders[2].value
  set_close_to!(sliders[2], 61)
  kMsx = sliders[3].value
  set_close_to!(sliders[3], 5e-6)
  kMo2 = sliders[4].value
  set_close_to!(sliders[4], 0.004)
  porosity = sliders[5].value
  set_close_to!(sliders[5], 0.7)
  sx = sliders[6].value
  set_close_to!(sliders[6], 0.05)
  Ts = sliders[7].value
  set_close_to!(sliders[7], 28)
  θ = sliders[8].value
  set_close_to!(sliders[8], 0.4)
  Q10Km = sliders[9].value
  set_close_to!(sliders[9], 1.0)

  r = 50
  x = collect(range(0, length=r, stop=40)) # T axis, °C from min to max
  y = @lift(collect(range(0, length=r, stop=$porosity))) # M axis, % from min to max
  X = repeat(1:r, inner=r) # X for DAMM matrix
  Y = repeat(1:r, outer=r) # Y for DAMM matrix
  X2 = repeat(x, inner=r) # T values to fit DAMM on
  Y2 = @lift(repeat($y, outer=r)) # M values to fit DAMM on
  xy = @lift(hcat(X2, $Y2)) # T and M matrix to create DAMM matrix

  params = @lift(($αsx, $Ea, $kMsx, $kMo2, $porosity, $sx, $Q10Km))
  DAMM_Matrix = @lift(Matrix(sparse(X, Y, DAMM($xy, $params))))
  Rₛ = @lift(DAMM(hcat($Ts, $θ), $params))
  Rₛᵣ = @lift(round.($Rₛ, sigdigits = 3)[1])

  vertical_sublayout = fig[1, 1] = hgrid!(vgrid!(
  Iterators.flatten(zip(texts[vcat(1:6, 9)], sliders[vcat(1:6, 9)]))...;
   ),  #width = 200, height = 1000);
  #vertical_sublayout2 = fig[1, 2] =
  vgrid!(
  Iterators.flatten(zip(texts[7:8], sliders[7:8]))...,
  Label(fig, text = lift(X -> string("\n", to_latex("R_{s} = "), X, to_latex(" (\\mumol m^{-2} s^{-1})")), Rₛᵣ), color = :red, justification = :left, halign = :left, textsize = fontsize, width = Auto(false))), halign = :left) #width = 200, height = 1000);

  s3D = surface!(ax3D, x, y, DAMM_Matrix, colormap = Reverse(:Spectral),
	transparency = true, alpha = 0.2, shading = false, colorrange = (0, 30))
  w3D = wireframe!(ax3D, x, y, DAMM_Matrix, overdraw = true,
	transparency = true, color = (:black, 0.1));
  point = @lift(DAMM(hcat($Ts, $θ), $params))
  point3D = @lift(Vec3f.($Ts, $θ, $point))
  scatter3D = scatter!(ax3D, point3D, markersize = 20, color = :black)
  cb = Colorbar(fig[1, 3], colormap = Reverse(:Spectral), limits = (0, 30),
	   label = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})"));

  ax2D = Axis(fig[2, 1])
  isoθ = @lift(collect(range($θ, length=length(x), stop=$θ)))
  isoy = @lift(DAMM(hcat(x, $isoθ), $params))
  lines!(ax2D, x, isoy, color = isoy, linewidth = 8,
	 colormap = Reverse(:Spectral), colorrange = (0, 30))
  pointTs2D = @lift(Point2f.($Ts, $point))
  scatter!(ax2D, pointTs2D, color = :black, markersize = 20)
  linep = @lift(repeat($point, r))
  linev = @lift(repeat([$θ], r))
  lineh = @lift(repeat([$Ts], r))
  allθ = collect(range(0, length=r, stop=1.0))
  allR = collect(range(0, length=r, stop=30.0))
  lines!(ax2D, x, linep, color = (:black, 0.4), linestyle = :dash)
  lines!(ax2D, lineh, x, color = (:black, 0.4), linestyle = :dash)

  ax2D2 = Axis(fig[2, 2:3])
  isoT = @lift(collect(range($Ts, length=length(x), stop=$Ts)))
  isox = @lift(DAMM(hcat($isoT, $y), $params))
  lines!(ax2D2, y, isox, color = isox, linewidth = 8,
	 colormap = Reverse(:Spectral), colorrange = (0, 30))
  pointθ2D = @lift(Point2f.($θ, $point))
  scatter!(ax2D2, pointθ2D, color = :black, markersize = 20)
  lines!(ax2D2, allθ, linep, color = (:black, 0.4), linestyle = :dash)
  lines!(ax2D2, linev, x, color = (:black, 0.4), linestyle = :dash)

  # isoline in the 3D figure
  lines!(ax3D, x, isoθ, isoy, color = isoy, linewidth = 8,
	 colormap = Reverse(:Spectral), colorrange = (0, 30))
  # lines!(ax3D, x .+1, isoθ, isoy, linewidth = 1, color = :black)
  # lines!(ax3D, x .-1, isoθ, isoy, linewidth = 1, color = :black)
  lines!(ax3D, isoT, y, isox, color = isox, linewidth = 8,
	 colormap = Reverse(:Spectral), colorrange = (0, 30))
  iso40 = repeat([40], r)
  iso1 = repeat([1], r)
  iso0 = repeat([0], r)
  lines!(ax3D, x, iso1, linep,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, iso40, allθ, linep,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, x, isoθ, iso0,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, iso40, isoθ, allR,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, isoT, allθ, iso0,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, isoT, iso1, allR,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, x, isoθ, linep,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, isoT, allθ, linep,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)
  lines!(ax3D, isoT, isoθ, allR,
	 color = (:black, 0.4), linestyle = :dash, transparency = true)

  lines!(ax3D, x, iso1, isoy, color = isoy, linewidth = 8, alpha = 0.2,
	 colormap = Reverse(:Spectral), colorrange = (0, 30), transparency = true)
  lines!(ax3D, iso40, y, isox, color = isox, linewidth = 8, alpha = 0.2,
	 colormap = Reverse(:Spectral), colorrange = (0, 30), transparency = true)

  point3D_T = @lift(Vec3f.($Ts, 1.0, $point))
  point3D_θ = @lift(Vec3f.(40.0, $θ, $point))
  scatter3D_T = scatter!(ax3D, point3D_T, markersize = 20, color = :black)
  scatter3D_θ = scatter!(ax3D, point3D_θ, markersize = 20, color = :black)

  ylims!(ax2D, 0.0, 30.0); xlims!(ax2D, 10.0, 40.0);
  ylims!(ax2D2, 0.0, 30.0); xlims!(ax2D2, 0.0, 1.0);
  xlims!(ax3D, 0.0, 40.0);
  ylims!(ax3D, 0.0, 1.0);
  zlims!(ax3D, 0.0, 30.0);

  ax3D.xlabel = to_latex("T_{soil} (°C)");
  ax3D.ylabel = to_latex("\\theta (m^3 m^{-3})");
  ax3D.zlabel = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})");
  ax3D.zlabeloffset = 80
  ax3D.xlabeloffset = 50
  ax3D.ylabeloffset = 60

  ax2D.xlabel = to_latex("T_{soil} (°C)");
  ax2D.ylabel = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})");
  # ax2D2.ylabel = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})");
  hideydecorations!(ax2D2, ticks = false)
  ax2D2.xlabel = to_latex("\\theta (m^3 m^{-3})");

  #FZ = 30; ax2D.xlabelsize = FZ; ax2D.ylabelsize = FZ; ax2D2.xlabelsize = FZ; ax2D2.ylabelsize = FZ;
  #ax2D.xticklabelsize = FZ; ax2D.yticklabelsize = FZ; ax2D2.xticklabelsize = FZ; ax2D2.yticklabelsize = FZ;
  colsize!(fig.layout, 1, Relative(1/2))
  #colsize!(fig.layout, 1, Aspect(1, 2.0))
  rowsize!(fig.layout, 1, Relative(1/2))
  #resize_to_layout!(fig)

  supertitle = Label(fig[0, :], "Dual Arrhenius and Michaelis-Menten (DAMM) interactive visualisation, v0.1.3", textsize = fontsize*4/3)

  cb.alignmode = Mixed(right = 0)
  #set_theme!(figure_padding = 30)
  DataInspector(fig)
  fig
  return fig
end


text:
_length(v::AbstractArray) = length(v)
_length(::Any) = 1
_map(f, v::AbstractArray) = map(f, v)
_map(f, v) = f(v)
function _searchsortedfirst(vals, t)
    rev = first(vals) > last(vals)
    searchsortedfirst(vals, t, rev = rev)
end

function format(x)
    io = IOBuffer()
    show(IOContext(io, :compact => true), x)
    String(take!(io))
end

for func in [:rangeslider, :slider]
    @eval begin
        function $func(WT::WidgetTheme, vals::AbstractArray, formatted_vals = format.(vec(vals)); value = medianelement(vals), kwargs...)

            T = Observables.to_value(value) isa Vector ? Vector{eltype(vals)} : eltype(vals)
            value isa AbstractObservable || (value = Observable{T}(value))

            vals = vec(vals)
            indices = axes(vals)[1]
            f = x -> _map(t -> _searchsortedfirst(vals, t), x)
            g = x -> vals[Int.(x)]
            index = ObservablePair(value, f = f, g = g).second
            wdg = Widget($func(WT, indices, formatted_vals; value = index, kwargs...), output = value)
            wdg["value"] = value
            wdg
        end
    end
end

"""
```
function slider(vals::AbstractArray;
                value=medianelement(vals),
                label=nothing, readout=true, kwargs...)
```

Creates a slider widget which can take on the values in `vals`, and updates
observable `value` when the slider is changed.
"""
function slider(::WidgetTheme, vals::AbstractUnitRange{<:Integer}, formatted_vals = format.(vals);
    className=getclass(:input, "range", "fullwidth"),
    readout=true, label=nothing, value=medianelement(vals), orientation = "horizontal", attributes = Dict(), kwargs...)

    min, max = extrema(vals)
    orientation = string(orientation)
    attributes = merge(attributes, Dict("orient" => orientation))
    (value isa AbstractObservable) || (value = convert(eltype(vals), value))
    format = js"""
        function(){
            return this.formatted_vals()[parseInt(this.index())-($min)];
        }
    """
    ui = input(value; bindto="index", attributes=attributes, extra_obs = ["formatted_vals" => formatted_vals], computed = ["formatted_val" => format],
               typ="range", min=min, max=max, step=1, className=className, kwargs...)
    if (label != nothing) || readout
        if orientation != "vertical"
            Widgets.scope(ui).dom = readout ?
                flex_row(wdglabel(label), Widgets.scope(ui).dom, node(:p, attributes = Dict("data-bind" => "text: formatted_val"))) :
                flex_row(wdglabel(label), Widgets.scope(ui).dom)
        else
            readout && (label = vbox(label, node(:p, attributes = Dict("data-bind" => "text: formatted_val"))))
            Widgets.scope(ui).dom  = hbox(wdglabel(label), dom"div[style=flex-shrink:1]"(Widgets.scope(ui).dom))
        end
    end
    Widget{:slider}(ui)
end

"""
```
function rangeslider(vals::AbstractArray;
                value=medianelement(vals),
                label=nothing, readout=true, kwargs...)
```

Creates a slider widget which can take on the values in `vals` and accepts several "handles".
Pass a vector to `value` with two values if you want to select a range.
"""
function rangeslider(theme::WidgetTheme, vals::AbstractUnitRange{<:Integer}, formatted_vals = format.(vals);
    style = Dict(), label = nothing, value = medianelement(vals), orientation = "horizontal", readout = true,
    className = "is-primary")

    T = Observables.to_value(value) isa Vector ? Vector{eltype(vals)} : eltype(vals)
    value isa AbstractObservable || (value = Observable{T}(value))

    index = value
    orientation = string(orientation)
    preprocess = T<:Vector ? js"unencoded.map(Math.round)" : js"Math.round(unencoded[0])"

    scp = Scope(imports = vcat([nouislider_min_js, nouislider_min_css], libraries(theme)))
    setobservable!(scp, "index", index)
    fromJS = Observable(scp, "fromJS", false)
    changes = Observable(scp, "changes", 0)
    connect = _length(index[]) > 1 ? js"true" : js"[true, false]"
    min, max = extrema(vals)
    s = step(vals)

    id = "slider"*randstring()
    start = JSExpr.@js $index[]
    updateValue = JSExpr.@js function updateValue(values, handle, unencoded, tap, positions)
        $fromJS[] = true
        $index[] = $preprocess
    end
    updateCount = JSExpr.@js function updateCount(values, handle, unencoded, tap, positions)
        $changes[] = $changes[]+1
    end
    tooltips = JSString("[" * join(fill(readout, _length(value[])), ", ") * "]")

    onimport(scp, js"""
        function (noUiSlider) {
            var vals = JSON.parse($(JSON.json(formatted_vals)));
            $updateValue
            $updateCount
            var slider = document.getElementById($id);
            noUiSlider.create(slider, {
            	start: $start,
                step: 1,
                tooltips: $tooltips,
                connect: $connect,
                orientation: $orientation,
                format: {
                    to: function ( value ) {
                        var ind = Math.round(value-($min));
                        return ind + 1 > vals.length ? vals[vals.length - 1] : vals[ind];
                    },
                    from: function ( value ) {
                        return parseInt(value);
                    }
                },
            	range: {
                        'min': ($min),
                        'max': ($max)
            	},})

            slider.noUiSlider.on("slide", updateValue);
            slider.noUiSlider.on("change", updateCount);
        }
        """)
    slap_design!(scp)
    onjs(index, @js function (val)
        if !$fromJS[]
            document.getElementById($id).noUiSlider.set(Array.isArray(val) ? val : [val])
        end
        $fromJS[] = false
    end)

    style = Dict{String, Any}(string(key) => val for (key, val) in style)

    haskey(style, "flex-grow") || (style["flex-grow"] = "1")
    !haskey(style, "height") && orientation == "vertical" && (style["height"] = "20em")
    scp.dom = node(:div, style = style, attributes = Dict("id" => id))
    layout = function (t)
        if orientation != "vertical"
            sld = t.scope
            sld = label !== nothing ?  flex_row(label, sld) : sld
            sld = readout ? vbox(vskip(3em), sld) : sld
            sld = div(sld, className = "field rangeslider rangeslider-horizontal interact-widget $className")
        else
            sld = t.scope
            sld = readout ? hbox(hskip(6em), sld) : sld
            sld = label !== nothing ?  vbox(label, sld) : sld
            sld = div(sld, className = "field rangeslider rangeslider-vertical interact-widget $className")
        end
        sld
    end
    Widget{:rangeslider}(["index" => index, "changes" => changes];
        scope = scp, output = value, layout = layout)
end

"""
```
function rangepicker(vals::AbstractArray;
                value=[extrema(vals)...],
                label=nothing, readout=true, kwargs...)
```

A multihandle slider with a set of spinboxes, one per handle.
"""
function rangepicker(::WidgetTheme, vals::AbstractRange{S}; value = [extrema(vals)...], readout = false, className = "is-primary") where {S}
    T = Observables.to_value(value) isa Vector ? Vector{eltype(vals)} : eltype(vals)
    value isa AbstractObservable || (value = Observable{T}(value))
    wdg = Widget{:rangepicker}(output = value)
    if !(T<:Vector)
        wdg["input"] = input(S, vals, value=value)
    else
        function newinput(i)
            f = t -> t[i]
            g = t -> (s = copy(value[]); s[i] = t; s)
            new_val = ObservablePair(value, f=f, g=g).second
            input(S, vals, value = new_val)
        end

        for i in eachindex(value[])
            wdg["input$i"] = newinput(i)
        end
    end
    inputs = t -> (val for (key, val) in components(t) if occursin(r"slider|input", string(key)))
    wdg.layout = t -> div(inputs(t)..., className = "interact-widget")
    wdg["slider"] = rangeslider(vals, value = value, readout = readout, className = className)
    wdg["changes"] = map(+, (val["changes"] for val in inputs(wdg))...)
    return wdg
end


text:
using Distributions

@enum DogState begin
    ACTIVE = 1
    OCCUPIED
    SLEEPY
end # enum
Base.length( ::DogState ) = 1
Base.iterate( ds::DogState, state = 1 ) = nothing

Base.@kwdef struct Dog{F<:Number}
    state::DogState = ACTIVE
    sleepyness::UInt32 = 0
    position::Point2{F}
    velocity::Point2{F}
end # struct

function Dog( position, velocity )
    T = eltype(position)
    @assert T == eltype(velocity)
    Dog{T}(;
        position = convert(Point2{T}, position),
        velocity = convert(Point2{T}, velocity),
        )
end # function

#TODO: this is slow and allocating
function Dog( v::D; kwargs... ) where D<:Dog
    nt = NamedTuple()
    for property in propertynames(v)
        if property in keys(kwargs)
            nt = merge( nt, NamedTuple{(property,)}((kwargs[property],)) )
            continue
        end
        value = getproperty( v, property )
        # if value isa Vector
        #     nt = merge( nt, NamedTuple{(property,)}(Ref(copy(value))) )
        #     continue
        # end # if
        nt = merge( nt, NamedTuple{(property,)}((value,)) )
    end
    return D(;nt...)
end # function


text:
################################################################################
#
#  AlgAssAbsOrd / AlgAssAbsOrdElem
#
################################################################################

# Orders in algebras over the rationals
@attributes mutable struct AlgAssAbsOrd{S, T} <: Ring
  algebra::S                       # Algebra containing the order
  dim::Int
  basis#::Vector{AlgAssAbsOrdElem{S, T}}
  basis_alg::Vector{T}             # Basis as array of elements of the algebra
  basis_matrix::FakeFmpqMat           # Basis matrix of order wrt basis of the algebra
  basis_mat_inv::FakeFmpqMat       # Inverse of basis matrix
  gen_index::fmpq                  # The det of basis_mat_inv as fmpq
  index::fmpz                      # The det of basis_mat_inv
                                   # (this is the index of the equation order
                                   #  in the given order)
  disc::fmpz                       # Discriminant

  ismaximal::Int                   # 0 Not known
                                   # 1 Known to be maximal
                                   # 2 Known to not be maximal

  #trace_mat::fmpz_mat              # The reduced trace matrix (if known)
  trred_matrix::fmpz_mat

  picard_group#::MapPicardGrp

  tcontain::FakeFmpqMat

  isnice::Bool
  nice_order#Tuple{AlgAssAbsOrd, T}
  nice_order_ideal::fmpz

  function AlgAssAbsOrd{S, T}(A::S) where {S, T}
    # "Default" constructor with default values.
    O = new{S, T}(A, dim(A))
    O.ismaximal = 0
    O.isnice = false
    O.tcontain = FakeFmpqMat(zero_matrix(FlintZZ, 1, dim(A)))
    return O
  end

  function AlgAssAbsOrd{S, T}(A::S, M::FakeFmpqMat, Minv::FakeFmpqMat, B::Vector{T}, cached::Bool = false) where {S, T}
    return get_cached!(AlgAssAbsOrdID, (A, M), cached) do
      O = AlgAssAbsOrd{S, T}(A)
      O.basis_alg = B
      O.basis_matrix = M
      O.basis_mat_inv = Minv
      return O
    end::AlgAssAbsOrd{S, T}
  end

  function AlgAssAbsOrd{S, T}(A::S, M::FakeFmpqMat, cached::Bool = false) where {S, T}
    return get_cached!(AlgAssAbsOrdID, (A, M), cached) do
      O = AlgAssAbsOrd{S, T}(A)
      d = dim(A)
      O.basis_matrix = M
      O.basis_alg = Vector{T}(undef, d)
      for i in 1:d
        O.basis_alg[i] = elem_from_mat_row(A, M.num, i, M.den)
      end
      return O
    end::AlgAssAbsOrd{S, T}
  end

  function AlgAssAbsOrd{S, T}(A::S, B::Vector{T}, cached::Bool = false) where {S, T}
    M = basis_matrix(B, FakeFmpqMat)
    return get_cached!(AlgAssAbsOrdID, (A, M), cached) do
      O = AlgAssAbsOrd{S, T}(A)
      O.basis_alg = B
      O.basis_matrix = M
      return O
    end::AlgAssAbsOrd{S, T}
  end
end

const AlgAssAbsOrdID = Dict{Tuple{AbsAlgAss, FakeFmpqMat}, AlgAssAbsOrd}()

mutable struct AlgAssAbsOrdElem{S, T} <: RingElem
  elem_in_algebra::T
  coordinates::Vector{fmpz}
  has_coord::Bool # needed for mul!
  parent::AlgAssAbsOrd{S, T}

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}) where {S, T}
    z = new{S, T}()
    z.parent = O
    z.elem_in_algebra = algebra(O)()
    z.coordinates = Vector{fmpz}(undef, degree(O))
    z.has_coord = false
    return z
  end

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}, a::T) where {S, T}
    z = new{S, T}()
    z.elem_in_algebra = a
    z.parent = O
    z.coordinates = Vector{fmpz}(undef, degree(O))
    z.has_coord = false
    return z
  end

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}, arr::Vector{fmpz}) where {S, T}
    z = new{S, T}()
    z.elem_in_algebra = dot(O.basis_alg, arr)
    z.coordinates = arr
    z.parent = O
    z.has_coord = true
    return z
  end

  function AlgAssAbsOrdElem{S, T}(O::AlgAssAbsOrd{S, T}, a::T, arr::Vector{fmpz}) where {S, T}
    z = new{S, T}()
    z.parent = O
    z.elem_in_algebra = a
    z.coordinates = arr
    z.has_coord = true
    return z
  end
end

################################################################################
#
#  AlgAssAbsOrdIdl
#
################################################################################

mutable struct AlgAssAbsOrdIdl{S, T}
  algebra::S

  basis::Vector{T} # Basis of the ideal as array of elements of the algebra
  # The basis matrix is in the BASIS of the ALGEBRA!
  basis_matrix::FakeFmpqMat
  basis_mat_inv::FakeFmpqMat

  # Basis matrices with respect to orders
  basis_matrix_wrt::Dict{AlgAssAbsOrd{S, T}, FakeFmpqMat}

  # Left and right order:
  # The largest orders of which the ideal is a left resp. right ideal.
  left_order::AlgAssAbsOrd{S, T}
  right_order::AlgAssAbsOrd{S, T}

  # Any order contained in the left or right order, that is, an order of which
  # the ideal is a (possibly fractional) ideal.
  order::AlgAssAbsOrd{S, T}
  gens::Vector{T}    # Generators of the ideal w. r. t. order

  # isleft and isright with respect to `order`
  isleft::Int                      # 0 Not known
                                   # 1 Known to be a left ideal
                                   # 2 Known not to be a left ideal
  isright::Int                     # as for isleft

  iszero::Int                      # 0: don't know, 1: known to be zero, 2: known to be not zero

  norm::Dict{AlgAssAbsOrd{S, T}, fmpq} # The ideal has different norms with respect
                                       # to different orders
  normred::Dict{AlgAssAbsOrd{S, T}, fmpq}

  function AlgAssAbsOrdIdl{S, T}(A::S) where { S <: AbsAlgAss, T <: AbsAlgAssElem }
    r = new{S, T}()
    r.algebra = A
    r.isleft = 0
    r.isright = 0
    r.iszero = 0
    r.basis_matrix_wrt = Dict{AlgAssAbsOrd{S, T}, FakeFmpqMat}()
    r.norm = Dict{AlgAssAbsOrd{S, T}, fmpq}()
    r.normred = Dict{AlgAssAbsOrd{S, T}, fmpq}()
    return r
  end

  function AlgAssAbsOrdIdl{S, T}(A::S, M::FakeFmpqMat) where { S <: AbsAlgAss, T <: AbsAlgAssElem }
    r = AlgAssAbsOrdIdl{S, T}(A)
    r.basis_matrix = M
    return r
  end
end

mutable struct AlgAssAbsOrdIdlSet{S, T}
  order::AlgAssAbsOrd{S, T}

  function AlgAssAbsOrdIdlSet{S, T}(O::AlgAssAbsOrd{S, T}) where {S, T}
    z = new{S, T}(O)
    return z
  end
end

function AlgAssAbsOrdIdlSet(O::AlgAssAbsOrd{S, T}) where {S, T}
  return AlgAssAbsOrdIdlSet{S, T}(O)
end

################################################################################
#
#  AlgAssAbsOrdQuoRing / AlgAssAbsOrdQuoRingElem
#
################################################################################

const AlgAssAbsOrdQuoRing{S, T} = AbsOrdQuoRing{AlgAssAbsOrd{S, T}, AlgAssAbsOrdIdl{S, T}} where {S, T}

const AlgAssAbsOrdQuoRingElem{S, T} = AbsOrdQuoRingElem{AlgAssAbsOrd{S, T}, AlgAssAbsOrdIdl{S, T}, AlgAssAbsOrdElem{S, T}} where {S, T}


text:
using Colors, Reactive, Quaternions, GeometryTypes
using GLFW, GLAbstraction, GLWindow, GLVisualize
import Base: push!, ==

export Signal, RGBA

include("Window.jl")
export init_window, camera

include("KeyBinds.jl")
export Key, Hotkeys
export init_screenshots!, init_QWEASD_controls!, init_angles

include("Real_space.jl")
export plot3D, glsphere


text:
#
# Bond energies taken from
# http://www.wiredchemist.com/chemistry/data/bond_energies_lengths.html
#

# Constants

const a0 = 4π * Unitful.ϵ0 * Unitful.ħ^2 / (Unitful.me * Unitful.q^2) |> Å

const HHssσ = -432 * Unitful.kJ / Unitful.mol / Unitful.Na |> eV

const HCssσ = HHssσ * exp(-(3 * 109u"pm" / 2 - 2 * 74u"pm") / a0)
const HCspσ = (- √3 * 411 * Unitful.kJ / Unitful.mol / Unitful.Na - HCssσ) / √2 |> eV

const CCssσ = HHssσ * exp(-(134u"pm" - 2 * 74u"pm") / a0)
const CCspσ = HCspσ * exp(-(134u"pm" - 3/2 * 109u"pm") / a0)
const CCppσ = (3 * (347 * Unitful.kJ / Unitful.mol / Unitful.Na |> eV) - CCssσ) / 2
const CCppπ = -(602 - 347) * Unitful.kJ / Unitful.mol / Unitful.Na |> eV

##################################################
# Onsite energies
##################################################

@inline onsiteenergy(::Hydrogen, ::Orbital"1s") = -13.59843449eV
@inline onsiteenergy(::Carbon, ::Orbital"2s") = -15.9625933eV
@inline onsiteenergy(::Carbon, ::Orbital"2p_{*}") = -11.2602880eV

##################################################
# Bond energies
##################################################

τ = 0//1

@inline bondenergy(::Val{:σ}, δr, ::Hydrogen, ::Hydrogen, ::Orbital"1s", ::Orbital"1s") = 2*HHssσ * exp(-(δr - 1Å)/a0*τ)

@inline bondenergy(::Val{:σ}, δr, ::Hydrogen, ::Carbon, ::Orbital"1s", ::Orbital"2s") = 2*HCssσ * exp(-8(δr - 1Å)/a0*τ)
@inline bondenergy(::Val{:σ}, δr, ::Carbon, ::Hydrogen, ::Orbital"2s", ::Orbital"1s") = 2*HCssσ * exp(-8(δr - 1Å)/a0*τ)
@inline bondenergy(::Val{:σ}, δr, ::Hydrogen, ::Carbon, ::Orbital"1s", ::Orbital"2p_{*}") = HCspσ * exp(-8(δr - 1Å)/a0*τ)
@inline bondenergy(::Val{:σ}, δr, ::Carbon, ::Hydrogen, ::Orbital"2p_{*}", ::Orbital"1s") = -HCspσ * exp(-8(δr - 1Å)/a0*τ)

@inline bondenergy(::Val{:σ}, δr, ::Carbon, ::Carbon, ::Orbital"2s", ::Orbital"2s") = 2*CCssσ * exp(-12(δr - 1Å)/a0*τ)
@inline bondenergy(::Val{:σ}, δr, ::Carbon, ::Carbon, ::Orbital"2s", ::Orbital"2p_{*}") = 2*CCspσ * exp(-12(δr - 1Å)/a0*τ)
@inline bondenergy(::Val{:σ}, δr, ::Carbon, ::Carbon, ::Orbital"2p_{*}", ::Orbital"2s") = -CCspσ * exp(-12(δr - 1Å)/a0*τ)
@inline bondenergy(::Val{:σ}, δr, ::Carbon, ::Carbon, ::Orbital"2p_{*}", ::Orbital"2p_{*}") = CCppσ * exp(-12(δr - 1Å)/a0*τ)
@inline bondenergy(::Val{:π}, δr, ::Carbon, ::Carbon, ::Orbital"2p_{*}", ::Orbital"2p_{*}") = CCppπ * exp(-12(δr - 1Å)/a0*τ)

##################################################
# Slater-Koster scheme
##################################################

function slaterkoster(f, δr::SVector{3}, atom1::Atom, atom2::Atom, s1::Orbital"*s", s2::Orbital"*s")
    return f(Val(:σ), norm(δr), atom1, atom2, s1, s2)
end

for (i, x) in [(1, :(:x)), (2, :(:y)), (3, :(:z))]
    @eval @inline function slaterkoster(f, δr::SVector{3}, atom1::Atom, atom2::Atom, s1::Orbital"*s", s2::Orbital{<:Any, :p, $x})
        d = normalize(δr)
        return d[$i] * f(Val(:σ), norm(δr), atom1, atom2, s1, s2)
    end

    @eval @inline function slaterkoster(f, δr::SVector{3}, atom1::Atom, atom2::Atom, s1::Orbital{<:Any,:p,$x}, s2::Orbital{<:Any,:s,Nothing})
        d = normalize(δr)
        return d[$i] * f(Val(:σ), norm(δr), atom1, atom2, s1, s2)
    end

    for (j, y) in [(1, :(:x)), (2, :(:y)), (3, :(:z))]
        if i == j
            @eval @inline function slaterkoster(f, δr::SVector{3}, atom1::Atom, atom2::Atom, s1::Orbital{<:Any,:p,$x}, s2::Orbital{<:Any,:p,$x})
                d = normalize(δr)
                return d[$i]^2 * f(Val(:σ), norm(δr), atom1, atom2, s1, s2) + (1 - d[$i]^2) * f(Val(:π), norm(δr), atom1, atom2, s1, s2)
            end
        else
            @eval @inline function slaterkoster(f, δr::SVector{3}, atom1::Atom, atom2::Atom, s1::Orbital{<:Any,:p,$x}, s2::Orbital{<:Any,:p,$y})
                d = normalize(δr)
                return d[$i] * d[$j] * (f(Val(:σ), norm(δr), atom1, atom2, s1, s2) - f(Val(:π), norm(δr), atom1, atom2, s1, s2))
            end
        end
    end
end


text:
# data.jl
# Takes in month of data, returns smaller (reduced) CSV with only relevant information for link time determination
# Authored by Brandon Zeng on 7/7/15

num = 12

using DataFrames, Geodesy, KDTrees
cd("/Users/bzeng/Dropbox (MIT)/7 Coding/UROP/Data")
function reduce(num::Int64)
	println("loading CSV")
	@time df = readtable("trip_data_$(num).csv");
	println("completed loading")
	rm("trip_data_$(num).csv")
	@time df = df[setdiff(names(df), [:medallion, :hack_license, :vendor_id, :rate_code, :store_and_fwd_flag, :passenger_count, :trip_time_in_secs, :trip_distance])];
	@time df = df[~isna(df[:,:dropoff_longitude]),:];

	#polygon representing manhattan
	city = [(40.87851017592601,-73.927903175354);(40.87743937506905,-73.92168045043945);(40.87516792197097,-73.91871929168701);(40.87393481479329,-73.91193866729736);(40.87182549927143,-73.90970706939697);(40.866016421491814,-73.91168117523193);(40.859655059077824,-73.91764640808105);(40.84855366761019,-73.9266586303711);(40.84394377141786,-73.92949104309082);(40.83550228531863,-73.9339542388916);(40.82816381252365,-73.9339542388916);(40.81972031701224,-73.93318176269531);(40.808677200957106,-73.93335342407227);(40.801920506109774,-73.92803192138672);(40.7958778790764,-73.92820358276367);(40.78210123234386,-73.94107818603516);(40.77547182731999,-73.93953323364258);(40.7701418259051,-73.94588470458984);(40.75831029512206,-73.95669937133789);(40.744395800976775,-73.96888732910156);(40.71746884168838,-73.97077560424805);(40.70927151739564,-73.9764404296875);(40.7066689811733,-73.99497985839844);(40.697299008636755,-74.01008605957031);(40.69899090674371,-74.02381896972656);(40.74907763805906,-74.01351928710938);(40.77950154452169,-73.99532318115234);(40.81822635589172,-73.96751403808594);(40.853163243121564,-73.94845962524414)]

	#keep points in the dataset only if inside manhattan
	function point_inside_polygon(x,y,poly)
	    n = length(poly)
	    inside =false

	    p1x,p1y = poly[1]
	    for i in 0:n
	        p2x,p2y = poly[i % n + 1]
	        if y > min(p1y,p2y)
	            if y <= max(p1y,p2y)
	                if x <= max(p1x,p2x)
	                    if p1y != p2y
	                        xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
	                    end
	                    if p1x == p2x || x <= xinters
	                        inside = !inside
	                    end
	                end
	            end
	        end
	        p1x,p1y = p2x,p2y
	    end
	    return inside
	end

	tic()

	mask = falses(nrow(df));
	for i = 1:nrow(df)
		mask[i] = point_inside_polygon(df[i, :pickup_latitude], df[i, :pickup_longitude], city) && point_inside_polygon(df[i, :dropoff_latitude], df[i, :dropoff_longitude], city)
	end
	df = df[mask, :];

	toc()

	tic()

	MANHATTAN_CENTER = LLA(40.782, -73.9706);

	px = Array(Float32,nrow(df));
	py = Array(Float32,nrow(df));
	dx = Array(Float32,nrow(df));
	dy = Array(Float32,nrow(df));

	for i = 1:nrow(df)
	    pENU = ENU(LLA(df[i,:pickup_latitude],df[i,:pickup_longitude]),MANHATTAN_CENTER)
	    dENU = ENU(LLA(df[i,:dropoff_latitude],df[i,:dropoff_longitude]),MANHATTAN_CENTER)
	    px[i] = pENU.east
	    py[i] = pENU.north
	    dx[i] = dENU.east
	    dy[i] = dENU.north
	end

	df[:pickup_x] = px;
	df[:pickup_y] = py;
	df[:dropoff_x] = dx;
	df[:dropoff_y] = dy;

	delete!(df,:pickup_latitude);
	delete!(df,:pickup_longitude);
	delete!(df,:dropoff_latitude);
	delete!(df,:dropoff_longitude);

	toc()

	@time writetable("reduced_trip_data_$(num).csv",df)
	println("complete reduced csv")
end

reduce(num)

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
struct UTCartesianIndices{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{CartesianIndex{N},N}
        indices::R
end


text:
# =========================================================================== #
# Version mise à jour du DM1 afin de pouvoir faire le DM2
#
# =========================================================================== #

#asked question:
#Mettre en place une heuristique de recherche locale (descente ou plus profonde descente) fondée sur deux voisinages (exemple : type “k-p exchange”).

#taking an array of lenght m
#returning another array of length m
include("DM1_1.jl")

# Notre structure de Solution
mutable struct Solution
    #""""état actuel de notre Solution"""
    x::Array{Bool,1}
    #"""fonction objectif de notre Solution"""
    objectif::Integer

    Solution(x,objectif) = new(x,objectif)
end

"""
    Calcule la fonction objectif d'une Solution
"""
function calculz(x,costs,m)
    objectif = 0
    for i in 1:m
        #checkbounds(x,m)
        objectif += x[i]*costs[i]
    end

    return objectif
end

"""
    fonction qui renvoies "vrai" si la variable j peu être entrée dans la Solution
"""

function select(crts,matrix,m,j)
    selec = true # juste une initialisation
    for l in 1:m
        selec = selec && (( crts[l] * matrix[l,j]) == 0)
    end
    return selec
end

"""
    fonction de kp-exchange avec k=1 et p=1
    fonction lourde : à optimiser
"""
function exchange1_1(solution,n,m,couts,crts,matrix)
    x_best = Solution(solution,calculz(solution,couts,n))
    # notre boucle k
    for i in 1:n
        if  x_best.x[i]  == 1
            x_prime = deepcopy(x_best)
            x_prime.x[i] = 0
            x_prime.objectif = calculz(x_prime.x,couts,n)
            # notre boucle p
            for j = 1:n
                if x_prime.x[j] == 0 && j != i && select(crts,matrix,m,j)# j != i car ça sers à rien de remettre une variable enlevée
                    x_seconde = deepcopy(x_prime)
                    x_seconde.x[j] = 1
                    x_seconde.objectif = calculz(x_seconde.x,couts,n)

                    #on met à jour notre meilleur Solution
                    if x_seconde.objectif > x_best.objectif
                        x_best = deepcopy(x_seconde)
                    end
                end
            end
        end
    end
    #x,objectif= Solution.x
    println(x_best.objectif, "\n")
    return(x_best.x,x_best.objectif)
end

"""
    fonction de kp-exchange avec k=1 et p=1
    fonction lourde : à optimiser
"""
function exchange1_2(solution,n,m,couts,crts,matrix)
    x_best = Solution(solution,calculz(solution,couts,n))
    # notre boucle k
    for i in 1:n
        if  x_best.x[i]  == 1
            x_prime = deepcopy(x_best)
            x_prime.x[i] = 0
            x_prime.objectif = calculz(x_prime.x,couts,n)
            # notre boucle p
            for j = 1:n
                if x_prime.x[j] == 0 && j != i select(crts,matrix,m,j) # j != i car ça sers à rien de remettre une variable enlevée
                    x_seconde = deepcopy(x_prime)
                    x_seconde.x[j] = 1
                    x_seconde.objectif = calculz(x_seconde.x,couts,n)
                    for p in j:n
                        if  x_seconde.x[p] == 0 && p != i && select(crts,matrix,m,p)
                            x_tierce = deepcopy(x_seconde)
                            x_tierce.x[p] = 1
                            x_tierce.objectif = calculz(x_tierce.x,couts,n)

                            if x_tierce.objectif > x_best.objectif
                                x_best = deepcopy(x_tierce)
                            end
                        end
                    end
                    #on met à jour notre meilleur Solution

                end
            end
        end
    end
    #x,objectif= Solution.x
    println(x_best.objectif, "\n")
    return(x_best.x,x_best.objectif)
end

#=
"""
    fonction de kp-exchange avec k=2 et p=1
"""
function exchange2_1(solution,n,m,couts,xbest,crts,matrix)
    x_best = Solution(solution,calculz(solution.x,couts,n))
    #premier k
    for i in 1:n
        if x_best.x[i]  == 1
            x_prime = deepcopy(x_best)
            x_prime.x[i] = 0
            x_prime.objectif = calculz(x_prime.x,couts,n)
            #deuxième k
            for j in 1:n
                if x_seconde[j] == 1
                    x_seconde = deepcopy(x_prime)
                    x_seconde.x[j] = 0
                    x_seconde.objectif = calculz(x_prime.x,couts,n)
                #notre p
                    for p in 1:n
                        if x_tierce[p] && p != i && p != j
                            x_tierce = deepcopy(x_seconde)
                            x_tierce.x[j] = 1
                            x_tierce.objectif = calculz(x_tierce.x,couts,n)
                            #on met a jour notre meilleur Solution
                            if x_tierce.objectif > x_best.objectif
                                x_best = deepcopy(x_tierce)
                            end
                        end
                    end
                end
            end
        end
    end
    #return(Solution.x,Solution.objectif)
end
=#


text:
function eki(y::Union{Array{R,1}, Array{Array{R,1},1}},
             σ::R,
             η::R,
             J::Integer,
             N::Integer,
             prior::Function,
             gmap::Function;
             ρ::R = convert(R,0.5),
             ζ::R = convert(R,2.0),
             γ0::R = convert(R,1e0),
             parallel::Bool = false,
             verbosity::Int=0,
             rerandomize::Bool = false,
             rerandom_fun::Function = prior,
             rerandom_coeff::R = convert(R,0.25),
             batched::Bool = false,
             batches::Int = 1) where R<:Real
    if batched
        @assert typeof(y) == Array{Array{R,1}, 1} "Must supply a list of data subsets for batched EKI, current type is $(typeof(y))"
        eki_batched(y, σ, η, J, N, prior, gmap, ρ, ζ, γ0, parallel, verbosity, rerandomize, rerandom_fun, rerandom_coeff, batches)
    else
        @assert typeof(y) == Array{R,1} "Must supply array of flat data for unbatched EKI, current type is $(typeof(y))"
        eki_nobatch(y, σ, η, J, N, prior, gmap, ρ, ζ, γ0, parallel, verbosity, rerandomize, rerandom_fun, rerandom_coeff)
    end
end

function ensemble_update!(u::Array{Array{R,1},1},
                          w::Array{Array{R,1},1},
                          wm::Array{R,1},
                          y::Array{R,1},
                          γ0::R,
                          σ::R,
                          ρ::R,
                          T::UniformScaling{R},
                          verbosity::Int) where R<:Real
        #Analysis
        Cuw = CrossCovarianceOperator(hcat(u...)', hcat(w...)')
        Cwwf = pheigfact(CovarianceOperator(hcat(w...)')) #this is a low-rank approx to Cww
        #set regularization
        γ = setγ(γ0, ρ, T, Cwwf, y, wm)
        if verbosity >= 2
            println("γ = $γ")
        end
        for j = 1:size(u)[1]
            yj = y.+σ.*randn(R, length(y))
            u[j] = u[j].+Cuw*(wbinv(T/γ,
            Cwwf[:vectors],
            Diagonal(convert(R,1.0)./Cwwf[:values]),
            (yj-w[j])))
        end
end

function eki_nobatch(y::Array{R,1},
                    σ::R,
                    η::R,
                    J::Integer,
                    N::Integer,
                    prior::Function,
                    gmap::Function,
                    ρ::R,
                    ζ::R,
                    γ0::R,
                    parallel::Bool,
                    verbosity::Int,
                    rerandomize::Bool,
                    rerandom_fun::Function,
                    rerandom_coeff::R) where {R<:Real}
    #Initialization
    Γ = σ^2 * I #this should use the efficient uniform scaling operator - can use both left & right division also
    T = σ^(-2) * I # Precision Matrix
    if parallel
        u = pmap(x->prior(), 1:J)
    else
        u = map(x->prior(), 1:J)
    end
    # u = [prior() for j = 1:J]
    #Main Optimization loop
    if verbosity > 0
        println("Starting up to $N iterations with $J ensemble members")
    end
    for i = 1:N
        if verbosity > 0
            println("Iteration # $i. Starting Forward Map")
        end
        #Prediction
        if parallel
            w = pmap(gmap, u)
        else
            w = map(gmap, u)
        end
        #Gum = gmap(mean(u))
        wm = mean(w)
        #Discrepancy principle
        convg = sqrt((y-wm)'*T*(y-wm))
        if verbosity > 0
            println("Iteration # $i. Discrepancy Check; Weighted Norm: $convg, Noise level: $(ζ*η)")
        end
        if convg <= ζ*η
            return mean(u)
        end
        ensemble_update!(u,w,wm,y,γ0,σ,ρ,T,verbosity)
        if rerandomize
            um = mean(u)
            if parallel
                u = pmap(x-> um .+ rerandom_coeff .* rerandom_fun(), 1:J)
            else
                u = map(x-> um .+ rerandom_coeff .* rerandom_fun(), 1:J)
            end
        end
    end
    mean(u)
end

function eki_batched(y::Array{Array{R,1},1},
                    σ::R,
                    η::R,
                    J::Integer,
                    N::Integer,
                    prior::Function,
                    gmap::Function,
                    ρ::R,
                    ζ::R,
                    γ0::R,
                    parallel::Bool,
                    verbosity::Int,
                    rerandomize::Bool,
                    rerandom_fun::Function,
                    rerandom_coeff::R,
                    batches::Int) where {R<:Real}
    #Initialization
    Γ = σ^2 * I #this should use the efficient uniform scaling operator - can use both left & right division also
    T = σ^(-2) * I # Precision Matrix
    if parallel
        u = pmap(x->prior(), 1:J)
    else
        u = map(x->prior(), 1:J)
    end
    #u = [prior() for j = 1:J]
    ya = vcat(y...)
    #Main Optimization loop
    if verbosity > 0
        println("Starting up to $N iterations with $J ensemble members")
    end
    for i = 1:N
        if verbosity > 0
            println("Iteration # $i. Starting Forward Map")
        end
        #Prediction
        if parallel
            w = pmap(x->vcat([gmap(x,b) for b in 1:length(y)]...), u)
        else
            w = map(x->vcat([gmap(x,b) for b in 1:length(y)]...), u)
        end
        #Discrepancy principle
        wm = mean(w)
        convg = sqrt((ya-wm)'*T*(ya-wm))
        if verbosity > 0
            println("Iteration # $i. Discrepancy Check; Weighted Norm: $convg, Noise level: $(ζ*η)")
        end
        if convg <= ζ*η
            return mean(u)
        end
        #set up batch ordering
        parts = kfoldperm(length(y), batches) #randomize batch groupings each iterations
        for (i, p) in enumerate(parts)
            if parallel
                w = pmap(x->vcat([gmap(x,b) for b in p]...), u)
            else
                w = map(x->vcat([gmap(x,b) for b in p]...), u)
            end
            wm = mean(w)
            ensemble_update!(u,w,wm,vcat(y[p]...),γ0,σ,ρ,T,verbosity)
            if rerandomize
                um = mean(u)
                if parallel
                    u = pmap(x-> um .+ rerandom_coeff .* rerandom_fun(), 1:J)
                else
                    u = map(x-> um .+ rerandom_coeff .* rerandom_fun(), 1:J)
                end
            end
        end
    end
    mean(u)
end


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

using Base.Test

@testset "Some other tests" begin
	@test 40+2 == 42
	@test 2*5 == 10
end


text:
export iauEpj
"""
**  Julian Date to Julian Epoch.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  support function.
**
**  Given:
**     dj1,dj2    double     Julian Date (see note)
**
**  Returned (function value):
**                double     Julian Epoch
**
**  Note:
**
**     The Julian Date is supplied in two pieces, in the usual SOFA
**     manner, which is designed to preserve time resolution.  The
**     Julian Date is available as a single number by adding dj1 and
**     dj2.  The maximum resolution is achieved if dj1 is 2451545.0
**     (J2000.0).
**
**  Reference:
**
**     Lieske, J.H., 1979, Astron.Astrophys. 73, 282.
**
**  This revision:  2013 August 7
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauEpj(double dj1, double dj2)
function iauEpj(dj1::Real, dj2::Real)
   return ccall((:iauEpj, libsofa_c), Cdouble,
               (Cdouble, Cdouble),
               convert(Float64, dj1), convert(Float64, dj2))
end


text:
import LinearAlgebra.tr
export tr

function tr(e::AbstractExpr)
    return sum(diag(e))
end


text:
module JETTest

let
    README = normpath(dirname(@__DIR__), "README.md")
    include_dependency(README)
    s = """$(read(README, String))

    The JETTest.jl toolset includes:
    - [Dispatch Analysis](@ref): automatically detects possible performance pitfalls, where the optimization was failed and/or runtime dispatch will happen
    """
    @doc s JETTest
end

const CC = Core.Compiler

# imports
# =======

using JET.JETInterfaces

import JET:
    get_cache_key

import Test:
    record

# usings
# ======

using JET

using JET:
    State,
    @invoke,
    @isexpr,
    gen_call_with_extracted_types_and_kwargs,
    get_reports,
    print_reports

using .CC:
    AbstractInterpreter,
    InferenceState,
    InferenceResult,
    OptimizationState,
    OptimizationParams,
    IRCode,
    find_throw_blocks,
    widenconst,
    argextype

using Core:
    Builtin,
    Const,
    CodeInfo,
    MethodInstance

using Base:
    isconcretedispatch

using Test:
    Test,
    Pass, Fail, Broken, Error,
    Threw,
    get_testset,
    TESTSET_PRINT_ENABLE,
    AbstractTestSet,
    DefaultTestSet,
    FallbackTestSetException

# filters
# =======

function module_filter(m)
    return function (frame::State)
        frame.mod === m
    end
end

# includes
# ========

include("dispatch.jl")

export
    analyze_dispatch,
    @analyze_dispatch,
    report_dispatch,
    @report_dispatch,
    @test_nodispatch,
    module_filter

end


text:
"""
    greenseltype(::Type{DQMC}, m::Model)

Returns the type of the elements of the Green's function matrix. Defaults to `ComplexF64`.
"""
greenseltype(::Type{DQMC}, m::Model) = ComplexF64

"""
    hoppingeltype(::Type{DQMC}, m::Model)

Returns the type of the elements of the hopping matrix. Defaults to `Float64`.
"""
hoppingeltype(::Type{DQMC}, m::Model) = Float64

"""
    energy(mc::DQMC, m::Model, conf)

Calculate bosonic part (non-Green's function determinant part) of energy for configuration `conf` for Model `m`.
"""
energy_boson(mc::DQMC, m::Model, conf) = 0.


text:
function ueberpruefePixelStrahl(PixelNummer,Strahl, Pixel1,Pixel2, Pixel3, Pixel4, Pixel5)

if PixelNummer[1,Strahl] == Pixel1
    if PixelNummer[2,Strahl] == Pixel2
        if PixelNummer[3,Strahl] == Pixel3
            if PixelNummer[4,Strahl] == Pixel4
                if PixelNummer[5,Strahl] == Pixel5
                    println("Ihr habt korrekt erkannt durch welche Pixel der Strahl läuft!")
                else
                    println("Überprüft noch einmal den Eintrag für Pixel5.")
                end
            else
                println("Überprüft noch einmal den Eintrag für Pixel4.")
            end
        else
            println("Überprüft noch einmal den Eintrag für Pixel3.")
        end
    else
        println("Überprüft noch einmal den Eintrag für Pixel2.")
    end
else
    println("Überprüft noch einmal den Eintrag für Pixel1.")
end

end


text:
#
# Creates phylogenetic tree of organisms data sets. These data sets will
# be stored in JSON files on HDD. One data set is stored after application
# reloaded. This approach will create small and separate JSON files for
# later visualization. Format of JSON:
#
# {
#   nodes: [{id: Number, org: String, mutations: [{org: String, amount: Number},...]},...],
#   edges: [{from:Number, to:Number},...]
# }
#
# where:
#   - id        - Organism unique id
#   - org       - Organism object in string representation
#   - mutations - All mutations of organism within his life
#     - org     - Same like org above
#     - amount  - Amount of mutations per one mutation session
#   - from      - Relation source organism id
#   - to        - Relation destination organism id
#
# @author DeadbraiN
#
module Phylogen
  import Helper
  import Creature
  import Event
  import Backup
  import ManagerTypes.ManagerData
  import ManagerTypes.Plugin

  export init
  #
  # Default postfix of phylogenetic tree files
  #
  const PHYLO_FILE_POSTFIX = "-jevo-phylo.json"
  #
  # Default name of the folder for phylogenetic tree JSON files
  #
  const PHYLO_FOLDER_NAME = "phylogen"
  #
  # Name of the current module
  #
  const MODULE_NAME = string(Phylogen)
  #
  # Describes one mutation of organism. Mutated organism - it's
  # original (clonned) organism + mutations.
  #
  type PhylogenMutation
    #
    # Mutated organism
    #
    org::Creature.Organism
    #
    # Amount of real mutations applied to original organism
    #
    mutations::Int
  end
  #
  # Describes all mutations of organism. Organism and mutations
  # should be separate copies from original organisms.
  #
  type PhylogenOrganism
    #
    # Organism object after birth
    #
    org::Creature.Organism
    #
    # Array of organisms based on org, but with mutations. This
    # is something like it's future versions.
    #
    mutations::Array{PhylogenMutation, 1}
  end
  #
  # Contains Phylogenetic tree of organisms. These data will
  # be used for Phylogenetic tree visualization.
  #
  type PhylogenData <: Plugin
    #
    # Map of organism copies. It's important to have full copy without
    # meta information about their code (funcs field should be empty).
    # Key is organism's unique id.
    #
    organisms::Dict{UInt, PhylogenOrganism}
    #
    # Array of relations. Relation it's two organism ids, that links
    # them together. For example parent and child organism ids. These
    # data will be used for phylogenetic graph creation.
    #
    relations::Array{UInt, 1}
    #
    # Set of killed organism ids. By killed, i mean organisms, who
    # were killed before they were saved to JSON file. We have to
    # store them in killedOrganisms, because JSON will have releation
    # to undefined (killed) organism and graph wisualizer will fail
    #
    killedIds::Set{UInt}
    #
    # Organisms, who were killed before they will be stored in JSON file.
    # See killedIds for details
    #
    killedOrganisms::Dict{UInt, PhylogenOrganism}
  end
  #
  # Module initializer
  #
  function init(man::ManagerData)
    #
    # We havr to add ourself to plugins map
    #
    man.plugins[MODULE_NAME] = PhylogenData(Dict{UInt, PhylogenOrganism}(), [], Set{UInt}(), Dict{UInt, PhylogenOrganism}())
    #
    # All event handlers should be binded here
    #
    Event.on(man.obs, "backup", _save)
    Event.on(man.obs, "backup", _clear)
    Event.on(man.obs, "killorganism", _delOrganism)
    Event.on(man.obs, "bornorganism", _addOrganism)
    Event.on(man.obs, "mutations", _addMutations)
    Event.on(man.obs, "clone", _addRelation)

  end
  #
  # Adds one organism to phylogenetic organisms pool. Creates it's
  # copy without meta information. Mutations array will be empty.
  # We have to call this method every time when new organism in a
  # world will be created. One organism in phylogenetic tree is a
  # node with relations. Every related node is a child or parent
  # organism.
  # @param man Manager related data object
  # @param org Organism we have to add to the pool
  #
  function _addOrganism(man::ManagerData, org::Creature.Organism)
    local orgCopy::Creature.Organism = _copyOrganism(man, org)

    if haskey(man.plugins[MODULE_NAME].organisms, org.id)
      Helper.error(string("Phylogen: Dublicate organism id: ", org.id))
    end
    man.plugins[MODULE_NAME].organisms[org.id] = PhylogenOrganism(orgCopy, [])
  end
  #
  # Opposite to addOrganism() function. Removes one organism
  # from phylogenetic tree pool. Should be called if one organism
  # has died in a world.
  # @param man Manager related data object
  # @param org Organism we have to add to the pool
  #
  function _delOrganism(man::ManagerData, org::Creature.Organism)
    #
    # We have unsaved relation for this organism, so we have to
    # store it in separate dictionary
    #
    if in(org.id, man.plugins[MODULE_NAME].killedIds)
      man.plugins[MODULE_NAME].killedOrganisms[org.id] = man.plugins[MODULE_NAME].organisms[org.id]
    end
    delete!(man.plugins[MODULE_NAME].organisms, org.id)
  end
  #
  # Adds specified amount of mutations to organism. org parameter
  # should be passed after all mutations are applied to him. Mutations
  # are not a nodes in a tree. They are additional information of
  # node (organism).
  # @param man Manager related data object
  # @param org Organism we have to add to the pool
  # @param amount Amount of real mutations
  # @param onClone true if current mutations were applied on clone
  #
  function _addMutations(man::ManagerData, org::Creature.Organism, amount::Int, onClone::Bool)
    local orgs::Dict{UInt, PhylogenOrganism} = man.plugins[MODULE_NAME].organisms
    if !haskey(orgs, org.id)
      Helper.error(string("Phylogen: Unknown organism id: ", org.id))
      return nothing
    end
    local orgCopy::Creature.Organism = _copyOrganism(man, org)

    push!(orgs[org.id].mutations, PhylogenMutation(orgCopy, amount))
  end
  #
  # Adds relation (link) between parent and child organisms. Should be
  # called every time, when clonning or crossover is called.
  # @param man Manager related data object
  # @param parentId Unique id of parent organism
  # @param childId Unique id of child organism
  #
  function _addRelation(man::ManagerData, parentId::UInt, childId::UInt)
    local data::PhylogenData = man.plugins[MODULE_NAME]

    push!(data.relations, parentId)
    push!(data.relations, childId)
    push!(data.killedIds, parentId)
    push!(data.killedIds, childId)
  end
  #
  # Clears phylogenetic tree arrays to start new portion of data. Only
  # relations may be cleared, because organisms pool should be actual
  # all the time.
  # @param man Manager related data object
  #
  function _clear(man::ManagerData)
    man.plugins[MODULE_NAME].relations = [];
    man.plugins[MODULE_NAME].killedIds = Set{UInt}();
    man.plugins[MODULE_NAME].killedOrganisms = Dict{UInt, PhylogenOrganism}()
  end
  #
  # Saves JSON data to the file on HDD.
  # @param man Manager related data object
  #
  function _save(man::ManagerData)
    Backup.save(_toJson(man), PHYLO_FOLDER_NAME, PHYLO_FILE_POSTFIX, true)
  end
  #
  # Converts phylogenetic tree arrays to JSON.
  # @param man Manager related data object
  #
  function _toJson(man::ManagerData)
    local org::PhylogenOrganism
    local i::Int
    local id::UInt
    local json::String
    local relations::Array{UInt, 1} = man.plugins[MODULE_NAME].relations
    local nodes::Array{String, 1} = []
    local edges::Array{String, 1} = []

    json = "{\"nodes\":["
    _fillBy(man.plugins[MODULE_NAME].organisms, nodes)
    _fillBy(man.plugins[MODULE_NAME].killedOrganisms, nodes)
    json *= join(nodes, ",")
    json *= "]"

    for i = 1:2:length(relations) push!(edges, string("{\"source\":\"", relations[i],"\",\"target\":\"", relations[i+1], "\"}")) end
    if length(edges) < 1 return json end
    json *= ",\"edges\":["
    json *= join(edges, ",")
    json *= "]}"

    json
  end
  #
  # Fills "nodes" array by "organisms" data
  # @param organisms
  # @param nodes
  # TODO: add mutations data here
  function _fillBy(orgs::Dict{UInt, PhylogenOrganism}, nodes::Array{String, 1})
    local org::PhylogenOrganism
    local id::UInt

    for (id, org) in orgs
      push!(nodes, string(
        "{\"id\": \"", id,
        "\",\"code\":\"", replace(replace(string(org.org.code), "\n", "\\n"), "\"", "\\\""),
        "\"}"))
    end
  end
  #
  # Makes full organism copy, but removes all meta information.
  # @param man Manager data object
  # @param org Organism we have to copy
  # @return {Organism}
  #
  function _copyOrganism(man::ManagerData, org::Creature.Organism)
    local orgCopy::Creature.Organism = Creature.create(org, man.cfg, org.id, org.pos)

    orgCopy.funcs  = []
    orgCopy.codeFn = Helper.emptyFn
    Event.clear(orgCopy.observer)

    orgCopy
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# # Export PubMed Citations

#md # [![nbviewer](https://img.shields.io/badge/jupyter_notebook-nbviewer-orange.svg)](http://nbviewer.jupyter.org/github/bcbi/BioMedQuery.jl/tree/master/docs/src/notebooks/4_pubmed_export_citations.ipynb)

# BioMedQuery has functions to search PubMed for PMIDs and save the xml data
# as either a BibTex or EndNote citation.
#
# Here we export EndNote/BibTex citations from a PMID or a list of PMIDs. If you need
# to search Entrez/PubMed and save the results as citations, refer to Examples / PubMed Search and Save.

# ### Set Up

using BioMedQuery.Processes

# The functions require a directory to save the citation files to
results_dir = ".";

if !isdir(results_dir)
     mkdir(results_dir)
end

# For this example, the below PMIDs are searched and saved as citations
pmid = 11748933;
pmid_list = [24008025, 24170597];

# ### Export as an EndNote library file

# Saving one PMID's citaiton as an EndNote file
enw_file = results_dir * "/11748933.enw";
export_citation(pmid, "endnote", enw_file);
sleep(1) # hide

# Saving two PMIDs' citations as an EndNote file
enw_file = results_dir * "/pmid_list.enw";
export_citation(pmid_list, "endnote", enw_file);
sleep(1) # hide

# #### Explore one of the output files

println(read(enw_file, String))

# ### Export as a Bibtex file

# Saving one PMID's citation as a BibTex file
bib_file = results_dir * "/11748933.bib";
export_citation(pmid, "bibtex", bib_file);
sleep(1) # hide

# Saving two PMIDs' citations as a BibTex file
bib_file = results_dir * "/pmid_list.bib";
export_citation(pmid_list, "bibtex", bib_file);
sleep(1) # hide

# #### Explore one of the output files

println(read(bib_file, String))


text:
using NNlibCPU
using Documenter

DocMeta.setdocmeta!(NNlibCPU, :DocTestSetup, :(using NNlibCPU); recursive=true)

makedocs(;
    modules=[NNlibCPU],
    authors="chriselrod <elrodc@gmail.com> and contributors",
    repo="https://github.com/chriselrod/NNlibCPU.jl/blob/{commit}{path}#{line}",
    sitename="NNlibCPU.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://chriselrod.github.io/NNlibCPU.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/chriselrod/NNlibCPU.jl",
)


text:
module TestSmoke

using Aqua
Aqua.test_all(Aqua)

end  # module