text:
module TestAudioIO

using FactCheck
using AudioIO
import AudioIO.AudioBuf

const TEST_SAMPLERATE = 44100
const TEST_BUF_SIZE = 1024

include("testhelpers.jl")

type TestAudioStream <: AudioIO.AudioStream
    root::AudioIO.AudioMixer
    info::AudioIO.DeviceInfo

    function TestAudioStream()
        root = AudioMixer()
        new(root, AudioIO.DeviceInfo(TEST_SAMPLERATE, TEST_BUF_SIZE))
    end
end

# render the stream and return the next block of audio. This is used in testing
# to simulate the audio callback that's normally called by the device.
function process(stream::TestAudioStream)
    out_array = zeros(AudioIO.AudioSample, stream.info.buf_size)
    in_array = zeros(AudioIO.AudioSample, stream.info.buf_size)
    rendered = AudioIO.render(stream.root, in_array, stream.info)
    out_array[1:length(rendered)] = rendered
    return out_array
end

#### Test playing back various vector types ####

facts("Array playback") do
    # data shared between tests, for convenience
    t = linspace(0, 2, 2 * 44100)
    phase = 2pi * 100 * t

    ## Test Float32 arrays, this is currently the native audio playback format
    context("Playing Float32 arrays") do
        f32 = convert(Array{Float32}, sin(phase))
        test_stream = TestAudioStream()
        player = play(f32, test_stream)
        @fact process(test_stream) => f32[1:TEST_BUF_SIZE]
    end

    context("Playing Float64 arrays") do
        f64 = convert(Array{Float64}, sin(phase))
        test_stream = TestAudioStream()
        player = play(f64, test_stream)
        @fact process(test_stream) => convert(AudioBuf, f64[1:TEST_BUF_SIZE])
    end

    context("Playing Int8(Signed) arrays") do
        i8 = Int8[-127:127]
        test_stream = TestAudioStream()
        player = play(i8, test_stream)
        @fact process(test_stream)[1:255] =>
                mse(convert(AudioBuf, linspace(-1.0, 1.0, 255)))
    end

    context("Playing Uint8(Unsigned) arrays") do
        # for unsigned 8-bit audio silence is represented as 128, so the symmetric range
        # is 1-255
        ui8 = Uint8[1:255]
        test_stream = TestAudioStream()
        player = play(ui8, test_stream)
        @fact process(test_stream)[1:255] =>
                mse(convert(AudioBuf, linspace(-1.0, 1.0, 255)))
   end
end

facts("AudioNode Stopping") do
    test_stream = TestAudioStream()
    node = SinOsc(440)
    play(node, test_stream)
    process(test_stream)
    stop(node)
    @fact process(test_stream) => zeros(AudioIO.AudioSample, TEST_BUF_SIZE)
end

facts("Audio Device Listing") do
    # there aren't any devices on the Travis machine so just test that this doesn't crash
    @fact get_audio_devices() => issubtype(Array)
end

end # module TestAudioIO


text:
ENV["JULIA_DEBUG"] = "all"
using HierarchicalTemporalMemory, Test

# Test iteration over a hypercube
startP= (50,80)
radius= 30
hc1= Hypercube(startP,radius,(200,200));
a= [convert(NTuple{2,Int},i) for i in hc1]
reference= [i for i in zip(
    repeat(startP[1]-radius : startP[1]+radius,outer=2*radius+1),
    repeat(startP[2]-radius : startP[2]+radius,inner=2*radius+1))]
@test a == reference

ext: jl


text:
using PermutationGroups
using Groups.KnuthBendix

@testset "Wajnryb presentation for Σ₄₁" begin

    genus = 4

    Fn = FreeGroup(2genus)
    G = SpecialAutomorphismGroup(Fn)

    T = Groups.mcg_twists(G)

    # symplectic pairing in the free Group goes like this:
    # f1 ↔ f5
    # f2 ↔ f6
    # f3 ↔ f7
    # f4 ↔ f8

    T = let G = G
        (Tas, Tαs, Tes) = Groups.mcg_twists(G)
        Ta = G.(Tas)
        Tα = G.(Tαs)
        Tes = G.(Tes)

        [Ta; Tα; Tes]
    end

    a1 = T[1]^-1 # Ta₁
    a2 = T[5]^-1 # Tα₁
    a3 = T[9]^-1 # Te₁₂
    a4 = T[6]^-1 # Tα₂
    a5 = T[12]^-1 # Te₂₃
    a6 = T[7]^-1 # Tα₃
    a7 = T[14]^-1 # Te₃₄
    a8 = T[8]^-1 # Tα₄

    b0 = T[2]^-1 # Ta₂
    a0 = (a1 * a2 * a3)^4 * b0^-1 # from the 3-chain relation
    X = a4 * a5 * a3 * a4 # auxillary, not present in the Primer
    b1 = X^-1 * a0 * X
    b2 = T[10]^-1 # Te₁₃

    As = T[[1, 5, 9, 6, 12, 7, 14, 8]] # the inverses of the elements a

    @testset "preserving relator" begin
        F = Groups.object(G)

        R = prod(commutator(gens(F,2i+1), gens(F,2i+2)) for i in 0:genus-1)

        for g in T
            @test g(R) == R
        end
    end

    @testset "commutation relations" begin
        for (i, ai) in enumerate(As) #the element ai here is actually the inverse of ai before. It does not matter for commutativity. Also, a0 is as defined before.
            for (j, aj) in enumerate(As)
                if abs(i - j) > 1
                    @test ai * aj == aj * ai
                elseif i ≠ j
                    @test ai * aj != aj * ai
                end
            end
            if i != 4
                @test a0 * ai == ai * a0
            end
        end
    end

    @testset "braid relations" begin
        for (i, ai) in enumerate(As) #the element ai here is actually the inverse of ai before. It does not matter for braid relations.
            for (j, aj) in enumerate(As)
                if abs(i - j) == 1
                    @test ai * aj * ai == aj * ai * aj
                end
            end
        end
        @test a0 * a4 * a0 == a4 * a0 * a4 # here, a0 and a4 are as before
    end

    @testset "3-chain relation" begin
        x = a4*a3*a2*a1*a1*a2*a3*a4 # auxillary; does not have a name in the Primer
        @test b0 == x*a0*x^-1
    end

    @testset "Lantern relation" begin

        @testset "b2 definition" begin
            @test b2 == (a2 * a3 * a1 * a2)^-1 * b1 * (a2 * a3 * a1 * a2)

            # some additional tests, checking what explicitly happens to the generators of the π₁ under b2
            d = Groups.domain(b2)
            img = evaluate(b2)
            z = img[3] * d[3]^-1

            @test img[1] == d[1]
            @test img[2] == d[2]
            @test img[3] == z * d[3]
            @test img[4] == z * d[4] * z^-1
            @test img[5] == z * d[5] * z^-1
            @test img[6] == z * d[6] * z^-1
            @test img[7] == d[7] * z^-1
            @test img[8] == d[8]
        end

        @testset "b2: commutation relations" begin
            @test b2 * a1 == a1 * b2
            @test b2 * a2 != a2 * b2
            @test b2 * a3 == a3 * b2
            @test b2 * a4 == a4 * b2
            @test b2 * a5 == a5 * b2
            @test b2 * a6 != a6 * b2
        end

        @testset "b2: braid relations" begin
            @test a2 * b2 * a2 == b2 * a2 * b2
            @test a6 * b2 * a6 == b2 * a6 * b2
        end

        @testset "lantern" begin
            u = (a6 * a5)^-1 * b1 * (a6 * a5)
            x = (a6 * a5 * a4 * a3 * a2 * u * a1^-1 * a2^-1 * a3^-1 * a4^-1) # yet another auxillary
            # x = (a4^-1*a3^-1*a2^-1*a1^-1*u*a2*a3*a4*a5*a6)
            @time evaluate(x)
            b3 = x * a0 * x^-1
            b3im = @time evaluate(b3)
            b3cim = @time let g = b3
                f = Groups.compiled(g)
                f(Groups.domain(g))
            end
            @test b3im == b3cim
            @test a0 * b2 * b1 == a1 * a3 * a5 * b3
        end
    end

    Base.conj(t::Groups.Transvection, p::Perm) =
        Groups.Transvection(t.id, t.i^p, t.j^p, t.inv)

    function Base.conj(elt::FPGroupElement, p::Perm)
        G = parent(elt)
        A = alphabet(elt)
        return G([A[conj(A[idx], p)] for idx in word(elt)])
    end

    @testset "Te₂₃ definition" begin
        Te₁₂, Te₂₃ = T[9], T[12]
        G = parent(Te₁₂)
        F₈ = Groups.object(G)
        (δ, d, γ, c, β, b, α, a) = Groups.gens(F₈)

        Groups.domain(Te₁₂)

        img_Te₂₃ = evaluate(Te₂₃)
        y = c * β^-1 * b^-1 * β
        @test img_Te₂₃ == (δ, d, y * γ, y * c * y^-1, β * y^-1, b, α, a)

        σ = perm"(7,5,3)(8,6,4)"
        Te₂₃_σ = conj(Te₁₂, σ)
        # @test word(Te₂₃_σ) == word(Te₂₃)

        @test evaluate(Te₂₃_σ) == evaluate(Te₂₃)
        @test Te₂₃ == Te₂₃_σ
    end

    @testset "Te₃₄ definition" begin
        Te₁₂, Te₃₄ = T[9], T[14]
        G = parent(Te₁₂)
        F₈ = Groups.object(G)
        (δ, d, γ, c, β, b, α, a) = Groups.gens(F₈)

        σ = perm"(7,3)(8,4)(5,1)(6,2)"
        Te₃₄_σ = conj(Te₁₂, σ)
        @test Te₃₄ == Te₃₄_σ
    end

    @testset "hyperelliptic τ is central" begin

        τ = Groups.rotation_element(G)
        τᵍ = τ^genus

        symplectic_gens = let genus = genus, G = G
            π₁Σ = Groups.SurfaceGroup(genus, 0)
            s_twists = Groups.symplectic_twists(π₁Σ)
            G.(word(t.autFn_word) for t in s_twists)
        end

        @test all(sg * τᵍ == τᵍ * sg for sg in symplectic_gens)
    end
end


text:
module PylonCameras

using Cameras
using ColorTypes
using FixedPointNumbers
using ImageCore
using ResourcePools

import Cameras:
    isopen,
    open!,
    close!,
    isrunning,
    start!,
    stop!,
    take!,
    trigger!,
    id,
    timestamp

export PylonCamera,
    isopen,
    open!,
    close!,
    isrunning,
    start!,
    stop!,
    take!,
    trigger!,
    id,
    image_number,
    timestamp

include("wrapper.jl")
include("pylon_acquired_image.jl")

mutable struct PylonCamera <: Camera
    device::Wrapper.IPylonDevice
    instant_camera::Wrapper.InstantCamera
    grab_result_ready_cond::Base.AsyncCondition
    grab_result_wait_timeout_ms::UInt32
    grab_result_retrieve_timeout_ms::UInt32
    feature_filename::Union{String, Nothing}
    function PylonCamera(;
            max_num_buffer = 10,
            grab_result_wait_timeout_ms = 5000,
            grab_result_retrieve_timeout_ms = 100,
            feature_filename::Union{String, Nothing} = nothing)
        Wrapper.pylon_initialize()
        transport_layer_factory = Wrapper.get_transport_layer_factory_instance()
        grab_result_ready_cond = Base.AsyncCondition()

        device = Wrapper.create_first_device(transport_layer_factory)
        instant_camera = Wrapper.InstantCamera(device)
        Wrapper.max_num_buffer!(instant_camera, UInt(max_num_buffer))
        new(device, instant_camera, grab_result_ready_cond, grab_result_wait_timeout_ms, grab_result_retrieve_timeout_ms, feature_filename)
    end
end

function info(c::PylonCamera)
    device_info = Wrapper.get_device_info(c.device)
    vendor_name = Wrapper.get_vendor_name(device_info)
    model_name = Wrapper.get_model_name(device_info)
    serial_number = Wrapper.get_serial_number(device_info)
    return vendor_name, model_name, serial_number
end

isopen(c::PylonCamera) = Wrapper.is_open(c.instant_camera)

function open!(c::PylonCamera)
    if c.feature_filename != nothing
        @debug "Removing default configuration from instance"
        Wrapper.register_configuration(c.instant_camera, C_NULL, Wrapper.RegistrationMode_ReplaceAll, Wrapper.Cleanup_None)
    end
    @debug "Opening camera"
    Wrapper.open(c.instant_camera)
    if c.feature_filename != nothing
        @debug "Getting camera node map"
        node_map = Wrapper.get_node_map(c.instant_camera)
        @debug "Loading node map from $(c.feature_filename)"
        Wrapper.load_features(c.feature_filename, node_map)
    end
end

close!(c::PylonCamera) = Wrapper.close(c.instant_camera)

isrunning(c::PylonCamera) = Wrapper.is_grabbing(c.instant_camera)

start!(c::PylonCamera) = Wrapper.start_grabbing_async(c.instant_camera, c.grab_result_wait_timeout_ms, Wrapper.notify_async_cond_safe_c, c.grab_result_ready_cond.handle)
start!(c::PylonCamera, images_to_grab::Int) = Wrapper.start_grabbing_async(c.instant_camera, UInt64(images_to_grab), c.grab_result_wait_timeout_ms, Wrapper.notify_async_cond_safe_c, c.grab_result_ready_cond.handle)
stop!(c::PylonCamera) = Wrapper.stop_grabbing(c.instant_camera)

function take!(c::PylonCamera)::AbstractAcquiredImage
    @debug "Waiting for result"
    wait(c.grab_result_ready_cond)
    @debug "Retrieving result for $(c.grab_result_retrieve_timeout_ms) ms"
    # Be careful not to pass grab_result on (e.g. to active logging like @info),
    # as it may cause the grab_result to be held for too long exhausting the pylon camera's buffer pool
    grab_result = Wrapper.retrieve_result(c.instant_camera, c.grab_result_retrieve_timeout_ms)
    @debug "Retrieved $(grab_result)"
    if Wrapper.grab_succeeded(grab_result)
        return PylonAcquiredImage(grab_result)
    else
        error("$(Wrapper.get_error_code(grab_result)) $(Wrapper.get_error_description(grab_result))")
    end
end

trigger!(camera::PylonCamera) = error("Not implemented yet")

end # module


text:
module DeepLife
using Genie, Genie.Router, Genie.Renderer, Genie.AppServer

function main()
  Base.eval(Main, :(const UserApp = DeepLife))

  include("../genie.jl")

  Base.eval(Main, :(const Genie = DeepLife.Genie))
  Base.eval(Main, :(using Genie))
end

main()

end


text:
ENV["PYTHON"] = ""

using PyCall

py"""
import numpy as np
import tensorflow as tf
import tensorflow_constrained_optimization as tfco
import time

def batch_provider(data, num_epochs=10, shuffle=True, batch_size=32):
    def make_batch():
        dataset = tf.data.Dataset.from_tensor_slices(data)
        if shuffle:
            dataset = dataset.shuffle(1000)
        dataset = dataset.batch(batch_size).repeat(num_epochs)
        return dataset
    return make_batch

def build_model_mnist():
    return tf.keras.Sequential([
        tf.keras.layers.Conv2D(20, kernel_size=(5,5), activation='relu', input_shape=(28,28,1), strides=(1, 1)),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Conv2D(50, kernel_size=(5,5), activation='relu', strides=(1, 1)),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Flatten(),
        tf.keras.layers.Dense(500),
        tf.keras.layers.Dense(1),
    ])

def build_model_cifar():
    return tf.keras.Sequential([
        tf.keras.layers.Conv2D(64, kernel_size=3, activation='relu', input_shape=(32,32,3), padding="valid"),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Conv2D(128, kernel_size=3, activation='relu', padding="valid"),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Conv2D(128, kernel_size=3, activation='relu', padding="valid"),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Flatten(),
        tf.keras.layers.Dense(1),
    ])

def build_model_svhn2():
    return tf.keras.Sequential([
        tf.keras.layers.Conv2D(64, kernel_size=3, activation='relu', input_shape=(32,32,3), padding="valid"),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Conv2D(128, kernel_size=3, activation='relu', padding="valid"),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Conv2D(128, kernel_size=3, activation='relu', padding="valid"),
        tf.keras.layers.MaxPooling2D(pool_size=(2, 2)),
        tf.keras.layers.Flatten(),
        tf.keras.layers.Dense(1),
    ])

def build_model_imagenet():
    return tf.keras.Sequential([
        tf.keras.layers.Dense(1, input_shape=(62720,)),
    ])

def build_model_molecules():
    return tf.keras.Sequential([
        tf.keras.layers.Dense(50, input_shape=(100,), activation='sigmoid'),
        tf.keras.layers.BatchNormalization(),
        tf.keras.layers.Dense(25, activation='sigmoid'),
        tf.keras.layers.BatchNormalization(),
        tf.keras.layers.Dense(1),
    ])

def make_model_fn(build_model, fpr, optimizer, steplength):
    def model_fn(features, labels, mode):
        model = build_model()
        scores = model(features)

        # output predictions
        if mode == tf.estimator.ModeKeys.PREDICT:
            return tf.estimator.EstimatorSpec(mode, predictions=scores)

        # Baseline cross-entropy loss.
        loss_fn = tf.keras.losses.BinaryCrossentropy(from_logits=True)
        loss = loss_fn(labels, scores)
        train_op = None

        # Set up precision at recall optimization problem.
        if mode == tf.estimator.ModeKeys.TRAIN:
            # Create optimization problem
            context = tfco.rate_context(scores, labels)
            problem = tfco.RateMinimizationProblem(
                tfco.false_negative_rate(context),
                [tfco.false_positive_rate(context) <= fpr]
            )

            lag_loss, update_ops_fn, multipliers = tfco.create_lagrangian_loss(problem)

            # Set up optimizer and the list of variables to optimize the loss.
            opt = optimizer(learning_rate=steplength)
            opt.iterations = tf.compat.v1.train.get_or_create_global_step()

            # Get minimize op and group with update_ops.
            var_list = (model.trainable_weights +
                        problem.trainable_variables + [multipliers])
            minimize_op = opt.get_updates(lag_loss(), var_list)
            update_ops = update_ops_fn()
            train_op = tf.group(*update_ops, minimize_op)

        return tf.estimator.EstimatorSpec(
            mode,
            loss=loss,
            train_op=train_op
        )
    return model_fn

def train_model(build_model, x_train, y_train, x_test, y_test, args, model_dir):
    # create model
    classifier = tf.estimator.Estimator(
        model_fn = make_model_fn(build_model, args["fpr"], args["optimiser"], args["steplength"]),
        config=tf.estimator.RunConfig(tf_random_seed=args["seed"]),
        model_dir = model_dir,
    )

    # train model
    make_batch = batch_provider((x_train, y_train), num_epochs=args["epochs"], shuffle=True, batch_size=args["batchsize"])

    start_time = time.time()
    classifier.train(make_batch)
    elapsed_time = time.time() - start_time

    # compute scores
    eval_data_train = batch_provider((x_train, y_train), num_epochs=1, shuffle=False)
    eval_data_test = batch_provider((x_test, y_test), num_epochs=1, shuffle=False)

    eval_train = classifier.evaluate(eval_data_train)
    eval_test = classifier.evaluate(eval_data_test)

    s_train = np.array(list(classifier.predict(input_fn=eval_data_train)))
    s_test = np.array(list(classifier.predict(input_fn=eval_data_test)))

    return {
        "s_train": s_train,
        "loss_train": eval_train['loss'],
        "s_test": s_test,
        "loss_test": eval_test['loss'],
        "tm": elapsed_time,
        "iters": eval_train['global_step'],
    }
"""

build_network_tfco(::Type{<:AbstractMNIST}) = py"build_model_mnist"
build_network_tfco(::Type{<:AbstractCIFAR}) = py"build_model_cifar"
build_network_tfco(::Type{<:AbstractSVHN2}) = py"build_model_svhn2"
build_network_tfco(::Type{<:Molecules}) = py"build_model_molecules"
build_network_tfco(::Type{<:ImageNet}) = py"build_model_imagenet"

function reshape_for_python(x, y)
    d = ndims(x)
    return Float32.(permutedims(x, vcat(d, 1:(d-1)))), Float32.(vec(y))
end

function run_simulations_tfco(Dataset_Settings, Train_Settings, Model_Settings)
    model_dir = mktempdir("."; cleanup=false)
    for dataset_settings in dict_list_simple(Dataset_Settings)
        @unpack dataset, posclass = dataset_settings
        @info "Dataset: $(dataset), positive class label: $(posclass)"

        labelmap = (y) -> y == posclass
        train, test = load(dataset; labelmap = labelmap)
        (x_train, y_train) = reshape_for_python(train...);
        (x_test, y_test) = reshape_for_python(test...);

        for train_settings in dict_list_simple(Train_Settings)
            @unpack batchsize, epochs, seed, optimiser, steplength = train_settings
            @info "Batchsize: $(batchsize), runid: $(seed)"

            for model_settings in dict_list_simple(Model_Settings)
                @unpack type, arg, surrogate, reg, buffer = model_settings
                model_settings[:seed] = seed

                if optimiser <: Descent
                    optm = py"tf.keras.optimizers.SGD"
                else
                    @error "unknown optimiser"
                end

                # create model
                settings = Dict(
                    :fpr => arg,
                    :seed => seed,
                    :optimiser => optm,
                    :steplength => steplength,
                    :batchsize => batchsize,
                    :epochs => epochs
                )

                build_net = build_network_tfco(dataset)
                res = py"train_model"(build_net, x_train, y_train, x_test, y_test, settings, model_dir)
                s_train, s_test, elapsed_time = res

                save_simulation_tfco(
                    dataset_settings,
                    train_settings,
                    model_settings,
                    res,
                    y_train,
                    y_test,
                )
            end
        end
    end
    rm(model_dir; recursive=true)
end


text:
function solver_OPF_dual(optimizer::Type{<:AbstractOptimizer},resultfile::String="")
#function solver_DeterministicEquivalent(optimizer::Type{<:AbstractOptimizer},resultfile::String="",load_st2::Type{<:AbstractArray{Int64,2}})
#Deterministic Equivalent------------------------------------------------------#
println("Loading inputs...")

(syscost_det,Dual_constraint2,Dual_constraint1,λ_value,f_value)=OPF_dual(optimizer,set_generators,set_nodes,set_demands,generators,demands,capacity,varcost,load_node,utility_node,links,links_rev,F_max_dict,B_dict,MapG,MapD)

resultfile != "" && open(resultfile, truncate = true) do f

    println(f,"*****************OPF Solution (Dual Problem)*****************")
    println(f,"")
    println(f,"Objective Function:",syscost_det)
    println(f,"")
    println(f,"Production level of generator g:",Dual_constraint2)
    println(f,"")
    println(f,"Consumption level of demand d:",Dual_constraint1)
    println(f,"")
    println(f,"Electricity Price:",λ_value)
    println(f,"")
    println(f,"Power Flow:", f_value)
    println(f,"")
    global j=1
    for k in links
    println(f,"Power Flow line $k:", f_value[j])
    global j=j+1
    end
end

return (syscost_det,Dual_constraint2,Dual_constraint1,λ_value,f_value)
end


text:
mutable struct SessionVars
    values::Dict{Symbol, Any}
    last_accessed::Dates.DateTime
end

struct Sessions
    sessions::Dict{String, SessionVars}
    fresh_session::Dict{Symbol, Any}
end

SessionVars(value::Dict) = SessionVars(value, Dates.now())

function Base.getproperty(s::SessionVars, name::Symbol)
    setfield!(s, :last_accessed, Dates.now())
    getfield(s, :values)[name]
end

function Base.setproperty!(s::SessionVars, name::Symbol, x)
    setfield!(s, :last_accessed, Dates.now())
    getfield(s, :values)[name] = x
end

function expire_sessions!(sessions)
    for (id, session) in sessions.sessions
        if Dates.now() - getfield(session, :last_accessed) > Dates.Minute(5)
            delete!(sessions.sessions, id)
        end
    end
    sessions
end

function Base.getindex(sessions::Sessions, session_id::AbstractString)
    expire_sessions!(sessions)
    if !haskey(sessions.sessions, session_id)
        sessions.sessions[session_id] = SessionVars(copy(sessions.fresh_session))
    end
    sessions.sessions[session_id]
end

function create_sessions(app)
    if :register_session_vars in names(app, all = true)
        fresh_session = app.register_session_vars()
    else
        fresh_session = Dict()
    end

    Sessions(Dict{String, SessionVars}(), fresh_session)
end


text:
using Revise
using DataBench, RCall
using FastGroupBy
using IterableTables
using DataFrames
using CSV

function bench_gen_string_vec_var_len(n, strlen, skipbase = n >= 100_000_000)
    x = DataBench.gen_string_vec_var_len(n, strlen);
    bench_string_vec(x, n, strlen, "variable", skipbase)
end

function bench_gen_string_vec_fixed_len(n, strlen, skipbase = n >= 100_000_000)
    x = DataBench.gen_string_vec_fixed_len(n, strlen);
    bench_string_vec(x, n, strlen, "fixed", skipbase)
end

function bench_gen_string_vec_id_len(n, strlen, grps, skipbase = n >= 100_000_000)
    x = DataBench.gen_string_vec_id_fixed_len(n, strlen, grps)
    bench_string_vec(x, n, strlen, "id", skipbase)
end

function bench_string_vec(x, n, strlen, strlentype, skipbase = n >= 100_000_000)
    if !skipbase
        copyx = copy(x);
        julia_sort_elapsed = @elapsed sort!(copyx);
        @assert issorted(copyx)
    end

    copyx = copy(x);
    julia_lsd_radixsort_elapsed = @elapsed radixsort_lsd!(copyx);
    @assert issorted(copyx)

    copyx = copy(x);
    R"""
    memory.limit(2^32-1)
    rcopyx = $copyx
    summary(rcopyx)
    pt = proc.time()
    sort(rcopyx, method="radix")
    pt2 = proc.time()
    """;

    @rget pt;
    @rget pt2;

    r_radixsort_elapsed = pt2[3] - pt[3]

    c_qsort_strcmp_elapsed = @elapsed str_qsort!(copyx)
    @assert issorted(copyx)

    copyx = copy(x)
    three_way_radix_qsort_elapsed = @elapsed three_way_radix_qsort!(copyx)
    @assert issorted(copyx)

    if skipbase
        DataFrame(
        test = ["julia_lsd_radixsort_elapsed", "r_radixsort_elapsed", "c_qsort_strcmp_elapsed", "three_way_radix_qsort_elapsed"]
        ,elapsed = [julia_lsd_radixsort_elapsed, r_radixsort_elapsed, c_qsort_strcmp_elapsed, three_way_radix_qsort_elapsed]
        ,n = [n for i = 1:4]
        ,strlen = [strlen for i = 1:4]
        ,strlentype = [strlentype for i = 1:4]
        ,grps = [n÷100 for i=1:4])
    else
        DataFrame(
            test = ["julia_sort_elapsed", "julia_lsd_radixsort_elapsed", "r_radixsort_elapsed", "c_qsort_strcmp_elapsed", "three_way_radix_qsort_elapsed"]
            ,elapsed = [julia_sort_elapsed, julia_lsd_radixsort_elapsed, r_radixsort_elapsed, c_qsort_strcmp_elapsed, three_way_radix_qsort_elapsed]
            ,n = [n for i = 1:5]
            ,strlen = [strlen for i = 1:5]
            ,strlentype = [strlentype for i = 1:5]
            ,grps = [n÷100 for i=1:5])
    end
end

# bench_gen_string_vec_var_len(1_000_000, 8)
# bench_gen_string_vec_var_len(1_000_000, 16)
# bench_gen_string_vec_var_len(1_000_000, 24)
# bench_gen_string_vec_var_len(1_000_000, 32)
# bench_gen_string_vec_var_len(1_000_000, 64)

# bench_gen_string_vec_fixed_len(1_000_000, 8)
# bench_gen_string_vec_fixed_len(1_000_000, 16)
# bench_gen_string_vec_fixed_len(1_000_000, 24)
# bench_gen_string_vec_fixed_len(1_000_000, 32)
function withrepeat(x,y)
    [repeat(x, inner=[size(y,1)]) repeat(y, outer=[size(x,1)])]
end

ns = Int.(10.^(6:8))
strlens = [3, 10]
nsx = withrepeat(ns, strlens)

for i in 1:size(nsx,1)
    n = nsx[i,1]
    strlen = nsx[i,2]
    df = bench_gen_string_vec_id_len(n, 3, 100)
    println(df)
    CSV.write("benchmark/bench_string_results/string_sort_perf $(replace(string(now()),":"," ")).csv", df);
    df = bench_gen_string_vec_id_len(n, 10, n÷100)
    println(df)
    CSV.write("benchmark/bench_string_results/string_sort_perf $(replace(string(now()),":"," ")).csv", df);
end

ns = Int.(10.^(8:-1:6))
strlens = 8.*(8:-1:1)
nsx = withrepeat(ns, strlens)
df = DataFrame()
for i in 1:size(nsx,1)
    n = nsx[i,1]
    strlen = nsx[i,2]
    df = bench_gen_string_vec_var_len(n, strlen)
    CSV.write("benchmark/bench_string_results/string_sort_perf $(replace(string(now()),":"," ")).csv", df);
    println(df)
    df = bench_gen_string_vec_fixed_len(n, strlen)
    println(df)

    CSV.write("benchmark/bench_string_results/string_sort_perf $(replace(string(now()),":"," ")).csv", df);
end


text:
module SegwayRobot

    include("src/dynamics.jl")
    #export segway_rk4, segway_rk4!

    rk4_func  = segway_rk4  # set standard Environment function
    cost_func = segway_cost # set standard Environment cost-function
    #s_init = [0.0, 0.0, 0.0, 0.10] # set standard Environment initial position

    #State vector is [x0, y0, phi, alpha, dalpha, v, dphi, ir, il]
    const s_init  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    const u_range = [11.1, 11.1]
    export rk4_func, cost_func, s_init, u_range
end # module SegwayRoboter

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
"""
    Segment(trial, src::Union{S,Type{S},String}; [start, finish, conditions])

Describes a portion of a source in `trial` from time `start` to `finish` with segment
specific `conditions`, if applicable.

If `src` is a String or AbstractSource, it must refer to a source that exists in `trial`.
If the `start` is omitted, the segment will start at the beginning of the source/trial.
If the `finish` is omitted, the segment will be from time `start` to the end of the
source/trial.

Any conditions present in `trial` will be merged into the conditions for the segment.

# Example:

```julia
t = Trial(1, "intervention")

struct MySource <: AbstractSource
    path::String
end

seg = Segment(t, MySource; start=0.0, finish=10.0, conditions=Dict(:group => "control"))

# Use a source that already exists in the trial
seg5 = Segment(t, "main")
```
"""
struct Segment{S<:AbstractSource,ID}
    trial::Trial{ID}
    source::S
    start::Union{Nothing,Float64} # Start time
    finish::Union{Nothing,Float64} # End time
    conditions::Dict{Symbol}

    function Segment{S,ID}(
        trial::Trial{ID}, source::S, start, finish, conditions
    ) where {ID,S<:AbstractSource}
        if !isnothing(finish)
            start ≤ finish || throw(DomainError(finish,
                "finish time must be ≥ start time; got $finish"))
        end

        return new(trial, source, start, finish, merge(trial.conditions, conditions))
    end
end

function Segment(
    trial::Trial{ID},
    source::S;
    start::Union{Nothing,Float64}=nothing,
    finish::Union{Nothing,Float64}=nothing,
    conditions::Dict{Symbol}=Dict{Symbol,Any}()
) where {ID, S <: AbstractSource}
    return Segment{S,ID}(trial, source, start, finish, conditions)
end

function Segment(
    trial::Trial{ID},
    sourcename::String;
    kwargs...
) where ID
    hassource(trial, sourcename) ||
        throw(ArgumentError("source \"$sourcename\" not found in trial"))
    source = sources(trial)[sourcename]

    if (start = get(kwargs, :start, nothing)) isa Function
        start = start(trial)
    end
    if (finish = get(kwargs, :finish, nothing)) isa Function
        finish = finish(trial)
    end

    return Segment(trial, source; kwargs..., start, finish)
end

function Segment(trial, src::Type{<:AbstractSource}; kwargs...)
    Segment(trial, getsource(trial, src); kwargs...)
end

function Base.show(io::IO, s::Segment{S,ID}) where {S,ID}
    print(IOContext(io, :limit=>true), "Segment{$S,$ID}($(s.trial), ", typeof(s.source), "(…), ")
    print(io, isnothing(s.start) ? "begin" : s.start, ":")
    print(io, isnothing(s.finish) ? "end" : s.finish)
    if s.conditions == conditions(trial(s))
        print(io, ")")
    else
        print(io, ", (")
        _io = IOContext(io, :typeinfo=>eltype(s.conditions))
        first = true
        for p in pairs(s.conditions)
            first || print(_io, ", ")
            first = false
            print(_io, p)
        end
        print(io, "))")
    end
end

function Base.show(io::IO, mimet::MIME"text/plain", s::Segment{S,ID}) where {S,ID}
    println(IOContext(io, :limit=>true), "Segment{$S,$ID}\n Trial: $(s.trial)")
    println(io, " Source: ", s.source)
    print(io, " Time: ", isnothing(s.start) ? "beginning" : s.start)
    println(io, " to ", isnothing(s.finish) ? "end" : s.finish)
    print(io, " Conditions: ")
    if s.conditions == conditions(trial(s))
        print(io, "(same as parent trial)")
    end
    for c in s.conditions
        print(io, "\n    ")
        print(io, repr(c.first), " => ", repr(c.second))
    end
end

"""
    readsegment(seg::Segment{S}; warn=true, kwargs...) where S <: AbstractSource

If defined for `S`, returns the portion of `seg.source` from `seg.start` to `seg.finish`.
Otherwise, equivalent to `readsource` (i.e. no trimming of time-series occurs). Warns by
default if the main method is called.
"""
function readsegment(seg::Segment; warn=true, kwargs...)
    if warn
        @warn "A method of `readsegment` has not been defined for `$(typeof(source(seg)))`;"*
            "this means that `seg.start` and `seg.finish` will be ignored" maxlog=1
    end
    readsource(source(seg); kwargs...)
end

"""
    trial(seg::Union{Segment,SegmentResults}) -> Trial

Return the parent `Trial` for the given `Segment` or `SegmentResult`
"""
trial(seg::Segment) = seg.trial

"""
    subject(seg::Union{Segment,SegmentResults})

Return the subject reference for the parent `Trial` of the given `Segment` or
`SegmentResult`
"""
subject(seg::Segment) = subject(seg.trial)

"""
    source(seg::Union{Segment,SegmentResults}) -> AbstractSource

Return the source for the parent `Trial` of the given `Segment` or
`SegmentResult`
"""
source(seg::Segment) = seg.source

"""
    conditions(seg::Union{Segment,SegmentResults}) -> Dict{Symbol}

Return the conditions for the given `Segment` or `SegmentResult`
"""
conditions(seg::Segment) = seg.conditions

"""
    SegmentResult(segment::Segment, results::Dict{Symbol)

Contains the results of any analysis/analyses performed on the trial segment in a `Dict`.

# Example:

```julia
segresult = SegmentResult(seg, Dict(:avg => 3.5))
```
"""
struct SegmentResult{S,ID}
    segment::Segment{S,ID}
    results::Dict{String}
end

function SegmentResult(segment::Segment{S,ID}, results) where {S,ID}
    SegmentResult{S,ID}(segment, results)
end
SegmentResult(seg::Segment) = SegmentResult(seg, Dict{String,Any}())

trial(sr::SegmentResult) = trial(sr.segment)
subject(sr::SegmentResult) = subject(sr.segment)
source(sr::SegmentResult) = source(sr.segment)

"Get the segment of a `SegmentResult`"
segment(sr::SegmentResult) = sr.segment
conditions(sr::SegmentResult) = conditions(sr.segment)

"Get the results of a `SegmentResult`"
results(sr::SegmentResult) = sr.results

"""
    resultsvariables(sr::Union{SegmentResult,Vector{SegmentResult}})

Get the unique variables for `SegmentResult`s.
"""
resultsvariables(sr::SegmentResult) = collect(keys(results(sr)))
function resultsvariables(srs::Vector{<:SegmentResult})
    ks = reduce(vcat, collect.(unique(keys.(results.(srs)))))
    sort!(ks)
    unique!(ks)
end

function Base.show(io::IO, sr::SegmentResult)
    print(IOContext(io, :compact=>true, :limit=>true), "SegmentResult(",sr.segment,", ")
    if isempty(results(sr))
        print(io, "No results)")
    else
        print(IOContext(io, :limit=>true), "Results keys: ", keys(sr.results), ')')
    end
end

function Base.show(io::IO, ::MIME"text/plain", sr::SegmentResult{S,ID}) where {S,ID}
    println(io, "SegmentResult{$S,$ID}")
    print(io, ' ', sr.segment)
    println(io)
    show(io, MIME("text/plain"), sr.results)
end


text:
"""
    abstract_eval(ctx::LintContext, ex) :: Union{Any, Nothing}

Like `eval`, but does it in the current context and without any dynamism.
Returns `nothing` if the result can't be evaluated.
"""
function abstract_eval(ctx::LintContext, ex::Symbol)
    let lu = lookup(ctx, ex)
        if lu !== nothing
            extractobject(lu)
        end
    end
end

"""
    abstract_eval(ctx::LintContext, ex::Expr) :: Union{Any, Nothing}

If the given expression is curly, and each component of the curly is a constant
object in the given `ctx`, construct the object `x` as would have been done in
the program itself, and return `x`.

Otherwise, if the given expression is `foo.bar`, and `foo` is a standard
library object with attribute `bar`, then construct `foo.bar` as would be done
in the program itself and return it.

Otherwise, return `nothing`.
"""
abstract_eval(ctx::LintContext, ex::Expr) = begin
    if isexpr(ex, :curly)
        objs = [abstract_eval(ctx, arg) for arg in ex.args]
        if all(e->e!==nothing, objs)
            try
                Core.apply_type(objs...)
            catch
                nothing
            end
        else
            nothing
        end
    elseif isexpr(ex, :(.))
        head = ex.args[1]
        tail = ex.args[2].value
        obj = abstract_eval(ctx, head)
        if obj !== nothing
            try
                getfield(obj, tail)
            catch
                nothing
            end
        else
            nothing
        end
    else
        nothing
    end
end

"""
    abstract_eval(ctx::LintContext, ex)

Return the literal embedded within a `Union{Any, Nothing}`.
"""
abstract_eval(ctx::LintContext, ex) = lexicalvalue(ex)


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

# Wavelength range [nm]
λ = 200:1100
# Angle of incidence [degrees]
θ = [0.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.fused_silica_uv(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident),
    ModelFit(:forouhibloomer),
    LayerTMMO(emergent),
]

# Measured absolute transmittance
Texp = SpectraDB.scandia_spectrum(beam.λ)
plot(beam.λ, Texp)
gui()

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^3, show_trace=true, store_trace=false,
);

seed = [vcat(90.0, # thickness
        [1.6, 5.0, # ϵinf, Eg
         0.15, 16, 65, # A1, B1, C1
         ],
    ),
]

solOptim = fit_tmm_optics(
    Transmittance(Texp), seed, beam, layers;
    options=options, alg=LBFGS(),
)

plot(FitSpectrum(), solOptim.beam.λ, solOptim.spectrumExp, solOptim.spectrumFit)
gui()


text:
## -----------------  Deviation vector

"""
    DevVector

Deviation vector
"""
mutable struct DevVector
    dv :: Vector{AbstractDeviation}
end

"""
    DevVector()

Constructs an empty deviation vector
"""
function DevVector()
    DevVector(Vector{AbstractDeviation}())
end

Base.isempty(d :: DevVector) = Base.isempty(d.dv);
Base.length(d :: DevVector) = Base.length(d.dv);
Base.getindex(d :: DevVector, j) = d.dv[j];

# Iteration
function Base.iterate(d :: DevVector, j)
    if isempty(d)
        return nothing
    elseif j > length(d)
        return nothing
    else
        return d.dv[j], j+1
    end
end

Base.iterate(d :: DevVector) = Base.iterate(d, 1);

"""
    append!(d :: DevVector, dev :: Deviation)

Append a deviation.
"""
function Base.append!(d :: DevVector, dev :: AbstractDeviation)
    @assert !dev_exists(d, dev.name)  "Deviation $(dev.name) already exists"
    Base.push!(d.dv, dev)
end

"""
    $(SIGNATURES)

Set model values for one deviation in a `DevVector` specified by `name`
"""
function set_model_values(d :: DevVector, name :: Symbol, modelV)
    dev = retrieve(d, name);
    @assert !isempty(dev)
    set_model_values(dev, modelV);
    return nothing
end

"""
	$(SIGNATURES)

Get model values from a `DevVector`.
"""
function get_model_values(devV :: DevVector, mName :: Symbol)
    d = retrieve(devV, mName);
    return get_model_values(d)
end

"""
	$(SIGNATURES)

Get data values from a `DevVector`.
"""
function get_data_values(devV :: DevVector, mName :: Symbol)
    d = retrieve(devV, mName);
    return get_data_values(d)
end

"""
	set_weights!
"""
function set_weights!(d :: DevVector, name :: Symbol, wtV)
    dev = retrieve(d, name);
    @assert !isempty(dev)
    @assert size(wtV) == size(dev.dataV)  "Size mismatch: $(size(wtV)) vs $(size(dev.dataV))"
    set_weights!(dev, wtV);
    return nothing
end

"""
    retrieve

If not found: return empty Deviation
"""
function retrieve(d :: DevVector, dName :: Symbol)
    outDev = empty_deviation();

    n = length(d);
    if n > 0
        dIdx = 0;
        for i1 in 1 : n
            #println("$i1: $(d.dv[i1].name)")
            if d.dv[i1].name == dName
                dIdx = i1;
                break;
            end
            #println("  not found")
        end
        if dIdx > 0
            outDev = d.dv[dIdx];
        end
    end
    return outDev
end

function dev_exists(d :: DevVector, dName :: Symbol)
    return !isempty(retrieve(d, dName))
end

Base.show(io :: IO, d :: DevVector) = show_deviations(io, d);

"""
    $(SIGNATURES)

Show all deviations. Each gets a short display with name and scalar deviation.
"""
function show_deviations(io :: IO,  d :: DevVector; sorted :: Bool = false)
    if length(d) < 1
        println(io, "No deviations");
    else
        lineV = Vector{String}();
        for i1 in 1 : length(d)
            dStr = short_display(d.dv[i1]);
            push!(lineV, dStr);
        end
    end
    if sorted
        lineV = sort(lineV);
    end
    show_string_vector(lineV, 80, io = io);
end

show_deviations(d :: DevVector; kwargs...) = show_deviations(stdout, d; kwargs...);

"""
	$(SIGNATURES)

Return vector of scalar deviations.

Returns empty vector if `DevVector` is empty.
"""
function scalar_devs(d :: DevVector; inclScalarWt :: Bool = true)
    sds = scalar_dev_dict(d; inclScalarWt = inclScalarWt);
    return collect(values(sds));
end

"""
	$(SIGNATURES)

Make a `Dict{Symbol, DevType}` that maps deviation names into scalar deviations.
Useful for saving to disk.
"""
function scalar_dev_dict(d :: DevVector; inclScalarWt :: Bool = true)
    sds = Dict{Symbol, DevType}();
    n = length(d);
    if n > 0
        for i1 in 1 : n
            dev,_ = scalar_dev(d.dv[i1],  inclScalarWt = inclScalarWt);
            sds[name(d.dv[i1])] = dev;
        end
    end
    return sds
end

"""
	$(SIGNATURES)

Overall scalar deviation. Weighted sum of the scalar deviations returned by all `Deviation` objects
"""
function scalar_deviation(d :: DevVector)
    scalarDev = 0.0;
    for dev in d.dv
        sDev, _ = scalar_dev(dev,  inclScalarWt = true);
        @assert sDev >= 0.0
        scalarDev += sDev;
    end
    return scalarDev
end

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


text:
# By Reuben Hill 2019, Gowerlabs Ltd, reuben@gowerlabs.co.uk
#
# Copyright (c) Gowerlabs Ltd.

module TestUtil

using Test
using LibFTD2XX.Util

@testset "util" begin
  @test "hello" == ntuple2string(Cchar.(('h','e','l','l','o')))
  @test "hello" == ntuple2string(Cchar.(('h','e','l','l','o','\0','x')))

  @test v"0.0.0" == versionnumber(0)
  @test v"99.99.99" == versionnumber(0x00999999)
  @test v"2.12.28" == versionnumber(0x00021228)
  @test_throws DomainError versionnumber(0x00999999 + 1)
  @test_throws DomainError versionnumber(0x000000AA)
  @test_throws DomainError versionnumber(0x0000AA00)
  @test_throws DomainError versionnumber(0x00AA0000)
end

end


text:
using RDatasets, DataArrays, Gadfly

plot(dataset("car", "Womenlf"), x="HIncome", y="Region", Geom.histogram2d)


text:
# test staggered grid things

function test_staggered()

  fname = "input_vals_staggered.jl"
  mesh, sbp, eqn, opts = run_solver(fname)
  test_staggered_same(mesh, sbp, eqn, opts)

  delete!(opts, "use_staggered_grid")
  opts["operator_type"] = "SBPOmega"
  opts["operator_type2"] = "SBPGamma"
  opts["order"] = 2
  opts["order2"] = 3

  fname = make_input(opts, "input_vals_staggered2")

  mesh, sbp, eqn, opts = run_solver(fname)

  test_staggered_different(mesh, sbp, eqn, opts)

  return nothing
end

"""
  Test staggered grid functions when operators on both meshes are the samea

  This function is called from test_staggered()
"""
function test_staggered_same(mesh, sbp, eqn, opts)

  # force it to not use the staggered grid
  EulerEquationMod.dataPrep(mesh, sbp, eqn, opts)
  opts["use_staggered_grid"] = false
  # evalute the residual on the original, second mesh, make sure they are
  # the same
  mesh2 = mesh.mesh2
  sbp2 = mesh.sbp2

  println("testing primary mesh")
  fill!(eqn.res, 0.0)
  evalResidual(mesh, sbp, eqn, opts)
  res_orig = copy(eqn.res)

  println("testing secondary mesh")
  fill!(eqn.res, 0.0)
  evalResidual(mesh2, sbp2, eqn, opts)

  @test isapprox( vecnorm(res_orig - eqn.res), 0.0) atol=1e-13

  # test a non-uniform condition
  EulerEquationMod.ICExp(mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  opts["use_staggered_grid"] = true
  EulerEquationMod.dataPrep(mesh, sbp, eqn, opts)
  opts["use_staggered_grid"] = false

  fill!(eqn.res, 0.0)
  evalResidual(mesh, sbp, eqn, opts)
  res_orig = copy(eqn.res)

  fill!(eqn.res, 0.0)
  evalResidual(mesh2, sbp2, eqn, opts)

  @test isapprox( vecnorm(res_orig - eqn.res), 0.0) atol=1e-13

  # test staggered grid volume integrals
  println("testing volume integrals")
  fill!(eqn.res, 0.0)
  EulerEquationMod.calcVolumeIntegralsSplitFormCurvilinear(mesh, sbp, eqn, opts,
                                                           eqn.flux_func)

  println("testing staggered volume integrals")
  res_orig = copy(eqn.res)
  fill!(eqn.res, 0.0)
  EulerEquationMod.calcVolumeIntegralsSplitFormCurvilinear(mesh, mesh2, sbp, sbp2, eqn, opts, eqn.flux_func)

  @test isapprox( vecnorm(eqn.res - res_orig), 0.0) atol=1e-13

  fill!(eqn.res, 0.0)
  EulerEquationMod.getFaceElementIntegral(mesh, sbp, eqn, eqn.face_element_integral_func, eqn.flux_func, mesh.sbpface, mesh.interfaces)

  res_orig = copy(eqn.res)
  fill!(eqn.res, 0.0)
  EulerEquationMod.getFaceElementIntegral(mesh, mesh2, sbp, sbp2, eqn, eqn.face_element_integral_func, eqn.flux_func, mesh2.sbpface, mesh.interfaces)

  @test isapprox( vecnorm(eqn.res - res_orig), 0.0) atol=1e-13

  fill!(eqn.res, 0.0)

  return nothing
end

"""
  Test staggered grid functions when the meshes are different

  This function is called from test_staggered()
"""
function test_staggered_different(mesh, sbp, eqn, opts)

  mesh2 = mesh.mesh2
  sbp2 = mesh.sbp2

  # test interpolation
  # the operators are
  q_ex = zeros(5, mesh.numNodesPerElement)
  q = zeros(q_ex)
  q2 = zeros(5, mesh.mesh2.numNodesPerElement)
  q2_ex = zeros(q2)

  # sbp_s is 2nd order, so make sure a 2nd order polynomial is interpolated
  # exactly
  for i=1:mesh.numNodesPerElement
    x = mesh.coords[1, i, 1]
    y = mesh.coords[2, i, 1]

    q_ex[1, i] = 1
    q_ex[2, i] = x + 1
    q_ex[3, i] = y + 1
    q_ex[4, i] = x + y + 1
    q_ex[5, i] = x*x + y*y + x + y + 1
  end

  for i=1:mesh2.numNodesPerElement
    x = mesh2.coords[1, i, 1]
    y = mesh2.coords[2, i, 1]

    q2_ex[1, i] = 1
    q2_ex[2, i] = x + 1
    q2_ex[3, i] = y + 1
    q2_ex[4, i] = x + y + 1
    q2_ex[5, i] = x*x + y*y + x + y + 1
  end

  smallmatmat!(q_ex, mesh.I_S2FT, q2)
  @test isapprox( norm(q2 - q2_ex), 0.0) atol=1e-13

  smallmatmat!(q2_ex, mesh.I_F2ST, q)
  @test isapprox( norm(q - q_ex), 0.0) atol=1e-13

  #=
  # test staggered grid volume integrals
  # S is skew-symmetric, so c.' S c = 0
  # this is basically the test of conservation
  EulerEquationMod.ICRhoe1E2U3(mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  q_node = eqn.q[:, 1, 1]

  res_orig = copy(eqn.res)
  =#

  EulerEquationMod.dataPrep(mesh, sbp, eqn, opts)

  fill!(eqn.res, 0.0)
  EulerEquationMod.calcVolumeIntegralsSplitFormCurvilinear(mesh, mesh2, sbp, sbp2, eqn, opts, eqn.flux_func)
  @test isapprox( sum(eqn.res), 0.0) atol=1e-13

  fill!(eqn.res, 0.0)
  applyPoly(mesh, sbp, eqn, opts, sbp.degree)
  penalty_functor = EulerEquationMod.FaceElementDict["ELFPenaltyFaceIntegral"](mesh, eqn)

  EulerEquationMod.getFaceElementIntegral(mesh, mesh2, sbp, sbp2, eqn, penalty_functor, eqn.flux_func, mesh2.sbpface, mesh.interfaces)

  @test isapprox( vecnorm(eqn.res), 0.0) atol=1e-13

  return nothing
end

add_func1!(EulerTests, test_staggered, [TAG_SHORTTEST, TAG_STAGGERED])


text:
using TerraDG

configfile = "./advection.yaml"
TerraDG.main(configfile)


text:
#Function for calculating matrix elements ⟨n₁,l₁,j₁|er|n₂,l₂,j₂⟩ between Rydberg states for a given atom.
#23/07/2017

function matrix_elements(atom,nn1,ll1,jj1,nn2,ll2,jj2)

  matrixelement = radiel(atom,nn1,ll1,jj1,nn2,ll2,jj2)
  matrixelementSI = matrixelement*bohrrad*eleccharge

  return matrixelement, matrixelementSI
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
using Flux
using Zygote

import Zygote:gradient

gradient(f, ::Val{:Zygote}, args...) = gradient(f, args...)

Zygote.@adjoint argmax(xs; dims = :) = argmax(xs;dims=dims), _ -> nothing

# ??? can safely removed now
Zygote.@adjoint function Base.broadcasted(::typeof(relu), x::Array{T}) where T<:Real
    y = relu.(x)
    return y, Δ -> begin
        res = similar(Δ)
        for i in 1:length(res)
            if y[i] > 0
                res[i] = Δ[i]
            else
                res[i] = zero(T)
            end
        end
        (nothing, res)
    end
end


text:
import POMDPs: MDP, Solver

type EVSarsaPolicySolver <: Solver
    epochs::Int64
    epsilon::Float64
    epsilon_discount::Float64
    learning_rate::Float64
    discount::Float64
    max_frame_iterations::Int64
    print_every_n::Int64
end

EVSarsaPolicySolver(epochs; epsilon = 0.25, epsilon_discount = 0.99, learning_rate = 0.8, discount = 0.8, max_frame_iterations = 100, print_every_n = 25) =
    EVSarsaPolicySolver(epochs, epsilon, epsilon_discount, learning_rate, discount, max_frame_iterations, print_every_n)

function solve(solver::EVSarsaPolicySolver, pomdp::MDP; verbose = true)

    policy = ValueTablePolicy(pomdp, solver.epsilon, solver.learning_rate, solver.discount)
    last_rewards = zeros(Float64, solver.epochs)
    actions_count = n_actions(pomdp)

    for i=1:solver.epochs

        if verbose && (i == 1 || i % solver.print_every_n == 0)
            println("Epoch: $i")
        end

        #
        # SARSA implementation that updates on every step
        #
        previous_state::Int64 = initial_state(pomdp)
        rewards = 0.0

        for j=1:solver.max_frame_iterations

            current_action_i = action(policy, previous_state)
            current_action = action_index(pomdp, current_action_i)
            current_state, current_reward = generate_sr(pomdp, previous_state, current_action)
            current_state_terminal = isterminal(pomdp, current_state)

            # expected value sarsa sum all values from a next state
            next_value = if current_state_terminal 0 else  sum(policy.value_map[:, current_state]) end
            new_value = current_reward + solver.discount * next_value / actions_count
            new_updated_value = solver.learning_rate * new_value + (1 - solver.learning_rate) * policy.value_map[current_action_i, previous_state]

            policy.value_map[current_action_i, previous_state] = new_updated_value
            previous_state = current_state
            rewards += current_reward

            if current_state_terminal
                break
            end
        end

        if 0 < solver.epsilon_discount < 1
            policy.epsilon *= solver.epsilon_discount
        end

        last_rewards[i] = rewards

        if verbose && (i > 100 && i % solver.print_every_n == 0)
            cum_sum_reward = sum(last_rewards[i-101:i])
            cum_avg_reward = mean(last_rewards[i-101:i])
            println("Epoch: $i, Cummulative reward: $(round(cum_avg_reward, 2))/$(round(cum_sum_reward, 2)), Epsilon: $(policy.epsilon)")
        end
    end

    return policy

end

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
######################################################################
# Functionals UTs
# -----
# Licensed under MIT License
module TestFunctionals
using Test
using KirUtil
import ..Helpers: Immutable

@testset "KirUtil Functionals" begin
    @testset "indexed" begin
        let ary = [1, 2, 3]
            @test indexed(ary) === ary
        end
        let tpl = (1, 2, 3)
            @test indexed(tpl) === tpl
        end
        let set = Set((1, 2, 3)), idxed = indexed(set)
            @test 1 ∈ idxed && 2 ∈ idxed && 3 ∈ idxed
            @test set !== idxed
        end
        let idxed = indexed(nothing)
            @test length(idxed) == 1
            @test idxed[1] === nothing
        end
        let idxed = indexed(42)
            @test length(idxed) == 1
            @test idxed[1] == 42
        end
    end

    @testset "iterable" begin
        @test iterable(42) === 42
        let ary = [1, 2, 3]
            @test iterable(ary) === ary
        end
        let tpl = (1, 2, 3)
            @test iterable(tpl) === tpl
        end
        let set = Set((1, 2, 3))
            @test iterable(set) === set
        end
        let itr = iterable(Ident{42}())
            @test itr isa Tuple
            @test hassignature(iterate, typeof(itr))
        end
    end

    @testset "split collections" begin
        @test split(iseven, collect(1:10)) == (collect(2:2:10), collect(1:2:9))
        @test split(iseven, tuple(1:10...)) == (tuple(2:2:10...), tuple(1:2:9...))
    end

    @testset "functional Base.insert!" begin
        let ary = [1, 3, 5]
            @test insert!(ary, 2, before=3) == [1, 2, 3, 5]
            @test insert!(ary, 4, after=3) == [1, 2, 3, 4, 5]
        end

        let ary = Immutable.([1, 4, 6])
            @test insert!(ary, Immutable(5), before=6, by=(im)->im.mutable[]) == Immutable.([1, 4, 5, 6])
        end
    end

    @testset "functionals" begin
        @test !isathing(nothing)
        @test isathing(42)

        @test negate(isodd)(42) == iseven(42)
        @test negate(iseven)(42) == isodd(42)

        @test truthy(true)
        @test truthy(42)
        @test truthy(69.69)
        @test truthy(1//2)
        @test falsy(false)
        @test falsy(0)
        @test falsy(0.0)
        @test !falsy(-1)
        @test !falsy(1//1)
    end

    @testset "bytes" begin
        @test bytes(UInt8(42)) == [0x2A]
        @test bytes(Int8(69)) == [0x45]
        @test bytes(hton(UInt16(42))) == [0x00, 0x2A]
        @test bytes(hton(Int16(69))) == [0x00, 0x45]
        @test bytes(hton(Int32(69696969))) == [0x04, 0x27, 0x7D, 0xC9]
        @test bytes(hton(Int64(6969696969))) == [0x00, 0x00, 0x00, 0x01, 0x9F, 0x6D, 0x22, 0xC9]

        @test bytes(0.42069) == [0x78, 0x45, 0xF0, 0xBF, 0x95, 0xEC, 0xDA, 0x3F]

        @test bytes([0x1, 0x2, 0x3]) == [0x01, 0x02, 0x03]
        @test bytes((0x1, 0x2, 0x3)) == [0x01, 0x02, 0x03]
    end
end
end # module TestFunctionals


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct StorageAccountPropertiesUpdateParameters <: SwaggerModel
    customDomain::Any # spec type: Union{ Nothing, CustomDomain } # spec name: customDomain
    encryption::Any # spec type: Union{ Nothing, Encryption } # spec name: encryption
    accessTier::Any # spec type: Union{ Nothing, String } # spec name: accessTier
    supportsHttpsTrafficOnly::Any # spec type: Union{ Nothing, Bool } # spec name: supportsHttpsTrafficOnly
    networkAcls::Any # spec type: Union{ Nothing, NetworkRuleSet } # spec name: networkAcls

    function StorageAccountPropertiesUpdateParameters(;customDomain=nothing, encryption=nothing, accessTier=nothing, supportsHttpsTrafficOnly=false, networkAcls=nothing)
        o = new()
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("customDomain"), customDomain)
        setfield!(o, Symbol("customDomain"), customDomain)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("encryption"), encryption)
        setfield!(o, Symbol("encryption"), encryption)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("accessTier"), accessTier)
        setfield!(o, Symbol("accessTier"), accessTier)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("supportsHttpsTrafficOnly"), supportsHttpsTrafficOnly)
        setfield!(o, Symbol("supportsHttpsTrafficOnly"), supportsHttpsTrafficOnly)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("networkAcls"), networkAcls)
        setfield!(o, Symbol("networkAcls"), networkAcls)
        o
    end
end # type StorageAccountPropertiesUpdateParameters

const _property_map_StorageAccountPropertiesUpdateParameters = Dict{Symbol,Symbol}(Symbol("customDomain")=>Symbol("customDomain"), Symbol("encryption")=>Symbol("encryption"), Symbol("accessTier")=>Symbol("accessTier"), Symbol("supportsHttpsTrafficOnly")=>Symbol("supportsHttpsTrafficOnly"), Symbol("networkAcls")=>Symbol("networkAcls"))
const _property_types_StorageAccountPropertiesUpdateParameters = Dict{Symbol,String}(Symbol("customDomain")=>"CustomDomain", Symbol("encryption")=>"Encryption", Symbol("accessTier")=>"String", Symbol("supportsHttpsTrafficOnly")=>"Bool", Symbol("networkAcls")=>"NetworkRuleSet")
Base.propertynames(::Type{ StorageAccountPropertiesUpdateParameters }) = collect(keys(_property_map_StorageAccountPropertiesUpdateParameters))
Swagger.property_type(::Type{ StorageAccountPropertiesUpdateParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_StorageAccountPropertiesUpdateParameters[name]))}
Swagger.field_name(::Type{ StorageAccountPropertiesUpdateParameters }, property_name::Symbol) =  _property_map_StorageAccountPropertiesUpdateParameters[property_name]

const _allowed_StorageAccountPropertiesUpdateParameters_accessTier = ["Hot", "Cool"]

function check_required(o::StorageAccountPropertiesUpdateParameters)
    true
end

function validate_property(::Type{ StorageAccountPropertiesUpdateParameters }, name::Symbol, val)
    if name === Symbol("accessTier")
        Swagger.validate_param(name, "StorageAccountPropertiesUpdateParameters", :enum, val, _allowed_StorageAccountPropertiesUpdateParameters_accessTier)
    end
end


text:
using SafeTestsets
using Test

module SampleVAS
using Combinatorics
using Distributions
function sample_transitions()
    all = [
        -1  0  1  0  0;
         1 -1 -1 -1  1;
         0  1  0  1 -1
    ]
    take = -1 * (all .* (all .< 0))
    give = all .* (all .> 0)
    rates = Function[x -> Exponential(1.0) for x in 1:size(all, 2)]
    (take, give, rates)
end

# cnt is the number of individuals
# vector is cnt of S, cnt of I, cnt of R
function sample_sir(cnt)
    s = 0
    i = cnt
    r = 2 * cnt
    i_rate, r_rate = (0.9, 0.01)
    infection = collect(combinations(1:cnt, 2))
    recover = cnt
    # Rows are s individuals, i individuals, r individuals
    # Columns are infections, then recoveries.
    take = zeros(Int, 3 * cnt, 2 * length(infection) + recover)
    give = zeros(Int, 3 * cnt, 2 * length(infection) + recover)
    rates = Function[]
    idx = 0
    for (a, b) in infection
        idx += 1  # b infects a
        take[s + a, idx] = 1
        take[i + b, idx] = 1
        give[i + a, idx] = 1
        give[i + b, idx] = 1
        push!(rates, state -> Exponential(i_rate))
        idx += 1  # a infects b
        take[i + a, idx] = 1
        take[s + b, idx] = 1
        give[i + a, idx] = 1
        give[i + b, idx] = 1
        push!(rates, state -> Exponential(i_rate))
    end
    for c in 1:recover
        idx += 1
        take[i + c, idx] = 1
        give[r + c, idx] = 1
        push!(rates, state -> Exponential(r_rate))
    end
    (take, give, rates)
end
end

@safetestset vas_creates = "VAS creates" begin
using Fleck: VectorAdditionSystem
using ..SampleVAS: sample_transitions
    take, give, rates = sample_transitions()
    vas = VectorAdditionSystem(take, give, rates)
    @test length(vas.rates) == size(vas.take, 2)
end

@safetestset vas_intializes = "VAS initializes state" begin
using Fleck: VectorAdditionSystem, vas_initial, zero_state
using ..SampleVAS: sample_transitions
    take, give, rates = sample_transitions()
    vas = VectorAdditionSystem(take, give, rates)
    initializer = vas_initial(vas, [1, 1, 0])
    state = zero_state(vas)
    initializer(state)
    @test state == [1, 1, 0]
end

@safetestset vas_reports_hazards = "VAS reports_hazards" begin
using Fleck: VectorAdditionSystem, vas_initial
using Fleck: zero_state, push!, fire!
using ..SampleVAS: sample_transitions
take, give, rates = sample_transitions()
vas = VectorAdditionSystem(take, give, rates)
initializer = vas_initial(vas, [1, 1, 0])
    state = zero_state(vas)
    disabled = zeros(Int, 0)
    enabled = zeros(Int, 0)
    track_hazards = (idx, dist, enable, rng) -> begin
        if enable == :Enabled
            push!(enabled, idx)
        elseif enable == :Disabled
            push!(disabled, idx)
        else
            @test enable in [:Enabled, :Disabled]
        end
    end
    fire!(track_hazards, vas, state, initializer, "rng")
    @test enabled == [1, 2, 3, 4]
    @test length(disabled) == 0
end

@safetestset vas_changes_state = "VAS hazards during state change" begin
using Fleck: VectorAdditionSystem, vas_input, vas_initial
using Fleck: zero_state, push!, fire!
using ..SampleVAS: sample_transitions
take, give, rates = sample_transitions()
vas = VectorAdditionSystem(take, give, rates)
initializer = vas_initial(vas, [1, 1, 0])
    state = zero_state(vas)
    initializer(state)
    disabled = zeros(Int, 0)
    enabled = zeros(Int, 0)
    track_hazards = (idx, dist, enable, rng) -> begin
        if enable == :Enabled
            push!(enabled, idx)
        elseif enable == :Disabled
            push!(disabled, idx)
        else
            @test enable in [:Enabled, :Disabled]
        end
    end
    fire_index = 2
    input_change = vas_input(vas, fire_index)
    fire!(track_hazards, vas, state, input_change, "rng")
    @test enabled == [5]
    @test disabled == [2, 3, 4]
end

@safetestset vas_loops = "VAS main loop" begin
using Fleck: VectorAdditionSystem, vas_input, vas_initial
using Fleck: zero_state, push!, fire!
using Random: MersenneTwister
using ..SampleVAS: sample_transitions
    rng = MersenneTwister(2930472)
    take, give, rates = sample_transitions()
    vas = VectorAdditionSystem(take, give, rates)
    state = zero_state(vas)
    next_transition = nothing
    disabled = zeros(Int, 0)
    enabled = zeros(Int, 0)
    newly_enabled = zeros(Int, 0)
    track_hazards = (idx, dist, enable, rng) -> begin
        if enable == :Enabled
            push!(newly_enabled, idx)
        elseif enable == :Disabled
            push!(disabled, idx)
        else
            @test enable in [:Enabled, :Disabled]
        end
    end
    for i in 1:10
        if isnothing(next_transition)
            input_change = vas_initial(vas, [1, 1, 0])
        else
            input_change = vas_input(vas, next_transition)
        end
        fire!(track_hazards, vas, state, input_change, rng)
        if any(state .< 0)
            println("idx $(i) state $(state)")
        end
        @test all(state .>= 0)

        filter!(x -> x ∉ disabled, enabled)
        append!(enabled, newly_enabled)

        if length(enabled) > 0
            next_transition = rand(rng, enabled)
        else
            println("Nothing enabled, $(state)")
            next_transition = nothing
        end
        empty!(newly_enabled)
        empty!(disabled)
    end
end

@safetestset vas_fsm_init = "VAS finite state init" begin
    using Fleck: VectorAdditionFSM, vas_initial, VectorAdditionSystem
    using Fleck: VectorAdditionModel, VectorAdditionFSM, zero_state
    using Fleck: DirectCall, simstep!, vas_input
    using Random: MersenneTwister
    using ..SampleVAS: sample_transitions
    rng = MersenneTwister(2930472)
    take, give, rates = sample_transitions()
    vas = VectorAdditionSystem(take, give, rates)
    state = zero_state(vas)
    vam = VectorAdditionModel(vas, state, 0.0)
    fsm = VectorAdditionFSM(vam, DirectCall(Int))
    initial_state = vas_initial(vas, [1, 1, 0])
    when, next_transition = simstep!(fsm, initial_state, rng)
    limit = 10
    while next_transition !== nothing && limit > 0
        token = vas_input(vas, next_transition)
        when, next_transition = simstep!(fsm, token, rng)
        limit -= 1
    end
end

@safetestset vas_fsm_sir = "VAS runs SIR to completion" begin
    using Fleck: VectorAdditionFSM, vas_initial, VectorAdditionSystem
    using Fleck: VectorAdditionModel, VectorAdditionFSM, zero_state
    using Fleck: DirectCall, simstep!, vas_input
    using Random: MersenneTwister
    using ..SampleVAS: sample_sir
    rng = MersenneTwister(979797)
    cnt = 10
    take, give, rates = sample_sir(cnt)
    vas = VectorAdditionSystem(take, give, rates)
    state = zero_state(vas)
    vam = VectorAdditionModel(vas, state, 0.0)
    fsm = VectorAdditionFSM(vam, DirectCall(Int))
    starting = zeros(Int, 3 * cnt)
    starting[2:cnt] .= 1
    starting[cnt + 1] = 1  # Start with one infected.
    initial_state = vas_initial(vas, starting)
    when, next_transition = simstep!(fsm, initial_state, rng)
    event_cnt = 0
    while next_transition !== nothing
        token = vas_input(vas, next_transition)
        when, next_transition = simstep!(fsm, token, rng)
        event_cnt += 1
    end
    println(event_cnt)
end


text:
module PandorasBoxColoredBigWaterfall

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/coloredBigWaterfall" ColoredBigWaterfall(x::Integer, y::Integer, width::Integer=16, height::Integer=64, layer::String="FG", color::String="LightSkyBlue")

const placements = Ahorn.PlacementDict(
    "Big Waterfall (FG) (Pandora's Box)" => Ahorn.EntityPlacement(
        ColoredBigWaterfall,
        "rectangle",
        Dict{String, Any}(
            "layer" => "FG"
        )
    ),
    "Big Waterfall (BG) (Pandora's Box)" => Ahorn.EntityPlacement(
        ColoredBigWaterfall,
        "rectangle",
        Dict{String, Any}(
            "layer" => "BG"
        )
    )
)

const segmentCache = Dict{String, Any}()

const waterSegmentLeftMatrix = [
    1 1 1 0 1 0;
    1 1 1 0 1 0;
    1 1 1 0 1 0;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 0 1 0;
    1 1 1 0 1 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
]

const waterSegmentRightMatrix = [
    0 1 0 1 1 1;
    0 1 0 1 1 1;
    0 1 0 1 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 1 0 1 1 1;
    0 1 0 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
]

function getColor(color)
    if haskey(Ahorn.XNAColors.colors, color)
        return Ahorn.XNAColors.colors[color]

    else
        try
            return ((Ahorn.argb32ToRGBATuple(parse(Int, replace(color, "#" => ""), base=16))[1:3] ./ 255)..., 1.0)

        catch

        end
    end

    return (1.0, 1.0, 1.0, 1.0)
end

function getSegments(color)
    if haskey(segmentCache, color)
        return segmentCache[color]

    else
        baseColor = getColor(color)
        fillColor = baseColor .* 0.3
        surfaceColor = baseColor .* 0.8

        waterSegmentRight = Ahorn.matrixToSurface(
            waterSegmentRightMatrix,
            [
                fillColor,
                surfaceColor
            ]
        )

        waterSegmentLeft = Ahorn.matrixToSurface(
            waterSegmentLeftMatrix,
            [
                fillColor,
                surfaceColor
            ]
        )

        segmentCache[color] = (waterSegmentLeft, waterSegmentRight)

        return segmentCache[color]
    end
end

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

Ahorn.minimumSize(entity::ColoredBigWaterfall) = 8, 8
Ahorn.resizable(entity::ColoredBigWaterfall) = true, true

Ahorn.selection(entity::ColoredBigWaterfall) = Ahorn.getEntityRectangle(entity)

Ahorn.editingOptions(entity::ColoredBigWaterfall) = Dict{String, Any}(
    "layer" => String["FG", "BG"],
    "color" => colors
)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ColoredBigWaterfall, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 16))
    height = Int(get(entity.data, "height", 64))

    color = get(entity.data, "color", "LightSkyBlue")
    baseColor = getColor(color)
    fillColor = baseColor .* 0.3

    segmentHeightLeft, segmentWidthLeft = size(waterSegmentLeftMatrix)
    segmentHeightRight, segmentWidthRight = size(waterSegmentRightMatrix)

    Ahorn.Cairo.save(ctx)

    Ahorn.rectangle(ctx, 0, 0, width, height)
    Ahorn.clip(ctx)

    waterSegmentLeft, waterSegmentRight = getSegments(color)

    for i in 0:segmentHeightLeft:ceil(Int, height / segmentHeightLeft) * segmentHeightLeft
        Ahorn.drawImage(ctx, waterSegmentLeft, 0, i)
        Ahorn.drawImage(ctx, waterSegmentRight, width - segmentWidthRight, i)
    end

    # Drawing a rectangle normally doesn't guarantee that its the same color as above
    if height >= 0 && width >= segmentWidthLeft + segmentWidthRight
        fillRectangle = Ahorn.matrixToSurface(fill(0, (height, width - segmentWidthLeft - segmentWidthRight)), [fillColor])
        Ahorn.drawImage(ctx, fillRectangle, segmentWidthLeft, 0)
    end

    Ahorn.restore(ctx)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
using ITensors
using PastaQ
using Test

@testset "productstate" begin
  s = siteinds("Qubit", 4)

  states0 = [state("0", s[n]) for n in 1:length(s)]
  states1 = [state("1", s[n]) for n in 1:length(s)]
  statesXp = [state("X+", s[n]) for n in 1:length(s)]
  statesYm = [state("Y-", s[n]) for n in 1:length(s)]

  ψ0 = productstate(s)
  @test ψ0 isa MPS
  @test eltype(ψ0) == ITensor
  @test PastaQ.promote_leaf_eltypes(ψ0) == Float64
  @test ψ0 ≈ prod(states0)

  @test productstate(s, 0) ≈ ψ0
  @test productstate(s, "0") ≈ ψ0

  ψXp = productstate(s, "X+")
  @test ψXp isa MPS
  @test PastaQ.promote_leaf_eltypes(ψXp) == Float64
  @test ψXp ≈ prod(statesXp)

  ψYm = productstate(s, "Y-")
  @test ψYm isa MPS
  @test PastaQ.promote_leaf_eltypes(ψYm) == ComplexF64
  @test prod(ψYm) ≈ prod(statesYm)

  @test productstate(s, fill("0", length(s))) ≈ ψ0
  @test productstate(s, fill(0, length(s))) ≈ ψ0
  @test productstate(s, fill("X+", length(s))) ≈ prod(statesXp)
  @test productstate(s, fill("Y-", length(s))) ≈ prod(statesYm)

  @test productstate(s, n -> isodd(n) ? 1 : 0) ≈
        states1[1] * states0[2] * states1[3] * states0[4]
  @test productstate(s, n -> isodd(n) ? "1" : "0") ≈
        states1[1] * states0[2] * states1[3] * states0[4]
  @test productstate(s, n -> isodd(n) ? "Y-" : "X+") ≈
        statesYm[1] * statesXp[2] * statesYm[3] * statesXp[4]

  ψ = runcircuit(s, ghz(4))
  ϕ = productstate(ψ, fill("0", length(s))) ≈ ψ0

  n = 4
  ψ = productstate(n; dim = 3)
  @test length(ψ) == n
  @test length(PastaQ.array(ψ)) == 3^n

  n = 2
  s = qudits(n; dim = 3)
  ψ = productstate(n, [0,2]; dim = 3)
  @test PastaQ.array(ψ) == [0,0,1,0,0,0,0,0,0]
  ψ = productstate(s, [2,1])
  @test PastaQ.array(ψ) == [0,0,0,0,0,0,0,1,0]

end

@testset "productoperator" begin
  s = siteinds("Qubit", 4)

  gatesI = [gate("Id", s, n) for n in 1:length(s)]

  I = productoperator(s)
  @test PastaQ.promote_leaf_eltypes(I) == Float64
  @test prod(I) ≈ prod(gatesI)
  @test I ≈ prod(gatesI)

  U = runcircuit(s, ghz(4); process = true)
  X = productoperator(U)
  @test X ≈ I

  s = qudits(3; dim = 3)
  X = productoperator(s)
  @test PastaQ.array(X) ≈ PastaQ.array(productoperator(3; dim = 3))
  @test PastaQ.array(X) ≈ Matrix(LinearAlgebra.I,27,27)
end


text:
include("DataFile.jl")
include("SolveBalances.jl")
using Debug

## Keep track of datafile indexes

# TGFB dosing
tgfb12_initial_cond_index = 76 # TGFB12
tgfb12_forward_rate_index = 98 # addition_TGFb12: [] =([])=> TGFB12
tgfb12_reverse_rate_index = 99 # addition_TGFb12: TGFB12 =([])=> [] (reverse)
tgfb12_binding_rate_index = 1 # TGFB12_binding: TGFB12+R =([])=> TGFB12_R

# TGFB signaling
tgfb3_binding_rate_index = 3 # TGFB3_binding: TGFB3+R =([])=> TGFB3_R
R_initial_cond_index = 75
R_deg_rate_index = 51 # R_degradation: R =([])=> []
TGFB12_R_activation_RAS_rate_index = 14 #TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
RAS_activation_RAF_rate_index  = 15 # RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
RAF_activation_MAPK_rate_index = 16 # RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
MAPK_activation_AP1_SP1_rate_index = 18 # MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
tgfb12_to_ap1_sp1_gain = 1.5

# Bcatenin complex formation

bcat_ecad_forward_index = 19 	# 19	Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
bcat_ecad_reverse_index = 20	# 20	Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
bcat_tcf4_forward_index = 21	# 21	TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
bcat_tcf4_reverse_index = 22	# 22	TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
bcat_lef1_forward_index = 23 	# 23	LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
bcat_lef1_reverse_index = 24	# 24	LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)

bcat_deg_rate_index = 37 # Bcatenin_degradation: Bcatenin =([])=> []
bcat_deg_rate = 0.3 # default 0.1

ecad_bais_factor = 2.0 #
bcat_complex_gain = 5.0

# VEGF dosing
vegf_initial_cond_index = 78 # VEGF
vegf_forward_rate_index = 100 # addition_VEGF: [] =([])=> VEGF
vegf_reverse_rate_index = 101 # addition_VEGF: VEGF =([])=> [] (reverse)
vegf_binding_rate_index = 5 # VEGF+VEGFR =([])=> VEGF_VEGFR

# VEGF signaing
PKC_activation_RAF_rate_index = 11 # PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P

# =============================================================================== #
# Solve.jl
# Function library to solve the EMT model for different experimental/test cases
# Author: J.Varner
#
# Version: 1.0
# Date: 12/17/15
# =============================================================================== #

# Simulate the addition of TGFB -
function SolveModelAddTGFB(TSTART,TSTOP,Ts,data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant

  # Modify rate constant array
  rate_constant_array=data_dictionary["RATE_CONSTANT_ARRAY"];
  # Kill TGFB signal for steady state
  rate_constant_array[tgfb12_forward_rate_index]=0;   # forward [] =([])=> TGFB12
  rate_constant_array[tgfb12_reverse_rate_index]=0;  # Reverse TGFB12 =([])=> [] (reverse)
  rate_constant_array[tgfb12_binding_rate_index]=3.0 # speed up TGFB12 binding
  rate_constant_array[tgfb3_binding_rate_index]=3.0 # speed up TGFB3 binding
  rate_constant_array[vegf_forward_rate_index]=0;   # forward [] =([])=> VEGF
  rate_constant_array[vegf_reverse_rate_index]=0;  # Reverse VEGF =([])=> [] (reverse)
  rate_constant_array[vegf_binding_rate_index]=3.0; # speed up VEGF binding
  rate_constant_array[TGFB12_R_activation_RAS_rate_index]=2.0 # speed up RAS activation due to Receptor delay
  rate_constant_array[PKC_activation_RAF_rate_index]=0.1 # bias towards NFATc1 (10x)
  rate_constant_array[bcat_ecad_forward_index]= bcat_complex_gain*ecad_bais_factor
  rate_constant_array[bcat_ecad_reverse_index]= 1.0
  rate_constant_array[bcat_tcf4_forward_index]= bcat_complex_gain
  rate_constant_array[bcat_tcf4_reverse_index]= 1.0
  rate_constant_array[bcat_lef1_forward_index]= bcat_complex_gain
  rate_constant_array[bcat_lef1_reverse_index]= 1.0
  rate_constant_array[bcat_deg_rate_index]= bcat_deg_rate
  rate_constant_array[RAS_activation_RAF_rate_index] = tgfb12_to_ap1_sp1_gain
  rate_constant_array[RAF_activation_MAPK_rate_index] = tgfb12_to_ap1_sp1_gain
  rate_constant_array[MAPK_activation_AP1_SP1_rate_index] = tgfb12_to_ap1_sp1_gain
  #sequestation Bcatenin rate
  data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array;

  # Modify CONTROL_PARAMETER_ARRAY
  GAIN_MASTER = 0.2; # default is 0.1
  control_parameter_array=data_dictionary["CONTROL_PARAMETER_ARRAY"];
  control_parameter_array[1,1]=GAIN_MASTER;
  control_parameter_array[2,1]=GAIN_MASTER;
  control_parameter_array[3,1]=GAIN_MASTER;
  control_parameter_array[4,1]=GAIN_MASTER;
  control_parameter_array[5,1]=GAIN_MASTER;
  control_parameter_array[6,1]=GAIN_MASTER;
  control_parameter_array[7,1]=GAIN_MASTER;
  control_parameter_array[8,1]=GAIN_MASTER;
  control_parameter_array[9,1]=GAIN_MASTER;
  control_parameter_array[10,1]=GAIN_MASTER;
  control_parameter_array[11,1]=GAIN_MASTER;
  control_parameter_array[12,1]=GAIN_MASTER;
  control_parameter_array[13,1]=GAIN_MASTER;
  data_dictionary["CONTROL_PARAMETER_ARRAY"] = control_parameter_array;

  # Run the model to steady-state *before* we add the TGFB -
  XSS = EstimateSteadyState(data_dictionary);
  println(XSS)
  # ok, we established a steady state -
  # reset the initial conditions -
  initial_condition_array = vec(XSS);
  println(initial_condition_array)

  # Add TGFB1/2 -
  initial_condition_array[tgfb12_initial_cond_index] = 6.0; #	TGFB12
  initial_condition_array[R_initial_cond_index] = 6.0; #	Let there be receptors!
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Add LOW LEVEL VEGF -
  # initial_condition_array[vegf_initial_cond_index] = 1.0; #	VEGF

  # RESUME TGFB signaling for simulation
  rate_constant_array=data_dictionary["RATE_CONSTANT_ARRAY"];
  rate_constant_array[tgfb12_forward_rate_index]=4.0;   # forward
  rate_constant_array[tgfb12_reverse_rate_index]=0.1;  # Reverse
  rate_constant_array[R_deg_rate_index]=0.01; # attenuate degradation rate
  rate_constant_array[vegf_forward_rate_index]=0;#0;   # forward [] =([])=> VEGF
  rate_constant_array[vegf_reverse_rate_index]=0;#0;  # Reverse VEGF =([])=> [] (reverse)
  data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array;

  # Run the simulation -
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

# Simulate the addition of VEGF -
function SolveModelAddVEGF(TSTART,TSTOP,Ts,data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant

  # Modify rate constant array
  rate_constant_array=data_dictionary["RATE_CONSTANT_ARRAY"];
  # Kill TGFB AND VEGF signal for steady state
  rate_constant_array[tgfb12_forward_rate_index]=0;   # forward [] =([])=> TGFB12
  rate_constant_array[tgfb12_reverse_rate_index]=0;  # Reverse TGFB12 =([])=> [] (reverse)
  rate_constant_array[vegf_forward_rate_index]=0;   # forward [] =([])=> VEGF
  rate_constant_array[vegf_reverse_rate_index]=0;  # Reverse VEGF =([])=> [] (reverse)
  rate_constant_array[tgfb12_binding_rate_index]=3.0 # speed up TGFB12 binding
  rate_constant_array[tgfb3_binding_rate_index]=3.0 # speed up TGFB3 binding
  rate_constant_array[vegf_binding_rate_index]=3.0; # speed up VEGF binding
  rate_constant_array[PKC_activation_RAF_rate_index]=0.1 # bias towards NFATc1 (10x)
  rate_constant_array[bcat_ecad_forward_index]= bcat_complex_gain*ecad_bais_factor
  rate_constant_array[bcat_ecad_reverse_index]= 1.0
  rate_constant_array[bcat_tcf4_forward_index]= bcat_complex_gain
  rate_constant_array[bcat_tcf4_reverse_index]= 1.0
  rate_constant_array[bcat_lef1_forward_index]= bcat_complex_gain
  rate_constant_array[bcat_lef1_reverse_index]= 1.0
  rate_constant_array[bcat_deg_rate_index]= bcat_deg_rate
  rate_constant_array[RAS_activation_RAF_rate_index] = tgfb12_to_ap1_sp1_gain
  rate_constant_array[RAF_activation_MAPK_rate_index] = tgfb12_to_ap1_sp1_gain
  rate_constant_array[MAPK_activation_AP1_SP1_rate_index] = tgfb12_to_ap1_sp1_gain
  data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array;

  # Modify CONTROL_PARAMETER_ARRAY
  GAIN_MASTER = 0.5; # default is 0.1
  control_parameter_array=data_dictionary["CONTROL_PARAMETER_ARRAY"];
  control_parameter_array[1,1]=GAIN_MASTER;
  control_parameter_array[2,1]=GAIN_MASTER;
  control_parameter_array[3,1]=GAIN_MASTER;
  control_parameter_array[4,1]=GAIN_MASTER;
  control_parameter_array[5,1]=GAIN_MASTER;
  control_parameter_array[6,1]=GAIN_MASTER;
  control_parameter_array[7,1]=GAIN_MASTER;
  control_parameter_array[8,1]=GAIN_MASTER;
  control_parameter_array[9,1]=GAIN_MASTER;
  control_parameter_array[10,1]=GAIN_MASTER;
  control_parameter_array[11,1]=GAIN_MASTER;
  control_parameter_array[12,1]=GAIN_MASTER;
  control_parameter_array[13,1]=GAIN_MASTER;
  data_dictionary["CONTROL_PARAMETER_ARRAY"] = control_parameter_array;

  # Run the model to steady-state *before* we add the TGFB -
  XSS = EstimateSteadyState(data_dictionary);

  # ok, we established a steady state -
  # reset the initial conditions -
  initial_condition_array = vec(XSS);

  # Add VEGF -
  initial_condition_array[vegf_initial_cond_index] = 6.0; #	VEGF 6.0
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # RESUME VEGF signaling for simulation
  rate_constant_array=data_dictionary["RATE_CONSTANT_ARRAY"];
  rate_constant_array[vegf_forward_rate_index]=4.0;   # forward 4.0
  rate_constant_array[vegf_reverse_rate_index]=0.1;  # Reverse
  data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array;

  # Run the simulation -
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

# Run the model to steady state -
function EstimateSteadyState(data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    # DB 02/27/2016
    TSTART = 0.0
    Ts=1.0
    TSTOP=100
    # DB
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Kill ALL signals for steady state
  rate_constant_array=data_dictionary["RATE_CONSTANT_ARRAY"];
#  rate_constant_array[tgfb12_forward_rate_index]=0;   # forward
#  rate_constant_array[tgfb12_reverse_rate_index]=0;  # Reverse
#  rate_constant_array[vegf_forward_rate_index]=0;   # forward [] =([])=> VEGF
#  rate_constant_array[vegf_reverse_rate_index]=0;  # Reverse VEGF =([])=> [] (reverse)
#  rate_constant_array[TGFB12_R_activation_RAS_rate_index]=2.0 # speed up RAS activation due to Receptor delay
#  rate_constant_array[PKC_activation_RAF_rate_index]=0.1 # bias towards NFATc1 (10x)

  rate_constant_array[tgfb12_forward_rate_index]=0;   # forward [] =([])=> TGFB12
  rate_constant_array[tgfb12_reverse_rate_index]=0;  # Reverse TGFB12 =([])=> [] (reverse)
  rate_constant_array[vegf_forward_rate_index]=0;   # forward [] =([])=> VEGF
  rate_constant_array[vegf_reverse_rate_index]=0;  # Reverse VEGF =([])=> [] (reverse)
  rate_constant_array[tgfb12_binding_rate_index]=3.0 # speed up TGFB12 binding
  rate_constant_array[tgfb3_binding_rate_index]=3.0 # speed up TGFB3 binding
  rate_constant_array[vegf_binding_rate_index]=3.0; # speed up VEGF binding
  rate_constant_array[PKC_activation_RAF_rate_index]=0.1 # bias towards NFATc1 (10x)
  rate_constant_array[bcat_ecad_forward_index]= bcat_complex_gain*ecad_bais_factor
  rate_constant_array[bcat_ecad_reverse_index]= 1.0
  rate_constant_array[bcat_tcf4_forward_index]= bcat_complex_gain
  rate_constant_array[bcat_tcf4_reverse_index]= 1.0
  rate_constant_array[bcat_lef1_forward_index]= bcat_complex_gain
  rate_constant_array[bcat_lef1_reverse_index]= 1.0
  rate_constant_array[bcat_deg_rate_index]= bcat_deg_rate
  rate_constant_array[RAS_activation_RAF_rate_index] = tgfb12_to_ap1_sp1_gain
  rate_constant_array[RAF_activation_MAPK_rate_index] = tgfb12_to_ap1_sp1_gain
  rate_constant_array[MAPK_activation_AP1_SP1_rate_index] = tgfb12_to_ap1_sp1_gain
  data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array;

  # Setup loop -
  EPSILON = 1e-3;
  TSTART = 0.0;
  Ts = 1.0;
  TSTOP = 100;
  did_reach_steady_state = false
  while (!did_reach_steady_state)

    # solve the balances -
    (TSIM,X1) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

    # Take a few additional steps -
    TNEXT_START = TSTOP+Ts;
    TNEXT_STOP = TNEXT_START+1.0;
    Ts = 0.1;

    # solve the balances again 0
    initial_condition_array = vec(X1[end,:])
    data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

    (TSIM,X2) = SolveBalances(TNEXT_START,TNEXT_STOP,Ts,data_dictionary);

    # Find the difference -
    DIFF = norm((X2[end,:] - X1[end,:]));

    # Should we stop -or- go around again?
    if (DIFF<EPSILON)
      did_reach_steady_state = true;
      return (X2[end,:]);
    else

      # No, we did *not* reach steady state ....
      TSTART = TSTOP+Ts
      TSTOP = 1000 + TSTART;
      Ts = 100.0;

      initial_condition_array = vec(X2[end,:])
      data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
    end
  end
  # return
  return XSS;

end


text:
# samplepath.jl

"""
    abstract type AbstractSamplePath <: Any

subtype: SamplePath.

Its subtype implements a continuous samplepath.
"""
abstract type AbstractSamplePath end

"""
    isincreasing(x::AbstractVector{T})::Bool where T <: Number
    isincreasing(x::Range{T})::Bool where T <: Number

Tests whether a vector of numbers is strictly increasing. Is internal to
NonparametricBayesForDiffusions.
"""
function isincreasing(x::S) where S<:AbstractArray{T} where T <: Number
  for i in 1:length(x)-1
    x[i+1] <= x[i] && return false
  end
  return true
end

isincreasing(x::S) where S<:Range{T} where T <: Number = step(x) > 0

"""
    SamplePath(timeinterval::S, samplevalues::T) where
        {S<:AbstractVector{Float64}, T<:AbstractVector{Float64}}

 Implements a continuous samplepath where the timeinterval is not necessarily
 equally spaced. Sample value samplevalues[k] is the value of the process at
 time timeinterval[k]. timeinterval is an increasing vector. timeinterval and
 samplevalues are of equal length.

 Performance hint: Use Range objects instead of vectors as much as possible, this
 will speed up computation time.

# Constructors:

```julia
    SamplePath(timeinterval, samplevalues)
    SamplePath(timeinterval, f::Function) =
        SamplePath(timeinterval, f.(timeinterval))
```

# Examples

```julia
t = 0.0:0.01:2.0
X = SamplePath(t, sinpi)
```

```julia
t = 0.0:0.01:1.0
v = map(x->floor(10*x), t)
X = SamplePath(t, v)
```
"""
struct SamplePath{S, T} <: AbstractSamplePath where {S<:AbstractVector{Float64},
        T<:AbstractVector{Float64}}
  timeinterval::S
  samplevalues::T

  function SamplePath(timeinterval::S, samplevalues::T) where {
      S<:AbstractVector{Float64}, T<:AbstractVector{Float64}}
    length(timeinterval) == length(samplevalues) ||
      throw(DimensionMismatch("Length of timeinterval should be equal
       to the length of the samplevalues vector"))
    isincreasing(timeinterval) ||
      throw(ArgumentError("Timeinterval should be increasing"))
    new{S,T}(timeinterval, samplevalues)
  end
end

step(X::SamplePath{S}) where S<:Range{T} where T<:Number = step(X.timeinterval)

#constructor
SamplePath(timeinterval, f::Function) =
    SamplePath(timeinterval, f.(timeinterval))

# Returns the length of timeinterval == samplevalues vectors, not the endtime!
# TO DO: Answers to examples
"""
    length(X::SamplePath)

Returns the length of the vector timeinterval == length vector samplevalues.

# Examples
```julia
X = SamplePath([0.,1.,2.], [3.,5., -1.])
length(X)
```

```julia
X = SamplePath(0.0:0.1:2π, sin)
length(X) == length(0.0:0.1:2π)
```
"""
length(X::SamplePath)=length(X.timeinterval)


text:
# ### m13.3t.jl

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

df = CSV.read(sr_datadir("UCBadmit.csv"), DataFrame);

dept_map = Dict(key => idx for (idx, key) in enumerate(unique(df.dept)))
df.male = ifelse.(df.gender .== "male", 1, 0)
df.dept_id = [dept_map[de] for de in df.dept]

@model ppl13_3(applications, dept_id, male, admit) = begin
    Rho ~ LKJ(2, 2.)
    sigma_dept ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 2)
    bm ~ Normal(0, 1)
    a ~ Normal(0, 1)

    dist_mu = [a, bm]
    dist_Sigma = sigma_dept .* Rho .* sigma_dept'
    dist_Sigma = (dist_Sigma' + dist_Sigma) / 2

    a_bm_dept ~ filldist(MvNormal(dist_mu, dist_Sigma), 6)

    a_dept, bm_dept = a_bm_dept[1, :], a_bm_dept[2, :]
    logit_p = a_dept[dept_id] + bm_dept[dept_id] .* male

    admit .~ BinomialLogit.(applications, logit_p)
end

m13_3t = ppl13_3(df.applications, df.dept_id, df.male, df.admit)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
nchains = 1
chns13_3t = mapreduce(c -> sample(m13_3t, sampler, nsamples), chainscat, 1:nchains)

m_13_3_rethinking = """
Inference for Stan model: f0d86ec689cbf7921aab4fc0f55616d2.
    4 chains, each with iter=5000; warmup=1000; thin=1;
    post-warmup draws per chain=4000, total post-warmup draws=16000.

                      mean se_mean   sd     2.5%      25%      50%      75%
    bm_dept[1]       -0.79    0.00 0.27    -1.32    -0.97    -0.79    -0.61
    bm_dept[2]       -0.21    0.00 0.33    -0.87    -0.42    -0.20     0.00
    bm_dept[3]        0.08    0.00 0.14    -0.18    -0.01     0.08     0.18
    bm_dept[4]       -0.09    0.00 0.14    -0.37    -0.19    -0.09     0.00
    bm_dept[5]        0.12    0.00 0.19    -0.24    -0.01     0.12     0.25
    bm_dept[6]       -0.12    0.00 0.27    -0.67    -0.29    -0.12     0.05
    a_dept[1]         1.30    0.00 0.25     0.82     1.13     1.30     1.47
    a_dept[2]         0.74    0.00 0.33     0.08     0.52     0.73     0.95
    a_dept[3]        -0.65    0.00 0.09    -0.82    -0.71    -0.65    -0.59
    a_dept[4]        -0.62    0.00 0.10    -0.83    -0.69    -0.62    -0.55
    a_dept[5]        -1.13    0.00 0.11    -1.36    -1.21    -1.13    -1.05
    a_dept[6]        -2.60    0.00 0.20    -3.01    -2.73    -2.60    -2.46
    a                -0.49    0.01 0.73    -1.96    -0.91    -0.49    -0.07
    bm               -0.16    0.00 0.24    -0.65    -0.29    -0.16    -0.02
    sigma_dept[1]     1.67    0.01 0.63     0.88     1.25     1.53     1.94
    sigma_dept[2]     0.50    0.00 0.26     0.16     0.33     0.45     0.61
    Rho[1,1]          1.00     NaN 0.00     1.00     1.00     1.00     1.00
    Rho[1,2]         -0.31    0.00 0.35    -0.87    -0.59    -0.35    -0.08
    Rho[2,1]         -0.31    0.00 0.35    -0.87    -0.59    -0.35    -0.08
    Rho[2,2]          1.00    0.00 0.00     1.00     1.00     1.00     1.00
    lp__          -2593.92    0.04 3.17 -2601.03 -2595.85 -2593.57 -2591.65
                     97.5% n_eff Rhat
    bm_dept[1]       -0.27  6918    1
    bm_dept[2]        0.46  9133    1
    bm_dept[3]        0.36 13016    1
    bm_dept[4]        0.19 12258    1
    bm_dept[5]        0.49 12953    1
    bm_dept[6]        0.40 12011    1
    a_dept[1]         1.81  6993    1
    a_dept[2]         1.39  9335    1
    a_dept[3]        -0.48 13489    1
    a_dept[4]        -0.42 12595    1
    a_dept[5]        -0.91 14659    1
    a_dept[6]        -2.23 13168    1
    a                 0.96  9541    1
    bm                0.31  8946    1
    sigma_dept[1]     3.27  9058    1
    sigma_dept[2]     1.15  6626    1
    Rho[1,1]          1.00   NaN  NaN
    Rho[1,2]          0.44 11358    1
    Rho[2,1]          0.44 11358    1
    Rho[2,2]          1.00 16069    1
    lp__          -2588.70  5310    1
"""

# End of m13.3.jl


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

name = "Catalyst"
version = v"0.1.0"
sources = [
    GitSource("https://gitlab.kitware.com/paraview/catalyst",
	      "16dc369855a7c29bb00829cf8ecb16fa3b7ebd4b")
]

script = raw"""
cd ${WORKSPACE}/srcdir/catalyst*
rm -rf build && mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      -DCATALYST_BUILD_SHARED_LIBS=ON \
      -DCATALYST_BUILD_STUB_IMPLEMENTATION=ON \
      -DCATALYST_BUILD_TESTING=OFF \
      -DCATALYST_BUILD_TOOLS=ON \
      -DCATALYST_USE_MPI=OFF \
      ..
make -j${nproc}
make install
# install catalyst-replay to bindir
mkdir -vp $bindir
cp -v ./bin/catalyst_replay* $bindir
"""

# Paraview / Catalyst only supports x86_64 builds for Linux, Windows
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
]
platforms = expand_cxxstring_abis(platforms)

# The stub impl has the wrong extension .so for darwin builds, need to fix upstream
products = [
    LibraryProduct("libcatalyst", :libcatalyst),
    ExecutableProduct("catalyst_replay", :catalyst_replay_exe),
    # LibraryProduct("libcatalyst-stub", Symbol("libcatalyst-stub"), ["lib/catalyst"]),
]

dependencies = Dependency[
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
	       julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil