text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.453305806680525 <= q <= 37.480198521751326)

                     add_NL_constraint(m, :(sigmoid(-0.5060529352779994 + 0.34605798298901425*sigmoid(-0.49627455608242954 + -0.46375151869975495*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.8504066842806357*sigmoid(-0.2321501407718678 + 0.5600091754803458*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.462211279142271*sigmoid(-0.9807289687619161 + -0.5968595287514287*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.0428411186394575*sigmoid(0.7899732054652109 + -0.5363355140855166*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + sigmoid(0.02992678336456356 + -0.8781954986928167*sigmoid(-0.49627455608242954 + -0.46375151869975495*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.84297522225665*sigmoid(-0.2321501407718678 + 0.5600091754803458*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.44937020244474857*sigmoid(-0.9807289687619161 + -0.5968595287514287*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.012040065223148932*sigmoid(0.7899732054652109 + -0.5363355140855166*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + sigmoid(-0.25761034404046734 + -0.19911102915318102*sigmoid(-0.49627455608242954 + -0.46375151869975495*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.7457041629864394*sigmoid(-0.2321501407718678 + 0.5600091754803458*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.7638178262546895*sigmoid(-0.9807289687619161 + -0.5968595287514287*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.07889233901952375*sigmoid(0.7899732054652109 + -0.5363355140855166*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + sigmoid(0.1153458350428731 + 0.9297971624965538*sigmoid(-0.49627455608242954 + -0.46375151869975495*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.7017891371439648*sigmoid(-0.2321501407718678 + 0.5600091754803458*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.8093911002605614*sigmoid(-0.9807289687619161 + -0.5968595287514287*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.43763294749909276*sigmoid(0.7899732054652109 + -0.5363355140855166*sigmoid(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*sigmoid(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*sigmoid(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*sigmoid(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module FluxTraining

#=
Refactoring:
- better serialization of Learners
=#

using LightGraphs
using BSON: @load, @save
using Flux
using Flux: Params, onecold
using Flux.Optimise: update!
using Glob
using DocStringExtensions
using LearnBase: nobs
import OnlineStats
using OnlineStats: EqualWeight, Mean, OnlineStat
using Parameters
using ProgressMeter: Progress, next!
using Statistics: mean
using UUIDs
using Zygote
using Zygote: Grads, gradient
using ValueHistories
using DataStructures: DefaultDict

# optimizers
include("./optimizers.jl")

# functional
include("./functional/anneal.jl")
include("./functional/metrics.jl")

# callback framework
include("./callbacks/protect.jl")
include("./callbacks/phases.jl")
include("./callbacks/events.jl")
include("./callbacks/callback.jl")
include("./callbacks/graph.jl")
include("./callbacks/execution.jl")

# utilities
include("./util/ioutils.jl")

# callback implementations
include("./callbacks/callbacks.jl")
include("./callbacks/customcallback.jl")
include("./callbacks/metrics.jl")
include("./callbacks/recorder.jl")
include("./callbacks/paramscheduler.jl")
include("./callbacks/checkpointer.jl")

# learner
include("./learner.jl")

include("./train.jl")
include("./util/trainutils.jl")

export AbstractCallback,
    AbstractMetric,
    Accuracy,
    AverageLoss,
    CancelBatchException,
    CancelEpochException,
    CancelFittingException,
    Checkpointer,
    CheckpointAny,
    CheckpointLowest,
    CustomCallback,
    DataBunch,
    DataLoader,
    EarlyStopping,
    ToGPU,
    GarbageCollect,
    Learner,
    LRFinderPhase,
    Metric,
    MeanMetric,
    Recorder,
    ProgressBarLogger,
    MetricsLogger,
    ParamSchedule,
    ParamScheduler,
    TrainingPhase,
    ValidationPhase,
    Schedule,
    Schedules,
    StopOnNaNLoss,

    accuracy,
    anneal_linear,
    anneal_cosine,
    anneal_exp,
    anneal_const,
    endtraining,
    fit!,
    loadmodel,
    loadweights,
    onecycleschedule,
    savemodel,
    saveweights,
    setschedule!,
    splitdataset,
    starttraining

end # module

ext: jl


text:
@testset "LinearBucket.jl" begin
    @testset "Constructors" begin
		model = LinearBucket()
		model2 = LinearBucket(model)
		@test model == model2
	end

    @testset "Check param or state values when creating new instance" begin
        @test_throws DomainError LinearBucket(-0.1, 0.0)
        @test_throws DomainError LinearBucket(1.1, 0.0)
        @test_throws DomainError LinearBucket(0.1, -0.5)
    end

    @testset "Check values when a parameter or state value is set" begin
        model = LinearBucket(0.1, 0.0)
        @test_throws DomainError model.K = -0.1
        @test_throws DomainError model.K = 1.1
        @test_throws DomainError model.S = -0.5
    end

    @testset "Stash and recover a model for hot start" begin
        model = LinearBucket(0.2, 0.5)
        d = stash(model)
        model2 = LinearBucket(d)
        @test model == model2
    end

    @testset "Print a model" begin
    	model = LinearBucket()
    	println(model)
	end
end


text:
using Test

include("dnd-character.jl")

function ischaracter(c)
    c.strength >= 3         && c.strength <= 18 &&
        c.dexterity >= 3    && c.dexterity <= 18 &&
        c.constitution >= 3 && c.constitution <= 18 &&
        c.intelligence >= 3 && c.intelligence <= 18 &&
        c.wisdom >= 3       && c.wisdom <= 18 &&
        c.charisma >= 3     && c.charisma <= 18 &&
        c.hitpoints == 10 + modifier(c.constitution)
end

@testset "ability modifier" begin
    @test modifier(3) == -4
    @test modifier(4) == -3
    @test modifier(5) == -3
    @test modifier(6) == -2
    @test modifier(7) == -2
    @test modifier(8) == -1
    @test modifier(9) == -1
    @test modifier(10) == 0
    @test modifier(11) == 0
    @test modifier(12) == 1
    @test modifier(13) == 1
    @test modifier(14) == 2
    @test modifier(15) == 2
    @test modifier(16) == 3
    @test modifier(17) == 3
    @test modifier(18) == 4
end

@testset "random ability is within range" begin
    abil = ability()
    @test abil >= 3 && abil <= 18
end

@testset "random character is valid" begin
    characters = DNDCharacter[]

    for i=1:10
        c = DNDCharacter()
        @testset "random character is valid and unique in history" begin
            @test ischaracter(c)
            @test !in(c, characters)
        end

        push!(characters, c)
    end
end


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

mutable struct ConnectionResetSharedKey <: SwaggerModel
    keyLength::Any # spec type: Union{ Nothing, Int32 } # spec name: keyLength

    function ConnectionResetSharedKey(;keyLength=nothing)
        o = new()
        validate_property(ConnectionResetSharedKey, Symbol("keyLength"), keyLength)
        setfield!(o, Symbol("keyLength"), keyLength)
        o
    end
end # type ConnectionResetSharedKey

const _property_map_ConnectionResetSharedKey = Dict{Symbol,Symbol}(Symbol("keyLength")=>Symbol("keyLength"))
const _property_types_ConnectionResetSharedKey = Dict{Symbol,String}(Symbol("keyLength")=>"Int32")
Base.propertynames(::Type{ ConnectionResetSharedKey }) = collect(keys(_property_map_ConnectionResetSharedKey))
Swagger.property_type(::Type{ ConnectionResetSharedKey }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionResetSharedKey[name]))}
Swagger.field_name(::Type{ ConnectionResetSharedKey }, property_name::Symbol) =  _property_map_ConnectionResetSharedKey[property_name]

function check_required(o::ConnectionResetSharedKey)
    (getproperty(o, Symbol("keyLength")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ConnectionResetSharedKey }, name::Symbol, val)
    if name === Symbol("keyLength")
        Swagger.validate_param(name, "ConnectionResetSharedKey", :maximum, val, 128, false)
        Swagger.validate_param(name, "ConnectionResetSharedKey", :minimum, val, 1, false)
    end
end


text:
#From http://home.online.no/~pjacklam/notes/invnorm/

function CTΦInv(p::Float64)::Float64

    a1::Float64 = -39.6968302866538
    a2::Float64 = 220.946098424521
    a3::Float64 = -275.928510446969
    a4::Float64 = 138.357751867269
    a5::Float64 = -30.6647980661472
    a6::Float64 = 2.50662827745924

    b1::Float64 = -54.4760987982241
    b2::Float64 = 161.585836858041
    b3::Float64 = -155.698979859887
    b4::Float64 = 66.8013118877197
    b5::Float64 = -13.2806815528857

    c1::Float64 = -7.78489400243029E-03
    c2::Float64 = -0.322396458041136
    c3::Float64 = -2.40075827716184
    c4::Float64 = -2.54973253934373
    c5::Float64 = 4.37466414146497
    c6::Float64 = 2.93816398269878

    d1::Float64 = 7.78469570904146E-03
    d2::Float64 = 0.32246712907004
    d3::Float64 = 2.445134137143
    d4::Float64 = 3.75440866190742

    #Rational approximation for lower region.

    if p < 0.02425
        q::Float64 = √(-2.0 * log(p))
        ans::Float64 = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1.0)
    elseif p < 0.97575
        q = (p - 0.5) * (p - 0.5)
        ans = (((((a1 * q + a2) * q + a3) * q + a4) * q + a5) * q + a6) * (p - 0.5) / (((((b1 * q + b2) * q + b3) * q + b4) * q + b5) * q + 1)
    else
        q = √(-2 * log(1 - p))
        ans = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
    end

    return ans
end

CTΦ(z::Real)::Real = 0.5+0.5*erf(z*sqrt2Inv)
CTϕ(z::Real)::Real = sqrt2PiInv*exp(-z^2.0/2.0)

#in-place binary search, returns only a boolean value
function inSorted(val::T, B::Vector{T})::Bool where T<:Any

  searched::Bool = false
  low::Int = 1
  high::Int = length(B)
  pivot::T = B[(high+low) ÷ 2]
  found::Bool = false

  #in-place binary sort
  while !found && !searched

    #println("high: $high low: $low mid: $((high+low)÷2) pivot: $pivot")

    if val > pivot
      low = (high+low) ÷ 2 + 1
    elseif val < pivot
      high = (high+low) ÷ 2 - 1
    else
      found = true
    end

    if high-low ≤ 1
      searched = true
      found = val == B[low] || val == B[high]
    end
    pivot = B[(high+low) ÷ 2]
  end

  return found
end

#=function testInSorted(;N::Int = 2000, K::Int = 10_000)

  for i ∈ 1:K
    A::Vector{Int} = rand(1:N, N)
    B::Vector{Int} = rand(1:N, N)

    sort!(B)
    for j ∈ 1:length(A)
      if inSorted(A[j],B) ≠ (A[j] ∈ B)
        error("inSorted:$(inSorted(A[j],B)) while A∈B:$(A[j] ∈ B)\nA: $(A[j])\nB:$B")
      end
    end
  end
end

@time testInSorted()=#

#returns a list of all elements in A
#compare in performance to setdiff. Maintains order unless sort=true
function ANotB(A::Vector{T}, B::Vector{T};
  sort::Bool=false, sorted::Bool=false) where T<:Any

  if !sort && !sorted
    A = deepcopy(A)
    B = deepcopy(B)
  end

  #Force mergesort to deal with sorted cases
  if !sorted
    if !issorted(A)
      sort!(A)
    end
    if !issorted(B)
      sort!(B)
    end
  end

  numA::Int = length(A)
  numB::Int = length(B)

  findings::Vector{Bool} = Vector{Bool}(numA)

  @fastmath @inbounds @simd for i ∈ 1:numA
    findings[i] = length(searchsorted(B,A[i]))>0
  end

  return A[findings]

end

function ANotBTest(N::Int = 1_000_000)
  println("\nStarting test of algorithm N=$N")
  A::Vector{Symbol} = (Symbol).(rand(1:N, N))
  B::Vector{Symbol} = (Symbol).(rand(1:N, N))

  @time setdiff(A,B) #for comparison
  @time ANotB(A,B, sort=true)
end

max_stars_repo_path: Fall 17 Ivo Project 2.2/histtaxes/CTCode/CTNumerical.jl

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

###############
# Prepare FMU #
###############

pathToFMU = joinpath(dirname(@__FILE__), "../model/IO.fmu")

myFMU = fmiLoad(pathToFMU)
@test fmiInstantiate!(myFMU; loggingOn=false) != 0

@test fmiSetupExperiment(myFMU, 0.0) == 0

@test fmiEnterInitializationMode(myFMU) == 0

#########################
# Testing Setter/Getter #
#########################

rndReal = 100 * rand()
rndInteger = round(Integer, 100 * rand())
rndBoolean = rand() > 0.5
rndString = Random.randstring(12)

@test fmiSetReal(myFMU, "p_real", rndReal) == 0
@test fmiGetReal(myFMU, "p_real") == rndReal
@test fmiSetReal(myFMU, "p_real", -rndReal) == 0
@test fmiGetReal(myFMU, "p_real") == -rndReal

@test fmiSetInteger(myFMU, "p_integer", rndInteger) == 0
@test fmiGetInteger(myFMU, "p_integer") == rndInteger
@test fmiSetInteger(myFMU, "p_integer", -rndInteger) == 0
@test fmiGetInteger(myFMU, "p_integer") == -rndInteger

@test fmiSetBoolean(myFMU, "p_boolean", rndBoolean) == 0
@test fmiGetBoolean(myFMU, "p_boolean") == rndBoolean
@test fmiSetBoolean(myFMU, "p_boolean", !rndBoolean) == 0
@test fmiGetBoolean(myFMU, "p_boolean") == !rndBoolean

@test fmiSetString(myFMU, "p_string", rndString) == 0
@test fmiGetString(myFMU, "p_string") == rndString

############
# Clean up #
############

@test fmiExitInitializationMode(myFMU) == 0

fmiUnload(myFMU)


text:
#fin. gen. submodules of Z^n and F_p^n (and possibly more)
#
#import Base.show, Base.inv, Nemo.inv, Base.solve, Hecke.solve,
#       Hecke.lift, Hecke.rational_reconstruction, Hecke.elementary_divisors,
#       Hecke.rank, Hecke.det

export det_mc, id, isupper_triangular, norm2, hadamard_bound2,
       hnf, hnf!, echelon_with_transform

add_verbose_scope(:HNF)

add_assert_scope(:HNF)
set_assert_level(:HNF, 0)

function show(io::IO, M::ModuleCtxNmod)
  print(io, "Sparse module over $(M.R) of (current) rank $(nrows(M.basis)) and $(nrows(M.gens))\n")
end

function show(io::IO, M::ModuleCtx_fmpz)
  print(io, "Sparse module over FlintZZ of (current) rank $(nrows(M.bas_gens)) and further $(nrows(M.rel_gens))\n")
  if isdefined(M, :basis_idx)
    print(io, "current index: $(M.basis_idx)\n")
  end
  if isdefined(M, :essential_elementary_divisors)
    print(io, "current structure: $(M.essential_elementary_divisors)\n")
  end
end

function add_gen!(M::ModuleCtxNmod, g::SRow{nmod})
  if M.basis.r == M.basis.c
    return false
  end
  h = reduce(M.basis, g)
  if !iszero(h)
    i = 1
    while i<= nrows(M.basis) && M.basis.rows[i].pos[1] < h.pos[1]
      i += 1
    end
    @hassert :HNF 2  i > nrows(M.basis) || M.basis[i].pos[1] > h.pos[1]
    insert!(M.basis.rows, i, h)
    M.basis.r += 1
    M.basis.nnz += length(h)
    M.basis.c = max(M.basis.c, h.pos[end])
    push!(M.gens, g)
    return true
  end
  return false
end

function add_gen!(M::ModuleCtx_fmpz, g::SRow{fmpz}, always::Bool = true)
  gp = change_base_ring(M.Mp.R, g)
  M.new = true
  if add_gen!(M.Mp, gp)
    push!(M.bas_gens, g)
    return true
  else
    always && push!(M.rel_gens, g)
  end
  return false
end

function check_index(M::ModuleCtx_fmpz)
  if nrows(M.Mp.basis) < ncols(M.Mp.basis)
    return fmpz(0)
  end

  if !M.new
    return M.basis_idx
  end

  if isdefined(M, :trafo)  #once we have trafo, we need to keep it!
    module_trafo_assure(M)
    return M.basis_idx
  end

  M.new = false

  if isdefined(M, :basis)
    C = copy(M.basis)
    @assert isupper_triangular(C)
    @assert M.basis_idx != 0
  else
    d = abs(det_mc(M.bas_gens))
    C = M.max_indep
    C.c = M.bas_gens.c
    for ii = M.bas_gens
      h = reduce(C, ii, 2*d) #to avoid problems with diag being 1...1 d
      @hassert :HNF 2  !iszero(h)
      i = 1
      while i<= nrows(C) && C.rows[i].pos[1] < h.pos[1]
        i += 1
      end
      @hassert :HNF 2  i > nrows(C) || C[i].pos[1] > h.pos[1]
      insert!(C.rows, i, h)
      C.r += 1
      C.nnz += length(h)
      C.c = max(C.c, h.pos[end])
    end
    M.max_indep = copy(C)
    @assert isupper_triangular(C)
    M.basis_idx = prod([C[i,i] for i=1:nrows(C)])
  end

  d = 2*M.basis_idx

  for i=length(M.rel_reps_p)+1:length(M.rel_gens)
    reduce(C, M.rel_gens[i], d)
    push!(M.rel_reps_p, solve_ut(M.Mp.basis, change_base_ring(M.Mp.R, M.rel_gens[i])))
  end

#=
  for l=1:5
    mis = find(i->C[i,i] != 1, 1:nrows(C))
    if length(mis) == 0
      break
    end
#    println("mis: $mis")
    for i = mis
      if C[i,i] == 1
#        println("miracle for $i")
        continue
      end
      r = find(x->i in M.rel_reps_p[x].pos, 1:length(M.rel_reps_p))
#      println("found $(length(r)) rows")
      if length(r) == 0
        break
      end
      g = M.rel_gens[rand(r)]
      for j=1:min(5, div(length(r), 2))
        g += rand(-10:10)*M.rel_gens[rand(r)]
      end
      reduce(C, g)
      if C[i,i] == 1
#        println("bingo for i=$i")
      end
    end
  end
=#
  M.basis = C
  M.basis_idx = prod([C[i,i] for i=1:nrows(C)])

  return M.basis_idx
end

function elementary_divisors(M::ModuleCtx_fmpz)
  if !isdefined(M, :basis)
    i = check_index(M)
    if i == 0
      return fmpz[]
    end
  end
  C = M.basis
  f = find(i -> C[i,i] != 1, 1:nrows(C))
  if length(f) == 0
    M.essential_elementary_divisors = fmpz[]
    return M.essential_elementary_divisors
  end
  e = minimum(f)
  m = fmpz_mat(sub(C, e:nrows(C), e:ncols(C)))
  s = snf(m)
  M.essential_elementary_divisors = [s[i,i] for i=1:nrows(s)]
  return M.essential_elementary_divisors
end

function missing_pivot(M::ModuleCtx_fmpz)
  C = M.Mp.basis
  return setdiff(BitSet(1:ncols(C)), [x.pos[1] for x=C])
end

function non_trivial_pivot(M::ModuleCtx_fmpz)
  h = check_index(M)
  if h == 0
    return missing_pivot(M)
  end
  C = M.basis
  @hassert :HNF 2  C.r == C.c
  return setdiff(BitSet(1:ncols(C)), findall(i->C[i].values[1] == 1, 1:C.c))
end

function rank(M::ModuleCtx_fmpz)
  return M.bas_gens.r
end

function rank(M::ModuleCtxNmod)
  return M.basis.r
end

function module_trafo_assure(M::ModuleCtx_fmpz)

  if !M.new && isdefined(M, :trafo)
    return nothing
  end
  if isdefined(M, :trafo)
    st = M.done_up_to + 1
    _, t = hnf_extend!(M.basis, sub(M.rel_gens, st:nrows(M.rel_gens), 1:ncols(M.rel_gens)), Val{true}, offset = st-1, truncate = true)
    append!(M.trafo, t)

  else
    z = vcat(M.bas_gens, M.rel_gens)
    h, t = hnf_kannan_bachem(z, Val{true}, truncate = true)
    M.trafo = t
    M.basis = h
  end
  M.done_up_to = nrows(M.rel_gens)
  M.basis_idx = det(M.basis) # h is upp_triangular, hence det is trivial
  @assert M.basis_idx > 0
  M.new = false

  return nothing
end


text:
include("../cosine-map.jl")

using Gadfly
using Compose

p = plot(x=[1:k], y=x,
    Geom.point,
    Geom.line,

    Guide.Title("Cosine Map"),
    Guide.YLabel("f(x)"),
    Guide.XLabel("Iteration"),
    Guide.XTicks(ticks=[1:2:k]),
    Guide.Annotation(
      compose(context(),
        compose(context(), text(fp, x[fp] + 0.04, "Fixed Point", hleft, vcenter)),
        compose(context(), circle(fp, x[fp], 2.5mm), fill(nothing), stroke("orange"))
      )
    ),
    Guide.Annotation(
      compose(context(),
        compose(context(),
          line([(1.4, x[fp] + 0.02), (fp - 0.4, x[fp] + 0.02)]),
          fill(nothing),
          stroke("steelblue"),
          strokedash([2mm, 1mm])
        ),
        compose(context(), text(fp / 2, x[fp] + 0.05, "Transient", hcenter, vcenter))
      )
    ))

draw(PNG("cosine-map.png", 600px, 600px / golden), p)


text:
# Rebuild shared library
try
    using CSDP
    isfile(CSDP.csdp) && rm(CSDP.csdp)
catch
    println("Lib unloadable")
end

Pkg.build("CSDP")


text:
module TestCore

using Atomix
using Atomix.Internal: referenceable
using Test

function test_indexableref()
    A = ones(Int, 3)
    ref = Atomix.IndexableRef(A, (1,))
    @test eltype(ref) === Int
    @test Atomix.get(ref) === 1
    Atomix.set!(ref, 123)
    @test Atomix.get(ref) === 123
    @test Atomix.modify!(ref, +, 1) === (123 => 124)
    @test Atomix.get(ref) === 124
    @test Atomix.swap!(ref, 345) == 124
    @test Atomix.get(ref) === 345
    @test Atomix.replace!(ref, 345, 567) === (old = 345, success = true)
    @test Atomix.replace!(ref, 345, 567) === (old = 567, success = false)
end

function test_referenceablearray()
    @testset for a in Any[
        ones(Int, 3),
        view(ones(Int, 3), 1:2),
        view(ones(Int, 2, 3), 1:1, 1:2),
        ones(Int, 3)',
    ]
        ra = referenceable(a)
        @test size(ra) == size(a)
        @test IndexStyle(ra) == IndexStyle(a)
    end
end

end  # module


text:
function randdom_search(objective::Function, search_space::Tuple{Real,Real};
    max_iter::Int=100, op = <)
    rnd() = search_space[1] + (search_space[2] - search_space[1]) * rand()
    best = objective(rnd())
    for i = 2:max_iter
        candidate = objective(rnd())
        op(candidate, best) && (best = candidate)
    end
    best
end


text:
# Testing gammatonegram response
using Compat.DelimitedFiles
s = read!(open(joinpath(dirname(@__FILE__), "data", "sa2.raw")), Array{Float64}(undef,42701))
fs = 16000
(W, _) = AuditoryFilters.fft2gammatonemx(512, fs, 23, 1, 150, 4000)
G = gammatonegram(s, fs, 0.025, 0.010, 23, 150, 4000, 1)

G_matlab = readdlm(open(joinpath(dirname(@__FILE__), "data", "G_matlab.csv")), ',')
F_matlab = readdlm(open(joinpath(dirname(@__FILE__), "data", "F_matlab.csv")), ',')
W_matlab = readdlm(open(joinpath(dirname(@__FILE__), "data", "W_matlab.csv")), ',')

@test W ≈ W_matlab
@test G.frequencies ≈ F_matlab
@test G.amplitude ≈ G_matlab


text:
module VideoIO

using Libdl
using Requires, Dates, ProgressMeter
using ImageCore: channelview, rawview
using ColorTypes: RGB, Gray, N0f8, N6f10, YCbCr, Normed, red, green, blue
using FileIO: File

using Base: fieldindex, RefValue, sigatomic_begin, sigatomic_end, cconvert
using Base.GC: @preserve
import Base:
    iterate,
    IteratorSize,
    IteratorEltype,
    setproperty!,
    convert,
    getproperty,
    unsafe_convert,
    propertynames,
    getindex,
    setindex!,
    parent,
    position,
    unsafe_wrap,
    unsafe_copyto!,
    write

const VIO_LOCK = ReentrantLock()

include("util.jl")
include("../lib/libffmpeg.jl")
using .libffmpeg
using FFMPEG: ffmpeg

include("avptr.jl")

include("info.jl")
include("avdictionary.jl")
include("avframe_transfer.jl")
include("frame_graph.jl")
include("avio.jl")
include("encoding.jl")
include("testvideos.jl")
using .TestVideos

if Sys.islinux()
    using Glob: Glob
    function init_camera_devices()
        append!(CAMERA_DEVICES, Glob.glob("video*", "/dev"))
        return DEFAULT_CAMERA_FORMAT[] = libffmpeg.av_find_input_format("video4linux2")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        return DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "" : CAMERA_DEVICES[1]
    end
elseif Sys.iswindows()
    function init_camera_devices()
        append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "dshow", "dummy"))
        return DEFAULT_CAMERA_FORMAT[] = libffmpeg.av_find_input_format("dshow")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        return DEFAULT_CAMERA_DEVICE[] = string("video=", isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1])
    end
elseif Sys.isapple()
    function init_camera_devices()
        try
            append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "avfoundation", "\"\""))
            DEFAULT_CAMERA_FORMAT[] = libffmpeg.av_find_input_format("avfoundation")
        catch
            try
                append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "qtkit", "\"\""))
                DEFAULT_CAMERA_FORMAT[] = libffmpeg.av_find_input_format("qtkit")
            catch
            end
        end
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        # Note: "Integrated" is another possible default value
        DEFAULT_CAMERA_OPTIONS["pixel_format"] = "uyvy422"
        return DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1]
    end
elseif Sys.isbsd()
    # copied loosely from apple above - needs figuring out
    function init_camera_devices()
        append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "avfoundation", "\"\""))
        return DEFAULT_CAMERA_FORMAT[] = libffmpeg.av_find_input_format("avfoundation")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        DEFAULT_CAMERA_OPTIONS["pixel_format"] = "uyvy422"
        return DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1]
    end
end

#Helper functions to explain about GLMakie load order requirement
function play(f; flipx = false, flipy = false)
    return error(
        "GLMakie must be loaded before VideoIO to provide video playback functionality. Try a new session with `using GLMakie, VideoIO`",
    )
end
function playvideo(video; flipx = false, flipy = false)
    return error(
        "GLMakie must be loaded before VideoIO to provide video playback functionality. Try a new session with `using GLMakie, VideoIO`",
    )
end
function viewcam(device = DEFAULT_CAMERA_DEVICE, format = DEFAULT_CAMERA_FORMAT)
    return error(
        "GLMakie must be loaded before VideoIO to provide camera playback functionality. Try a new session with `using GLMakie, VideoIO`",
    )
end

## FileIO interface
fileio_load(f::File; kwargs...) = load(f.filename; kwargs...)
fileio_save(f::File, video; kwargs...) = save(f.filename, video; kwargs...)

function __init__()
    # Always check your dependencies from `deps.jl`
    # TODO remove uncessary ENV["LD_LIBRARY_PATH"] from check_deps, so that
    # it doesn't mess with LD_LIBRARY_PATH, which was causing CI download issues due to issues with julia's curl
    # since check_deps is optional, I hope this is ok for now

    #check_deps()

    loglevel!(libffmpeg.AV_LOG_FATAL)
    # @info "VideoIO: Low-level FFMPEG reporting set to minimal (AV_LOG_FATAL). See `? VideoIO.loglevel!` for options"

    read_packet[] = @cfunction(_read_packet, Cint, (Ptr{AVInput}, Ptr{UInt8}, Cint))

    av_register_all()

    libffmpeg.avdevice_register_all()

    @require GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin
        # Define read and retrieve for Images
        function play(f; flipx = false, flipy = false, pixelaspectratio = nothing)
            eof(f) && error("VideoReader at end of file. Use `seekstart(f)` to rewind")
            # if user did not specify the aspect ratio we'll try to use the one stored in the video file
            if pixelaspectratio === nothing
                pixelaspectratio = aspect_ratio(f)
            end
            h = height(f)
            w = round(typeof(h), width(f) * pixelaspectratio) # has to be an integer
            flips_to_dims = Dict(
                (true, true) => (1, 2),
                (true, false) => 1,
                (false, true) => 2,
                (false, false) => nothing,
            )
            flipping_dims = flips_to_dims[(flipx, flipy)]
            flipping = i -> i
            if flipping_dims !== nothing
                flipping = i -> reverse(i, dims = flipping_dims)
            end
            flip_and_rotate = i -> begin
                rotated = GLMakie.rotr90(i)
                flipping(rotated)
            end
            img = read(f)
            obs_img = GLMakie.Observable(flip_and_rotate(img))
            scene =
                GLMakie.Scene(camera = GLMakie.campixel!, resolution = reverse(size(img)))

            GLMakie.image!(scene, obs_img)
            display(scene)
            # issue 343: camera can't run at full speed on MacOS
            fps = Sys.isapple() ? min(framerate(f), 24) : framerate(f)
            while isopen(scene) && !eof(f)
                read!(f, img)
                obs_img[] = flip_and_rotate(img)
                sleep(1 / fps)
            end
        end

        function playvideo(video; flipx = false, flipy = false, pixelaspectratio = nothing)
            f = VideoIO.openvideo(video)
            try
                play(f, flipx=flipx, flipy=flipy, pixelaspectratio=pixelaspectratio)
            finally
                close(f)
            end
        end

        function viewcam(device=nothing, format=nothing, options=nothing, pixelaspectratio=nothing)
            camera = opencamera(device, format, options)
            try
                play(camera; flipx=true, pixelaspectratio)
            finally
                close(camera)
            end
        end
    end
end

"""
VideoIO supports reading and writing video files.

  - `VideoIO.load` to load an entire video into memory as a vector of images (a framestack)
  - `openvideo` and `opencamera` provide access to video files and livestreams
  - `read` and `read!` allow reading frames
  - `seek`, `seekstart`, `skipframe`, and `skipframes` support access of specific frames
  - `VideoIO.save` for encoding an entire framestack in one step
  - `open_video_out`, `write` for writing frames sequentially to a file
  - `gettime` and `counttotalframes` provide information

Here's a brief demo reading through each frame of a video:

```julia
using VideoIO
r = openvideo(filename)
img = read(r)
while !eof(r)
    read!(r, img)
end
```

An example of encoding one frame at a time:

```julia
using VideoIO
framestack = map(x -> rand(UInt8, 100, 100), 1:100) #vector of 2D arrays
encoder_options = (crf = 23, preset = "medium")
open_video_out("video.mp4", framestack[1], framerate = 24, encoder_options = encoder_options) do writer
    for frame in framestack
        write(writer, frame)
    end
end
```
"""
VideoIO

end # VideoIO


text:
import DandelionWebSockets: AbstractWSClient, stop, send_text, send_binary

type MockThrottlingWSClient <: AbstractWSClient
    sent::Vector{UTF8String}
    bin_send::Vector{Vector{UInt8}}
    closed_called::Int

    function MockThrottlingWSClient()
        new(Vector{UTF8String}(), Vector{Vector{UInt8}}(), 0)
    end
end

function stop(c::MockThrottlingWSClient)
    c.closed_called += 1
end

send_text(c::MockThrottlingWSClient, s::UTF8String) = push!(c.sent, s)

function expect(c::MockThrottlingWSClient, s::UTF8String)
    if isempty(c.sent)
        error("Expected sent to have at least one element, with text '$s'")
    end
    actual = shift!(c.sent)
    @fact actual --> s
end

facts("Throttling") do
    context("ThrottledWSClient") do
        throttle_time = 0.05
        mock_ws = MockThrottlingWSClient()
        throttled_ws = ThrottledWSClient(mock_ws, throttle_time)

        n = 5
        for i in 1:n
            send_text(throttled_ws, utf8("Hello"))
        end

        # Sleep a short amount of time
        sleep(0.1)
        # All events should not have reached the mock WebSocket client.
        @fact length(mock_ws.sent) < n --> true

        # After this sleep all events should have reached the WebSocket client
        sleep(throttle_time * n)
        @fact length(mock_ws.sent) --> n

        expect(mock_ws, utf8("Hello"))
        expect(mock_ws, utf8("Hello"))
        expect(mock_ws, utf8("Hello"))
        expect(mock_ws, utf8("Hello"))
        expect(mock_ws, utf8("Hello"))
    end

    context("ThrottledWSClient does not throttle close") do
        throttle_time = 0.05
        mock_ws = MockThrottlingWSClient()
        throttled_ws = ThrottledWSClient(mock_ws, throttle_time)

        n = 5
        for i in 1:n
            stop(throttled_ws)
        end

        # Sleep a short amount of time, less than what it would have taken had all requests been
        # throttled.
        sleep(throttle_time)
        # Check that all close requests have reached the mock WSClient.
        @fact mock_ws.closed_called --> n
    end
end


text:
using Test, LinearAlgebra, FFTW
using BSE_k_ISDF

# %% bse problem type
example_path = "diamond_333_20_test/"

N_1d = 20
N_rs = (N_1d, N_1d, N_1d)
N_r = prod(N_rs)
N_core = 0
N_v = 4
N_c = 5
N_ks = (3, 3, 3)
N_k = 27
prob = BSE_k_ISDF.BSEProblemExciting(N_core, N_v, N_c, N_ks, N_rs, example_path)

@test typeof(prob) <: BSE_k_ISDF.AbstractBSEProblem

# test energies and orbitals
E_v, E_c = BSE_k_ISDF.energies(prob)
u_v, u_c = BSE_k_ISDF.orbitals(prob)
Ω0_vol = BSE_k_ISDF.Ω0_volume(prob)

@test BSE_k_ISDF.size(prob) == (N_v, N_c, N_k)
@test BSE_k_ISDF.size_k(prob) == N_ks
@test BSE_k_ISDF.size_r(prob) == N_rs
@test size(E_v) == (N_v, N_k)
@test size(E_c) == (N_c, N_k)
@test size(u_v) == (N_r, N_v, N_k)
@test size(u_c) == (N_r, N_c, N_k)
@test sqrt(Ω0_vol / N_r) * norm(u_v[:, end, 1]) ≈ 1 atol=0.001 # relatively large error, should they be normalized by hand?
@test sqrt(Ω0_vol / N_r) * norm(u_c[:, 1, 1]) ≈ 1 atol=0.001

# test periodicity, test if boundary value is close to mean over the neighbors
for u_i in [u_v, u_c]
    u_i_reshaped = reshape(u_i[:, 1, 1], N_rs)
    N_rs_position = div.(N_rs, 2)

    @test isapprox(u_i_reshaped[1, N_rs_position[2], N_rs_position[3]], 0.5 * (u_i_reshaped[2, N_rs_position[2], N_rs_position[3]] + u_i_reshaped[end, N_rs_position[2], N_rs_position[3]]), rtol = 1e-1)
    @test isapprox(u_i_reshaped[N_rs_position[1], 1, N_rs_position[3]], 0.5 * (u_i_reshaped[N_rs_position[1], 2, N_rs_position[3]] + u_i_reshaped[N_rs_position[1], end, N_rs_position[3]]), rtol = 1e-1)
    @test isapprox(u_i_reshaped[N_rs_position[1], N_rs_position[2], 1], 0.5 * (u_i_reshaped[N_rs_position[1], N_rs_position[2], 2] + u_i_reshaped[N_rs_position[1], N_rs_position[2], end]), rtol = 2e-1)
end

# test lattice
a_mat = BSE_k_ISDF.lattice_matrix(prob)
rl = BSE_k_ISDF.r_lattice(prob)
rc = BSE_k_ISDF.r_cartesian(prob)

@test size(a_mat) == (3, 3)
@test size(rl) == (3, N_r)
@test size(rc) == (3, N_r)
@test rc[:, 1] == a_mat * rl[:, 1]

shift = [0.016666666666666, 0.05, 0.0833333333333]
b_mat = 2 * π * inv(a_mat)
kl = BSE_k_ISDF.k_lattice(prob, shift)
kc = BSE_k_ISDF.k_cartesian(prob, b_mat * shift)

@test size(kl) == (3, N_k)
@test size(kc) == (3, N_k)
@test kc[:, 1] ≈ b_mat * kl[:, 1]

fft_frequencies = BSE_k_ISDF.fftfreq(N_rs...)

@test size(fft_frequencies) == (3, N_r)
@test fft_frequencies[:, 1] == [0, 0, 0]
@test fft_frequencies[:, N_rs[1]] == [-1, 0, 0]
@test fft_frequencies[:, N_rs[1] + 1] == [0, 1, 0]
@test fft_frequencies[:, end] == [-1, -1, -1]
for iG in 1:N_r
    @test BSE_k_ISDF.G_vector_to_index(Int.(fft_frequencies[:, iG]), N_rs) == iG
end

# ISDF
N_μ_vv, N_μ_cc, N_μ_vc = 77, 194, 78

isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)
ζ_vv, ζ_cc, ζ_vc = BSE_k_ISDF.interpolation_vectors(isdf)
u_v_vv, u_c_cc, u_v_vc, u_c_vc = BSE_k_ISDF.interpolation_coefficients(isdf)
N_μ_vv, N_μ_cc, N_μ_vc = BSE_k_ISDF.size(isdf)

@test BSE_k_ISDF.size(isdf) == (N_μ_vv, N_μ_cc, N_μ_vc)
@test size(ζ_vv) == (N_r, N_μ_vv)
@test size(ζ_cc) == (N_r, N_μ_cc)
@test size(ζ_vc) == (N_r, N_μ_vc)
@test size(u_v_vv) == (N_μ_vv, N_v, N_k)
@test size(u_c_cc) == (N_μ_cc, N_c, N_k)
@test size(u_v_vc) == (N_μ_vc, N_v, N_k)
@test size(u_c_vc) == (N_μ_vc, N_c, N_k)

# test the convolution with padding
function w_conv_reference!(b, w, a)
    dim = ndims(a)
    sa = size(a)
    sw = size(w)

    for i in CartesianIndices(sa)
        b[i] = zero(eltype(b))
        for j in CartesianIndices(sa)
            b[i] += w[CartesianIndex(mod1.((i - j).I .+ 1, sw))] * a[j]
        end
    end

    return b
end

conv_size = (9, 11, 13)
double_conv_size = (18, 22, 26)
a = rand(Complex{Float64}, conv_size)
w = rand(Complex{Float64}, double_conv_size)
b = rand(Complex{Float64}, conv_size)
ap = zeros(Complex{Float64}, double_conv_size)
bp = zeros(Complex{Float64}, double_conv_size)
cp = zeros(Complex{Float64}, double_conv_size)
w_hat = zeros(Complex{Float64}, double_conv_size)
p = plan_fft(bp)
p_back = plan_bfft(bp)
@test w_conv_reference!(b, w, a) ≈ BSE_k_ISDF.w_conv!(b, w, a, ap, bp, cp, w_hat, p, p_back)

conv_size = (9, 11, 1) # convolution for 2d k point set
double_conv_size = (18, 22, 1)
a = rand(Complex{Float64}, conv_size)
w = rand(Complex{Float64}, double_conv_size)
b = rand(Complex{Float64}, conv_size)
ap = zeros(Complex{Float64}, double_conv_size)
bp = zeros(Complex{Float64}, double_conv_size)
cp = zeros(Complex{Float64}, double_conv_size)
w_hat = zeros(Complex{Float64}, double_conv_size)
p = plan_fft(bp)
p_back = plan_bfft(bp)
@test w_conv_reference!(b, w, a) ≈ BSE_k_ISDF.w_conv!(b, w, a, ap, bp, cp, w_hat, p, p_back)

# %% Hamiltonian
D = BSE_k_ISDF.setup_D(prob)
@test size(D) == (N_v * N_c * N_k, N_v * N_c * N_k)
@test isdiag(D)
@test D[1, 1] == E_c[1, 1] - E_v[1, 1]
@test D[2, 2] == E_c[1, 1] - E_v[2, 1]

V = BSE_k_ISDF.setup_V(prob, isdf)
V_dense = Matrix(V)
V_exact = BSE_k_ISDF.assemble_exact_V(prob)

# the tolerances are chosen to catch regressions in accuracy
@test size(V) == (N_v * N_c * N_k, N_v * N_c * N_k)
@test ishermitian(V)
@test V_dense ≈ V_dense'
@test isapprox(V_dense, V_exact, atol = 2e-2)
@test isapprox(V_dense, V_exact, rtol = 6e-2)

W = BSE_k_ISDF.setup_W(prob, isdf)
W_dense = Matrix(W)
W_exact = BSE_k_ISDF.assemble_exact_W(prob)

@test size(W) == (N_v * N_c * N_k, N_v * N_c * N_k)
@test ishermitian(W)
@test isapprox(W_dense, W_dense', atol = 9e-3)
@test isapprox(W_dense, W_dense', rtol = 2e-2)
@test isapprox(W_dense, W_exact, atol=5e-2)
@test isapprox(W_dense, W_exact, rtol=6e-2)

H = BSE_k_ISDF.setup_H(prob, isdf)
H_dense = Matrix(H)
H_exact = D + 2 * V_exact - W_exact

@test size(H) == (N_v * N_c * N_k, N_v * N_c * N_k)
@test ishermitian(H)
@test isapprox(H_dense, H_dense', atol = 9e-3)
@test isapprox(H_dense, H_dense', rtol = 5e-4)
@test H_dense ≈ D + 2 * V_dense - W_dense
@test isapprox(H_dense, H_exact, atol=5e-2)
@test isapprox(H_dense, H_exact, rtol=3e-3)

# absorption spectrum

direction = 1
N_iter = 50
σ = 0.0055
g = ω -> 1 / π * σ / (ω^2 + σ^2)
E_min, E_max = 0.0, 1.0
Erange = E_min:0.001:E_max

d = BSE_k_ISDF.optical_absorption_vector(prob, direction)
scaling = 8 * π^2 / (Ω0_vol * N_k)

optical_absorption_lanc = BSE_k_ISDF.lanczos_optical_absorption(prob, isdf, direction, N_iter, g, Erange)
optical_absorption_dense_lanc = BSE_k_ISDF.lanczos_optical_absorption(H_dense, d, N_iter, g, Erange, norm(d)^2 * scaling)
optical_absorption_exact_lanc = BSE_k_ISDF.lanczos_optical_absorption(H_exact, d, N_iter, g, Erange, norm(d)^2 * scaling)

ev_dense, ef_dense = eigen(H_dense)
ev_dense = real.(ev_dense)
weights_dense = [abs2(ef_dense[:, i]' * d) for i in 1:size(ef_dense, 2)]
optical_absorption_dense_eig = BSE_k_ISDF.lanczos_optical_absorption(ev_dense, weights_dense, g, Erange, scaling)

ev_exact, ef_exact = eigen(H_exact)
ev_exact = real.(ev_exact)
weights_exact = [abs2(ef_exact[:, i]' * d) for i in 1:size(ef_exact, 2)]
optical_absorption_exact_eig = BSE_k_ISDF.lanczos_optical_absorption(ev_exact, weights_exact, g, Erange, scaling)

# the tolerances are chosen to catch regressions in accuracy
@test optical_absorption_lanc ≈ optical_absorption_dense_lanc
@test isapprox(optical_absorption_lanc, optical_absorption_exact_lanc, rtol = 1e-1)
@test isapprox(optical_absorption_lanc, optical_absorption_dense_eig, rtol = 7e-2)
@test isapprox(optical_absorption_dense_eig, optical_absorption_exact_eig, rtol = 1e-1)

# resonant and antiresonant peaks, cancellation close to zero
@test BSE_k_ISDF.lanczos_optical_absorption([0.01], [1.0], g, Erange, 1.0) ≈ g.(Erange .- 0.01) .- g.(Erange .+ 0.01)
@test BSE_k_ISDF.lanczos_optical_absorption([0.01], [1.0], g, Erange, 1.0)[1] == 0.0


text:
using VortexDistributions, Distributions, Statistics, BenchmarkTools

include("../util.jl")

function euclid(x, y)
    return sqrt((x[1]-y[1])^2 + (x[2]-y[2])^2)
end

##
Nx = 4000; Ny = 4000
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi = Torus(psi0,x,y)

##

rand_vorts = rand_pointvortex(2, psi)
# p1 = PointVortex(0., 0., 1)
# p2 = PointVortex(0.001, 0.002, 1)

vortex!(psi, rand_vorts)
# vortex!(psi, [p1, p2])

@time vfound1 = findvortices(psi)
@time vfound2 = naivePlaquette(psi.ψ, [psi.x, psi.y], 0)
@time vfound3 = remove_vortices_edge(findvortices_jumps(psi), psi)
@time vfound4 = remove_vortices_edge(findvortices_grid(psi), psi)

# sort!(rand_vorts, by = x -> x.xv)
sort!(vfound1, by = x -> x.xv)
sort!(vfound2, by = x -> x.xv)
sort!(vfound3, by = x -> x.xv)
sort!(vfound4, by = x -> x.xv)

vf1diff = [euclid([rand_vorts[i].xv, rand_vorts[i].yv], [vfound1[i].xv, vfound1[i].yv]) for i in 1:length(rand_vorts)]
vf2diff = [euclid([rand_vorts[i].xv, rand_vorts[i].yv], [vfound2[i].xv, vfound2[i].yv]) for i in 1:length(rand_vorts)]
vf2_3diff = [euclid([vfound2[i].xv, vfound2[i].yv], [vfound3[i].xv, vfound3[i].yv]) for i in 1:length(rand_vorts)]

vf1avg = mean(vf1diff)
vf2avg = mean(vf2diff)
vf2_3avg = mean(vf2_3diff)

naive = benchmark2Dnaive(2)

naive = @benchmark naivePlaquette(psi.ψ, [psi.x, psi.y], 0)
optimised = @benchmark remove_vortices_edge(findvortices_grid(psi), psi)

benchmarks = benchmark2Dnaive2(2)

using JLD2

@save "data/naiveBenchmarks.jld2" benchmarks;


text:
module Dec1Part1
function main(file_path::String)
    lines = readlines(file_path)
    lines = parse.(Int32, lines)
    n_lines = length(lines)
    for i in 1:n_lines
        for j in i + 1:n_lines
            val1 = lines[i]
            val2 = lines[j]
            if val1 + val2 == 2020
                return val1 * val2
            end
        end
    end
end
end


text:
function ThreadedScans.simple!(
    op,
    xs;
    ntasks::Integer = Threads.nthreads(),
)
    if length(xs) <= 2
        if length(xs) == 2
            xs[begin+1] = op(xs[begin], xs[begin+1])
        end
        return xs
    end
    if ntasks < 2
        scan!(op, xs)
        return xs
    end
    nchunks = min(length(xs), ntasks + 1)
    ntasks = nchunks - 1
    @assert ntasks > 1

    chunks = partitionto(xs, ntasks + 1)
    accs = similar(xs, ntasks)
    spawn_foreach(1:ntasks) do i
        if i == 1
            accs[i] = scan!(op, chunks[i])
        else
            accs[i] = foldl(op, chunks[i])
        end
    end
    scan!(op, accs)
    spawn_foreach(1:ntasks) do i
        scan!(op, accs[i], chunks[i+1])
    end

    return xs
end


text:
# this file is to run JWAS
using DataFrames,CSV,Random,Distributions,JWAS,DelimitedFiles, LinearAlgebra
Random.seed!(123)

method            = ARGS[1]
rep               = parse(Int, ARGS[2])  #rep1,...,rep20
pct               = parse(Float64, ARGS[3])

@show method,rep,pct

mainpath="/group/qtlchenggrp/tianjing/res_andrew_1000snp/missing_pattern_legarra/"
cd(mainpath*"$pct/$method/rep$rep/")

############ READ DATA ##########
#read missing index
missingIndex=Int.(vec(readdlm(mainpath*"missing_index/missingIndex.pct$pct.rep$rep.txt")))  #same as before
@show missingIndex[1:3]
@show length(missingIndex)

#read phenotype
phenofile="/group/qtlchenggrp/tianjing/omics_data/andrew_1000snp/data1/y.csv"
phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=["NA"])

#read omics
omicsfile  = "/group/qtlchenggrp/tianjing/omics_data/andrew_1000snp/data1/omics.csv"
omics  = CSV.read(omicsfile,DataFrame,delim = ',',header=true,missingstrings=["NA"]);
insertcols!(omics,2,:y => phenotypes[:,:y]);
for i in 2:size(omics,2) #skip ID
    omics[!,i]= convert(Vector{Union{Missing,Float64}}, omics[!,i]);
end
omics[956:1055,:y] .= missing  #testing dataset
omics[missingIndex,3:end].=missing; #set missing patterns
@show [count(ismissing,col) for col in eachcol(omics)]

omics_name=["gene$i" for i in 1:20]

@show size(omics)
@show omics[1:5,1:5];
@show omics[end-4:end,1:5];

#read genotype    # methodAll=("omics_bayescpi" "omics_gblup")
if occursin("gblup",method)  #GRM for gblup_fixedvar
	genofile   = "/group/qtlchenggrp/tianjing/omics_data/andrew_1000snp/data1/GRM.csv"
	genotypes  = get_genotypes(genofile;separator=',',method="GBLUP",header=false);
else #genotype for rrblup, bayescpi
	genofile   = "/group/qtlchenggrp/tianjing/omics_data/andrew_1000snp/data1/geno1k.csv"
	if occursin("rrblup",method)
		genotypes  = get_genotypes(genofile,separator=',',method="RR-BLUP");
	elseif occursin("bayescpi",method)
		genotypes  = get_genotypes(genofile,separator=',',method="BayesC");
	elseif occursin("bayesa",method)
		genotypes  = get_genotypes(genofile,separator=',',method="BayesA");
	elseif occursin("bayesb",method)
		genotypes  = get_genotypes(genofile,separator=',',method="BayesB");
	elseif occursin("bayesl",method)
		genotypes  = get_genotypes(genofile,separator=',',method="BayesL");
	end
end

model_equation = "y = intercept + genotypes";

model = build_model(model_equation;
                    num_hidden_nodes=20,
                    latent_traits=omics_name,
                    nonlinear_function="linear");
out   = runMCMC(model,omics,chain_length=5000,printout_model_info=false);

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
function elapsed_time_reporting_criterion(report_after::Number, channel::RemoteChannel)
    previous_time = time()
    while true
        current_time = time()
        try
            if current_time - previous_time >= report_after
                previous_time = time()
                put!(channel, true)
            else
                put!(channel, false)
            end
        catch
            break
        end
    end
end

function iterations_reporting_criterion(report_after::Number, channel::RemoteChannel)
    previous_iteration = 0
    current_iteration  = previous_iteration
    while true
        current_iteration += 1
        try
            if current_iteration - previous_iteration >= report_after
                previous_iteration = current_iteration
                put!(channel, true)
            else
                put!(channel, false)
            end
        catch
            break
        end
    end
end


text:
function run_function(f, timeline, modules, name::Union{Symbol,AbstractString}, args::Vector)
    # Find the function ip
    (h, base, sym)  = Gallium.lookup_sym(timeline, modules, name)
    addr = Gallium.compute_symbol_value(h, base, sym)
    run_function(f, timeline, addr, args)
end

run_function(f, timeline, modules, name, args::Union{Integer,Ptr}) = run_function(f, timeline, modules, name, [args])

const args_regs = [:rdi, :rsi, :rdx, :rcx, :r8, :r9]
function run_function(f, timeline, addr::Union{Integer,RemotePtr}, args::Vector)
    diversion = prepare_remote_execution(timeline)
    icxx"$diversion->set_visible_execution(true);"

    addr = UInt64(addr)

    # Pick an arbitrary task with the same tgid as our current task to run our
    # expression
    task = icxx"""
        for (auto &t : $diversion->tasks())
            if (t.second->tgid() == $(current_task(timeline))->tgid())
                return t.second;
        return (rr::Task *)nullptr;
    """
    regs = icxx"$task->regs();"

    # Set up the call frame
    Gallium.set_ip!(regs, addr)
    new_rsp = Gallium.get_dwarf(regs, :rsp)-250
    new_rsp -= (new_rsp % 16)
    new_rsp += 8
    Gallium.set_dwarf!(regs, :rsp, new_rsp)

    # Set return address to 0
    Gallium.store!(task, Gallium.RemotePtr{UInt64}(new_rsp), UInt64(0))

    # Set up arguments
    for (i,val) in enumerate(args)
        i > length(args_regs) && error("Too many arguments")
        Gallium.set_dwarf!(regs, args_regs[i], UInt64(val))
    end

    # Writes registers back to task
    icxx"$task->set_regs($regs);"

    # Alright, let's go
    while true
        res = icxx"$diversion->diversion_step($task);"
        if icxx"$res.status != rr::DiversionSession::DIVERSION_CONTINUE;" ||
            icxx"$res.break_status.signal != 0;"
            break
        end
    end

    f(task)
end

# Remote C execution
include("remoteir.jl")

@cxxm "uint64_t GalliumCallbacks::allocateMem(uint32_t kind, uint64_t Size, uint64_t Align)" begin
    global code_mem, ro_mem, rw_mem
    if kind == icxx"llvm::sys::Memory::MF_EXEC | llvm::sys::Memory::MF_READ;"
        code_mem += code_mem % Align
        ret = code_mem
        code_mem += Size
        return UInt64(ret)
    elseif kind == icxx"llvm::sys::Memory::MF_READ;"
        ro_mem += ro_mem % Align
        ret = ro_mem
        ro_mem += Size
        return UInt64(ret)
    elseif kind == icxx"llvm::sys::Memory::MF_READ | llvm::sys::Memory::MF_WRITE;"
        rw_mem += rw_mem % Align
        ret = rw_mem
        rw_mem += Size
        return UInt64(ret)
    else
        error("Unknown kind")
    end
end

@cxxm "void GalliumCallbacks::writeMem(uint64_t remote, uint8_t *localaddr, size_t size)" begin
    session = unsafe_pointer_to_objref(icxx"$this->session;")
    Gallium.store!(session,RemotePtr{UInt8}(remote),
        unsafe_wrap(Array, localaddr, size, false))
end

function lookup_external_symbol(modules, name)::UInt64
    global data_buffer_start
    name == "data_buffer_start" && return data_buffer_start
    h,base,sym = Gallium.lookup_sym(timeline, modules, name)
    ret = ObjFileBase.compute_symbol_value(h, base, sym)
    return ret
end

# Now allocate some memory for the JIT
function create_remote_jit(timeline, near_addr)
    always_free_addresses = icxx"""
        rr::TraceReader reader{$(current_session(timeline))->trace_reader()};
        reader.rewind();
        rr::ReplaySession::always_free_address_space(reader);
    """

    start_addr = icxx"""
        for (auto range : $always_free_addresses) {
            if (std::abs((intptr_t)(range.start().as_int() - $(UInt64(near_addr)))) < INT32_MAX/2) {
                return range.start().as_int();
            } else if (std::abs((intptr_t)(range.end().as_int() - 0x1000 - $(UInt64(near_addr)))) < INT32_MAX/2) {
                return range.end().as_int() -  0x40000;
            }
        }
        return (uint64_t)0;
    """
    global code_mem, ro_mem, rw_mem
    region_size = 0x10000 # 16 pages
    code_mem = start_addr
    icxx"""
    rr::AutoRemoteSyscalls remote($(current_task(current_session(timeline))));
    remote.infallible_mmap_syscall($code_mem, $region_size,
        PROT_EXEC | PROT_READ,
        MAP_ANONYMOUS | MAP_PRIVATE | MAP_FIXED, -1, 0);
    """
    ro_mem = start_addr + region_size
    icxx"""
    rr::AutoRemoteSyscalls remote($(current_task(current_session(timeline))));
    remote.infallible_mmap_syscall($ro_mem, $region_size,
        PROT_READ,
        MAP_ANONYMOUS | MAP_PRIVATE | MAP_FIXED, -1, 0);
    """
    rw_mem = start_addr + 2region_size
    icxx"""
    rr::AutoRemoteSyscalls remote($(current_task(current_session(timeline))));
    remote.infallible_mmap_syscall($rw_mem, $region_size,
        PROT_READ | PROT_WRITE,
        MAP_ANONYMOUS | MAP_PRIVATE | MAP_FIXED, -1, 0);
    """
    stack_mem = 0x700001000
    icxx"""
    rr::AutoRemoteSyscalls remote($(current_task(current_session(timeline))));
    remote.infallible_mmap_syscall($stack_mem, $region_size,
        PROT_READ | PROT_WRITE,
        MAP_ANONYMOUS | MAP_PRIVATE | MAP_FIXED, -1, 0);
    """

    memm = icxx"""
    GalliumCallbacks callbacks;
    callbacks.session = $(pointer_from_objref(timeline));
    new RCMemoryManager(std::move(callbacks));
    """
    callbacks = icxx"&$memm->Client;"
    jit = icxx"new RemoteJIT(*llvm::EngineBuilder().selectTarget(), $memm);"
    jit, callbacks
end

cxxinclude(joinpath(dirname(@__FILE__),"FunctionMover.cpp"))
TargetClang = Cxx.new_clang_instance(false)

function allocate_code(callbacks, code)
    addr = icxx"$callbacks->allocateMem(
        llvm::sys::Memory::MF_EXEC | llvm::sys::Memory::MF_READ,
        $(sizeof(code)), 0x16);"
    icxx"$callbacks->writeMem($addr,(uint8_t*)$(pointer(code)),$(sizeof(code)));"
    addr
end

function rewrite_instruction(inst, end_ip)
    # Adjust RIP-relative mov
    length(inst) == 7 || return inst
    # REX.W(R)
    ((inst[1] & 0b11111011) == 0b01001000) || return inst
    # opcode
    (inst[2] == 0x8b) || return inst
    # MODRM, mod = 0b00, r/m = 0b101
    ((inst[3] & 0b11000111) == 0b00000101) || return inst
    # Ok, we have a RIP-relative MOV
    addr = UInt64(end_ip+reinterpret(Int32, inst[4:7])[])
    # Instead encode an absolute MOV
    reg = (inst[3]&0b111000) >> 3
    [
        # mov $abs, %reg
        inst[1]; 0xb8+reg; reinterpret(UInt8,[addr]);
        # move (%reg), %reg
        inst[1]; 0x8b; (reg<<3)|reg
    ]
end

function rewrite_instructions(insts, start_ip)
    rewritten_insts = UInt8[]
    current_end_ip = start_ip
    while current_end_ip < start_ip+length(insts)
        next_inst = extract_next_inst(insts[(1+current_end_ip-start_ip):end])
        current_end_ip += length(next_inst)
        next_inst = rewrite_instruction(next_inst, current_end_ip)
        append!(rewritten_insts, next_inst)
    end
    rewritten_insts
end

Base.unsafe_string(ref::vcpp"llvm::StringRef") =
    unsafe_string(icxx"$ref.data();", icxx"$ref.size();")
function run_func(timeline, jit, callbacks,
        fname::Union{pcpp"clang::Decl",pcpp"clang::FunctionDecl"}, retT=Void, TargetClang = TargetClang, args = Any[])
    run_func(timeline, jit, callbacks,
        unsafe_string(icxx"cast<clang::NamedDecl>($fname)->getName();"),
        retT, TargetClang, args)
end

function run_func(timeline, jit, callbacks,
        fname::pcpp"llvm::Function", retT=Void, TargetClang = TargetClang, args = Any[])
    run_func(timeline, jit, callbacks,
        unsafe_string(icxx"$fname->getName();"), retT, TargetClang, args)
end

function run_func(timeline, jit, callbacks, fname, retT=Void, TargetClang = TargetClang, args=Any[])
    shadowmod = Cxx.instance(TargetClang).shadow
    targetmod = icxx"""new llvm::Module("Target Module", $shadowmod->getContext());"""
    icxx"""$targetmod->setDataLayout($shadowmod->getDataLayout());"""

    mover = icxx"new FunctionMover2($targetmod);"

    F = icxx"$shadowmod->getFunction($(pointer(fname)));"
    @assert F != C_NULL
    icxx"MapFunction($F, $mover);"

    icxx"""
    $jit->addModule(std::unique_ptr<llvm::Module>($targetmod));
    """

    addr = icxx"""$jit->findSymbol($(pointer(fname)), false).getAddress();"""
    @assert UInt64(addr) != 0

    run_function(timeline, UInt64(addr), args) do task
        regs = icxx"$task->regs();"
        x = Gallium.ip(regs)
        @assert UInt64(x) == 0
        sizeof(retT) == 0 && return retT.instance
        reinterpret(retT, [Gallium.get_dwarf(regs, :rax)])[]
    end
end

function trace_func(jit, callbacks, fname, entry_func, exit_func = "")
    h,base,sym = Gallium.lookup_sym(timeline, modules, fname)
    hook_addr = ObjFileBase.compute_symbol_value(h, base, sym)

    shadowmod = Cxx.instance(TargetClang).shadow
    targetmod = icxx"""new llvm::Module("Target Module", $shadowmod->getContext());"""
    icxx"""$targetmod->setDataLayout($shadowmod->getDataLayout());"""

    mover = icxx"new FunctionMover2($targetmod);"

    F = icxx"$shadowmod->getFunction($(pointer(entry_func)));"
    @assert F != 0
    icxx"MapFunction($F, $mover);"

    if !isempty(exit_func)
        F = icxx"$shadowmod->getFunction($(pointer(exit_func)));"
        @assert F != 0
        icxx"MapFunction($F, $mover);"
    end

    icxx"""
    $jit->addModule(std::unique_ptr<llvm::Module>($targetmod));
    """

    entry_hook = icxx"""$jit->findSymbol($(pointer(entry_func))).getAddress();"""
    exit_hook = isempty(exit_func) ? UInt64(0) :
        icxx"""$jit->findSymbol($(pointer(exit_func))).getAddress();"""

    hook_template = Gallium.Hooking.hook_asm_template(UInt64(0),
        UInt64(0); call = false)

    orig_bytes = Gallium.load(task, RemotePtr{UInt8}(hook_addr), length(hook_template)+15)
    nbytes = Gallium.Hooking.determine_nbytes_to_replace(length(hook_template), orig_bytes)

    ret_jmp_addr = isempty(exit_func) ? UInt64(0) : allocate_code(callbacks, [
        Gallium.Hooking.return_hook_template(0x700011000, exit_hook);
    ])

    jmp_addr = allocate_code(callbacks, [
        Gallium.Hooking.instrument_jmp_template(0x700011000,entry_hook,ret_jmp_addr);
        Gallium.Hooking.hook_tail_template(
            rewrite_instructions(orig_bytes[1:nbytes],UInt(hook_addr)),
            UInt(hook_addr)+nbytes)
    ])

    hook_template = Gallium.Hooking.hook_asm_template(UInt64(hook_addr),
        UInt64(jmp_addr); call = false)

    replacement = [hook_template; zeros(UInt8,nbytes-length(hook_template))]
    Gallium.store!(task, RemotePtr{UInt8}(hook_addr), replacement)

end


text:
# load requried packages
using IncrementalInference
using Test

## parameters

lm_prior_noise = 0.01
meas_noise = 0.25
odom_noise = 0.1
n_samples = 200

# initialize mean landmark locations
l0 = 0.0
l1 = 10.0
l2 = 20.0
l3 = 40.0

# "Ground-truth" robot poses
x0 = 0.0
x1 = 10.0
x2 = 20.0
x3 = 40.0

##

@testset "Basic 3 door, trinary multihypothesis tests..." begin

## Initialize empty factor graph
fg = initfg()
getSolverParams(fg).N = n_samples
getSolverParams(fg).gibbsIters = 5

# Place strong prior on locations of three "doors"
addVariable!(fg, :l0, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l0], Prior(Normal(l0, lm_prior_noise)))

addVariable!(fg, :l1, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l1], Prior(Normal(l1, lm_prior_noise)))

addVariable!(fg, :l2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l2], Prior(Normal(l2, lm_prior_noise)))

addVariable!(fg, :l3, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l3], Prior(Normal(l3, lm_prior_noise)))

# Add first pose
addVariable!(fg, :x0, ContinuousScalar, N=n_samples)

# Make first "door" measurement
addFactor!(fg, [:x0; :l0; :l1; :l2; :l3], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; (1/4 for _=1:4)...])

# Add second pose
addVariable!(fg, :x1, ContinuousScalar, N=n_samples)

# Gaussian transition model
addFactor!(fg, [:x0; :x1], LinearRelative(Normal(x1-x0, odom_noise)))

# Make second "door" measurement
addFactor!(fg, [:x1; :l0; :l1; :l2; :l3], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; (1/4 for _=1:4)...])

##

solveGraph!(fg)

##

# check there is enough likelihood in the right places
@test 0.1 < getBelief(fg, :x0)([l0])[1]
@test 0.1 < getBelief(fg, :x0)([l1])[1]
@test getBelief(fg, :x0)([l2])[1] < 0.3

@test getBelief(fg, :x1)([l0])[1] < 0.3
@test 0.1 < getBelief(fg, :x1)([l1])[1]
@test 0.1 < getBelief(fg, :x1)([l2])[1]

##

for i in 1:3
  solveGraph!(fg);
end

##

# check there is enough likelihood in the right places
@test 0.1 < getBelief(fg, :x0)([l0])[1]
@test 0.1 < getBelief(fg, :x0)([l1])[1]
@test getBelief(fg, :x0)([l2])[1] < 0.03

@test getBelief(fg, :x1)([l0])[1] < 0.03
@test 0.1 < getBelief(fg, :x1)([l1])[1]
@test 0.1 < getBelief(fg, :x1)([l2])[1]

## Add one more pose/odometry to invoke issue #236

# Add third pose
addVariable!(fg, :x2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(x2-x1, odom_noise)))

## Solve graph

tree = solveTree!(fg)
# drawGraph(fg)
# drawTree(tree, show=true)

##

# check there is enough likelihood in the right places
@test 0.1 < getBelief(fg, :x0)([l0])[1]
@test 0.1 < getBelief(fg, :x0)([l1])[1]
@test getBelief(fg, :x0)([l2])[1] < 0.03

@test getBelief(fg, :x1)([l0])[1] < 0.03
@test 0.1 < getBelief(fg, :x1)([l1])[1]
@test 0.1 < getBelief(fg, :x1)([l2])[1]

dx = x2-x1
@test getBelief(fg, :x2)([l0 + dx])[1] < 0.03
@test 0.1 < getBelief(fg, :x2)([l1 + dx])[1]
@test 0.1 < getBelief(fg, :x2)([l2 + dx])[1]

##

# Add third pose
addVariable!(fg, :x3, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x2; :x3], LinearRelative(Normal(x3-x2, odom_noise)))

# Make third "door" measurement
addFactor!(fg, [:x3; :l0; :l1; :l2; :l3], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; (1/4 for _=1:4)...])

##

solveGraph!(fg)

##

@test isapprox(mean(getBelief(fg, :x0))[1], x0; atol = 2.0)
@test isapprox(mean(getBelief(fg, :x1))[1], x1; atol = 2.0)
@test isapprox(mean(getBelief(fg, :x2))[1], x2; atol = 2.0)
@test isapprox(mean(getBelief(fg, :x3))[1], x3; atol = 2.0)

@test isapprox(mean(getBelief(fg, :l0))[1], l0; atol = 3.0)
@test isapprox(mean(getBelief(fg, :l1))[1], l1; atol = 3.0)
@test isapprox(mean(getBelief(fg, :l2))[1], l2; atol = 3.0)
@test isapprox(mean(getBelief(fg, :l3))[1], l3; atol = 3.0)

##

# check the PPEs are the same
@test isapprox(getPPE(fg, :x0).suggested[1], x0; atol = 2.0)
@test isapprox(getPPE(fg, :x1).suggested[1], x1; atol = 2.0)
@test isapprox(getPPE(fg, :x2).suggested[1], x2; atol = 2.0)
@test isapprox(getPPE(fg, :x3).suggested[1], x3; atol = 2.0)

@test isapprox(getPPE(fg, :l0).suggested[1], l0; atol = 3.0)
@test isapprox(getPPE(fg, :l1).suggested[1], l1; atol = 3.0)
@test isapprox(getPPE(fg, :l2).suggested[1], l2; atol = 3.0)
@test isapprox(getPPE(fg, :l3).suggested[1], l3; atol = 3.0)

##

end

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)
#
# plotBelief(fg, sortDFG(ls(fg, r"x")))

#


text:
# In this example we show how to pass a DMRGObserver to
# the dmrg function which allows tracking energy convergence and
# convergence of local operators.
using ITensors

"""
  Get MPO of transverse field Ising model Hamiltonian with field strength h
"""
function tfimMPO(sites,
                 h::Float64)
  # Input operator terms which define a Hamiltonian
  N = length(sites)
  ampo = AutoMPO()
  for j=1:N-1
    ampo += -1,"Sz",j,"Sz",j+1
  end
  for j=1:N
    ampo += h,"Sx",j
  end
  # Convert these terms to an MPO tensor network
  return MPO(ampo,sites)
end

let
  N = 100
  sites = siteinds("S=1/2",N)
  psi0 = randomMPS(sites,10)

  # define parameters for DMRG sweeps
  sweeps = Sweeps(15)
  maxdim!(sweeps, 10,20,100,100,200)
  cutoff!(sweeps, 1E-10)

  #=
  create observer which will measure Sᶻ at each
  site during the dmrg sweeps and track energies after each sweep.
  in addition it will stop the computation if energy converges within
  1E-7 tolerance
  =#
  let
    Sz_observer = DMRGObserver(["Sz"],sites,energy_tol=1E-7)

    # we will now run DMRG calculation for different values
    # of the transverse field and check how local observables
    # converge to their ground state values

    println("Running DMRG for TFIM with h=0.1")
    println("================================")
    H = tfimMPO(sites,0.1)
    energy, psi = dmrg(H,psi0,sweeps,observer=Sz_observer)

    for (i,Szs) in enumerate(measurements(Sz_observer)["Sz"])
      println("<Σ Sz> after sweep $i = ", sum(Szs)/N)
    end
  end

  let
    println("\nRunning DMRG for TFIM with h=0.5 (critical point)")
    println("================================")
    Sz_observer= DMRGObserver(["Sz"],sites,energy_tol=1E-7)
    H = tfimMPO(sites,0.5)
    energy, psi = dmrg(H,psi0,sweeps,observer=Sz_observer)

    for (i,Szs) in enumerate(measurements(Sz_observer)["Sz"])
      println("<Σ Sz> after sweep $i = ", sum(Szs)/N)
    end
  end

  let
    println("\nRunning DMRG for TFIM with h=5.")
    println("================================")
    Sz_Sx_observer= DMRGObserver(["Sz","Sx"],sites,energy_tol=1E-7)
    H = tfimMPO(sites,5.0)
    energy, psi = dmrg(H,psi0,sweeps,observer=Sz_Sx_observer)

    for (i,Szs) in enumerate(measurements(Sz_Sx_observer)["Sz"])
      println("<Σ Sz> after sweep $i = ", sum(Szs)/N)
    end
    println()
    for (i,Sxs) in enumerate(measurements(Sz_Sx_observer)["Sx"])
      println("<Σ Sx> after sweep $i = ", sum(Sxs)/N)
    end
  end
end


text:
using Test
using ModelingToolkit, OrdinaryDiffEq

@isdefined(t) || @parameters t
@connector function Pin(; name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name = name)
end

function Ground(; name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name = name), g)
end

function OnePort(; name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i]
    compose(ODESystem(eqs, t, sts, []; name = name), p, n)
end

function Resistor(; name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R = R
    eqs = [
        v ~ i * R,
    ]
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

function Capacitor(; name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C = C
    D = Differential(t)
    eqs = [
        D(v) ~ i / C,
    ]
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

function ConstantVoltage(; name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V = V
    eqs = [
        V ~ v,
    ]
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

function Inductor(; name, L = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters L = L
    D = Differential(t)
    eqs = [
        D(i) ~ v / L,
    ]
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

@connector function HeatPort(; name)
    @variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow]
    ODESystem(Equation[], t, [T, Q_flow], [], name = name)
end

function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0)
    @named p = Pin()
    @named n = Pin()
    @named h = HeatPort()
    @variables v(t) RTherm(t)
    @parameters R=R TAmbient=TAmbient alpha=alpha
    eqs = [RTherm ~ R * (1 + alpha * (h.T - TAmbient))
           v ~ p.i * RTherm
           h.Q_flow ~ -v * p.i # -LossPower
           v ~ p.v - n.v
           0 ~ p.i + n.i]
    compose(ODESystem(eqs, t, [v, RTherm], [R, TAmbient, alpha],
                      name = name), p, n, h)
end

function HeatCapacitor(; name, rho = 8050, V = 1, cp = 460, TAmbient = 293.15)
    @parameters rho=rho V=V cp=cp
    C = rho * V * cp
    @named h = HeatPort()
    D = Differential(t)
    eqs = [
        D(h.T) ~ h.Q_flow / C,
    ]
    compose(ODESystem(eqs, t, [], [rho, V, cp],
                      name = name), h)
end


text:
# Several functions in this file are adapted from Frank T. Bergman
# Date: 2019
# Availability: https://groups.google.com/forum/#!topic/sbml-discuss/inS4Lzp3Ri8 or
# https://www.dropbox.com/s/2bfpiausejp0gd0/convert_reactions.py?dl=0
###################################################################################

@ModelingToolkit.parameters t
_Differential = ModelingToolkit.Differential(t)
export _Differential, t

"""
    simulatesbml(sbmlfile::String,tspan;saveat::Real=1.,jac::Bool=true,
                 solver=OrdinaryDiffEq.Tsit5())

Simulate model from SBML file and return an ODESolution.

# Arguments
- `sbmlfile::String`: path to SBML file.
- `tspan`: The solution u(t) will be computed for tspan[1] ≤ t ≤ tspan[2].
- `saveat::Real=1.`: time intervals of u(t) output.
- `jac::Bool=true`: Jacobian matrix.
- `solver=OrdinaryDiffEq.Tsit5()`: Integration algorithm.

# Example
```julia-repl
julia> sol = simulatesbml("mymodel.sbml",(0.0,1.0))
retcode: Success
Interpolation: 1st order linear
t: 2-element Array{Float64,1}:
 0.0
 1.0
u: 2-element Array{Array{Float64,1},1}:
 [0.0, 1.0]
 [0.4, 0.6]
```
"""
function simulatesbml(sbmlfile::String,tspan;saveat::Real=1.,jac::Bool=true,
                      solver=OrdinaryDiffEq.Tsit5())
    prob = sbml2odeproblem(sbmlfile,tspan,jac=jac)
    OrdinaryDiffEq.solve(prob,solver,saveat=saveat)
end

"""
    sbml2odeproblem(sbmlfile::String,tspan;jac::Bool=true)

Convert model from SBML file to an `ODEProblem` and return it.

# Arguments
- `sbmlfile::String`: path to SBML file.
- `tspan`: The solution u(t) will be computed for tspan[1] ≤ t ≤ tspan[2].
- `jac::Bool=true`: Jacobian matrix.

# Example
```julia-repl
julia> prob = sbml2odeproblem(SBML_FILE,(0.0,10.0))
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: [0.0, 1.0]
```
"""
function sbml2odeproblem(sbmlfile::String,tspan;jac::Bool=true)
    model = SbmlInterface.getmodel(sbmlfile)
    p = getparameters(model)
    u0 = getinitialconditions(model)
    sys= sbml2odesystem(sbmlfile)
    ModelingToolkit.ODEProblem(sys, Pair[], tspan)
end

"""
    sbml2odesystem(sbmlfile::String)

Given an `sbmlfile` return an `ODESystem` with `default_u0` and `default_p` set.

# Example
```julia-repl
julia> sys = sbml2odesystem(SBML_FILE)
Model ##ODESystem#270 with 2 equations
States (2):
  B(t) [defaults to 0.0]
  A(t) [defaults to 1.0]
Parameters (3):
  k2 [defaults to 0.6]
  compartment [defaults to 1.0]
  k1 [defaults to 0.8]
```
"""
function sbml2odesystem(sbmlfile::String)
    if !isfile(sbmlfile)
        throw(DomainError("`sbmlfile` is not a file"))
    elseif !(splitext(sbmlfile)[2] in (".xml", ".sbml"))
        throw(DomainError("`sbmlfile` must have `.xml` or `.sbml` ending"))
    end
    model = getmodel(sbmlfile)
    p = getparameters(model)
    u0 = getinitialconditions(model)
    rxs = getreactions(model)
    rs  = ModelingToolkit.ReactionSystem(rxs, t, [item.first for item in u0], [item.first for item in p])
    odesys = convert(ModelingToolkit.ODESystem, rs)
    eqs = odesys.eqs
    sys = ModelingToolkit.ODESystem(eqs; default_p=p, default_u0=u0)
    sys
end

"""
    getmodel(sbmlfile::String)

Given an `sbmlfile` return a PyObject (libsbml.Model).

# Example
```julia-repl
julia> model = getmodel(SBML_FILE)
PyObject <Model conversion_reaction_0 "Conversion Reaction 0">
```
"""
function getmodel(sbmlfile::String)
    doc = libsbml.readSBMLFromFile(sbmlfile)
    if doc.getNumErrors(libsbml.LIBSBML_SEV_FATAL) == 1
            @warn("Encountered serious errors while reading file")
            @error(doc.getErrorLog().toString())
    end

    # clear errors
    doc.getErrorLog().clearLog()

    # perform conversions
    props = libsbml.ConversionProperties()
    props.addOption("promoteLocalParameters", true)

    if doc.convert(props) != libsbml.LIBSBML_OPERATION_SUCCESS
        @warn("The document could not be converted")
        @warn(doc.getErrorLog().toString())
    end

    props = libsbml.ConversionProperties()
    props.addOption("expandFunctionDefinitions", true)  # Todo: ask PEtab developers set this to `True` when creating `petab.problem.Problem()`

    if doc.convert(props) != libsbml.LIBSBML_OPERATION_SUCCESS
        @warn("The document could not be converted")
        @warn(doc.getErrorLog().toString())
    end

    mod = doc.getModel()
end

"""
    getreactions(model)::Array

Create array of Reactions from a libsbml model.

# Example
```julia-repl
julia> rxs = getreactions(model)
2-element Array{ModelingToolkit.Reaction,1}:
 ModelingToolkit.Reaction{Any,Float64}(compartment*k1*A(t), SymbolicUtils.Term{Real}[A(t)],
                                       SymbolicUtils.Term{Real}[B(t)], [1.0], [1.0],
                                       Pair{Any,Float64}[B(t) => 1.0, A(t) => -1.0], true)
 ModelingToolkit.Reaction{Any,Float64}(compartment*k2*B(t), SymbolicUtils.Term{Real}[B(t)],
                                       SymbolicUtils.Term{Real}[A(t)], [1.0], [1.0],
                                       Pair{Any,Float64}[B(t) => -1.0, A(t) => 1.0], true)
```
"""
function getreactions(model)::Array
    assignments = []
    for a in model.getListOfRules()
        if (a.getMath().getName() == "piecewise") # || (a.getMath().getName() == nothing)
            @error("Assignment rules are not implemented.")
        end
        push!(assignments, a.getId() => a.getMath().getName())
    end

    getparameters(model)
    getinitialconditions(model)

    rxs = ModelingToolkit.Reaction[]
    # reactions = Dict()  # Dict of reaction => kinetic formula
    for i in 0:model.getNumReactions()-1
        reaction = model.getReaction(i)
        kineticlaw = reaction.getKineticLaw()
        kineticformula = kineticlaw.getFormula()
        kineticformula = eval(Meta.parse(replace(kineticformula, "pow"=>"^")))
        substrates = [eval(Meta.parse(substrate.getSpecies())) for substrate in reaction.getListOfReactants()]
        substoich = [substrate.getStoichiometry() for substrate in reaction.getListOfReactants()]
        products = [eval(Meta.parse(product.getSpecies())) for product in reaction.getListOfProducts()]
        prodstoich = [product.getStoichiometry() for product in reaction.getListOfProducts()]
        mtkreaction = ModelingToolkit.Reaction(kineticformula, substrates, products, substoich, prodstoich, only_use_rate=true)
        push!(rxs, mtkreaction)
    end
    rxs
end

"""
    getparameters(model)

Create array of parameters from a libsbml Model.

# Example
```julia-repl
julia> pars = getparameters(model)
3-element Array{Pair{Num,Float64},1}:
          k1 => 0.8
          k2 => 0.6
 compartment => 1.0
```
"""
function getparameters(model)
    parameters = Pair{ModelingToolkit.Num,Float64}[]
    for i in 0:model.getNumParameters()-1
        par = model.getParameter(i)
        parname = par.getId()
        parval = par.getValue()
        eval(Meta.parse("@ModelingToolkit.parameters $parname"))
        eval(Meta.parse("export $parname"))
        push!(parameters, eval(Meta.parse("$parname => $parval")))
    end
    for i in 0:model.getNumCompartments()-1
        comp = model.getCompartment(i)
        compname = comp.getId()
        compval = comp.getSize()
        eval(Meta.parse("@ModelingToolkit.parameters $compname"))
        eval(Meta.parse("export $compname"))
        push!(parameters, eval(Meta.parse("$compname => $compval")))
    end
    parameters
end

"""
    getinitialconditions(model)

Create array of initial conditions from a libsbml Model.

# Example
```julia-repl
julia> u0 = getinitialconditions(model)
2-element Array{Pair{Num,Float64},1}:
 A(t) => 1.0
 B(t) => 0.0
```
"""
function getinitialconditions(model)
    eval(Meta.parse("@ModelingToolkit.parameters t"))

    initialassignments = Dict()
    for var in model.getListOfInitialAssignments()
        varname = var.getId()
        varval = var.getMath().getName()
        if varval isa String
            varval = model.getParameter(varval).getValue()
        end
        initialassignments[varname] = varval
    end

    initialconditions = Pair{ModelingToolkit.Num,Float64}[]
    for i in 0:model.getNumSpecies()-1
        var = model.getSpecies(i)
        varname = var.getId()
        if var.isSetInitialConcentration()
            varval = var.getInitialConcentration()
        elseif var.isSetInitialAmount()
            varval = var.getInitialAmount()
        else
            varval = initialassignments[varname]
        end
        eval(Meta.parse("@ModelingToolkit.variables $varname(t)"))
        eval(Meta.parse("export $varname"))
        push!(initialconditions, eval(Meta.parse("$varname => $varval")))
    end
    initialconditions
end


text:
u01  = GpuArray(zeros(M,N))
u02  = GPUArray(zeros(M,N))
u03  = GPUArray(0.2ones(M,N))
u04  = GPUArray(zeros(M,N))

u0 = ArrayPartition(u01,u02,u03,u04)
function hox_krox(du,u,p,t)
  RAout,RAin,gₕ,gₖ = u.x
  dRAout,dRAin,dgₕ,dgₖ = du.x
  Cyp .= !NCYP_kmax.*(kdeg.*γ.*RAin.^n)./(1.+γ.*RAin.^n.+f₀.*exp.(-λ.*(xf.-X))) .+
         NCYP_kmax.*kmax
  dRAout .= DRA.*(RAout*Ax + Ay*RAout) .+ VRA .- (1.+β).*kA.*RAout .+ kA.*RAin
  dRAin  .= kA.*RAout .- (kA .+ Cyp).*RAin
  if t > 500
    dgₕ .= (cₕ.*gₕ.^nₕ .+ (κₕ*RAin).^m) ./ (1 .+ cₕ.*gₕ.^nₕ .+ cₖ.*gₖ.^nₖ .+ (κₕ*RAin).^m) - dₕ.*gₕ
    dgₖ .= (cₖ.*gₖ.^nₖ .+ (κₖ*RAin).^m) ./ (1 .+ cₕ.*gₕ.^nₕ .+ cₖ.*gₖ.^nₖ .+ (κₖ*RAin).^m) - dₖ.*gₖ
  end
  nothing
end
function hox_krox_noise(du,u,p,t)
  RAout,RAin,gₕ,gₖ = u.x
  dRAout,dRAin,dgₕ,dgₖ = du.x
  dRAout .= ϵout .* RAout
  dRAin  .= ϵin  .* RAin
  if t > 500
    dgₕ .= aₕ .* gₕ
    dgₖ .= aₖ .* gₖ
  end
  nothing
end

prob = SDEProblem(hox_krox,hox_krox_noise,u0,(0.0,1000.0))
sol = solve(prob,SRIW1())

max_stars_repo_name: ChrisRackauckas/RAMeanVar

max_issues_repo_name: ChrisRackauckas/RAMeanVar

max_forks_repo_name: ChrisRackauckas/RAMeanVar


text:
module chapters_callbacks

using Dash, DashHtmlComponents, DashCoreComponents

include("../../utils.jl")

include("../../reusable_components/Header.jl")
include("../../reusable_components/Example.jl")

export examples

examples_path = joinpath(@__DIR__, "examples")

hello_slider = LoadExampleCode(string(examples_path, "/hello-slider.jl"))

multi_inputs = LoadExampleCode(string(examples_path, "/multi-inputs.jl"))

multi_outputs = LoadExampleCode(string(examples_path, "/multi-outputs.jl"))

simple_callback = LoadExampleCode(string(examples_path, "/simple-callback.jl"))

basic_input = LoadExampleCode(string(examples_path, "/basic-input.jl"))

basic_state = LoadExampleCode(string(examples_path, "/basic-state.jl"))

callback_chain = LoadExampleCode(string(examples_path, "/getting-started-callback-chain.jl"))

examples = [hello_slider, multi_inputs, multi_outputs, simple_callback, basic_input, basic_state, callback_chain]

app =  dash()
basic_input.callback!(app)
basic_state.callback!(app)
callback_chain.callback!(app)
hello_slider.callback!(app)
multi_inputs.callback!(app)
multi_outputs.callback!(app)
simple_callback.callback!(app)

app.layout = html_div() do
    Header("Basic Dash Callbacks"),
    html_blockquote(dcc_markdown("This is the 3rd chapter of the [Dash Tutorial](/).
    The previous chapter covered the Dash app [layout](/getting-started) and the [next chapter](/state) covers an additional concept of callbacks
    known as `state`.")),
    dcc_markdown("
    In the [previous chapter](/getting-started) on the app `layout` we learned that the
    `app.layout` describes what the app looks like and is a hierarchical tree of components.
    The `DashHtmlComponents` package provides classes for all of the HTML tags, and the keyword
    arguments describe the HTML attributes like `style`, `className`, and `id`. The `DashCoreComponents`
    package generates higher level components like controls and graphs.

    This chapter describes how to make your Dash apps interactive.

    Let's get started with a simple example.
    "),

    simple_callback.source_code,
    Example(simple_callback.layout),

    dcc_markdown("""
    Try typing in the text box. The `children` property of the output
    component updates right away. Let's break down what's happening here:

    1. The "inputs" and "outputs" of our application interface are described
    declaratively by the `callback!` function definition.

    2. In Dash, the inputs and outputs of our application are simply the
    properties of a particular component. In this example, our input is the `value`
    property of of the component with the ID `input`. Our output is the `children` property
    of the compnoent with the ID `output`.

    3. Whenever an input property changes, the callback function will be executed automatically.
    Dash provides the callback function with the new value of the input property as an input argument
    and Dash updates the property of the output component with whatever was returned by the callback
    function.

    4. Don't confuse the `Dash.Input` object with the `DashCoreComponents.Input` object. The former is just used to declare
    inputs of callback functions while the latter is an UI component which is used to render HTML input elements.

    5. Notice how we don't set a value for the `children` property of the `input` component in the `layout`. When
    the Dash app starts, it automatically calls all the callbacks with the initial values of the input componets in
    order to populate the initial state of the output components. In this example, if you specified something like
    `html_div(id=\"output\", children=\"hello world\")`, it would get overwritten when the app starts by what is
    returned by the callback function.

    It's sort of like programming with Microsoft Excel: whenever an input cell changes, all of the cells that
    depend on that cell will get updated automatically. This is called \"Reactive Programming\".

    Remember how every component was described entirely through its set of keyword arguments? Those properties are important now.
    With Dash's interactivity, we can dynamically update any of those properties through a callback function.
    Frequently we'll update the `children` of a component to display new text or the `figure` of a `Graph` component to
    display new data, but we could also update the `style` of a component or even the available `options` of a `Dropdown`
    component!
    """),
    html_hr(),
    dcc_markdown("""
    Let's take a look at another example where a `Slider` updates a `Graph`.

    ## Dash App Layout With Figure and Slider
    """),
    hello_slider.source_code,

    Example(hello_slider.layout),

    dcc_markdown("""
    In this example, the `value` property of the `Slider` component is the input of the app and the output of the app is the `figure` property of the `Graph` component. When the `value` of the `Slider` changes, Dash calls the callback function with new value. The function filters the dataframe with this new value, constructs a `figure` object, and returns it to the Dash app.

    There are a few nice patterns in this example:

    1. We're using the `CSV` and `DataFrames` libraries to import and filter datasets in memory.

    2. We load our dataframe at the start of the app: `df = CSV.read("...")`. This dataframe `df` is in the global state of the app and can be read inside callback functions.

    3.Loading data into memory can be expensive. By loading data at the start of the app instead of inside callback functions, we ensure that this operation is only done with the app server starts. When a user visits the app or interacts with the app, the data (the `df`) is already in memory. If possible, expensive initialization (like downloading or querying data) should be done in the global scope of the app instead of within callback functions.

    4, The callback does not modify the original data, it just creates copies of the dataframe by filtering. This is important: *your callbacks should never mutate variables outside of their scope*. If your callbacks modify global state, then one user's session might affect the next user's session and when the app is deployed on a server with multiple processes or threads, those modifications will not be shared across sessions.

    5. We are turning on transitions with `layout.transition` to give an idea of how the dataset evolves with time: transitions allow the chart to update from one state to the next smoothly, as if it were animated.
    """),
    Header("Dash App With Multiple Inputs"),
    dcc_markdown("""
    In Dash, any `Output` can have multiple `Input` components.
    Here's a simple example that binds five Inputs (the `value` property of
    2 `Dropdown` components, 2 `RadioItems` components, and
    1 `Slider` component) to 1 Output component (the `figure` property of the `Graph` component).
    """),
    multi_inputs.source_code,

    Example(multi_inputs.layout),
    dcc_markdown("""

    In this example, the callback function gets called whenever the `value`
    property of the `Dropdown`, `Slider`, or `RadioItems` components change.

    The input arguments of the callback function are the new or current value of
    each of the `Input` properties, in the order that they were specified.

    Even though only a single `Input` changes at a time (a user can only
    change the value of a single `Dropdown` at a given moment), Dash collects
    the current state of all of the specified `Input` properties and passes them
    into your callback function for you. Your callback functions are always
    guaranteed to have be passed the representative state of the app.

    Let's extend our example to include multiple outputs.
    """),
    Header("Dash App With Multiple Outputs"),
    dcc_markdown("""
    So far all the callbacks we've written only update a single Output property. We can also
    update several at once: put all the properties you want to update as inputs to the callback,
    and return that many items from the callback. This is particularly nice if two outputs depend
    on the same computationaly intense itermediate result, such as a slow database query.
    """),
    multi_outputs.source_code,

    Example(multi_outputs.layout),
    dcc_markdown("""

    A word of caution: it's not always a good idea to combine Outputs, even if you can:
    - If the Outputs depend on some but not all of the same inputs, keeping them separate
    can avoid unnecessary updates.
    - If they have the same inputs but do independent computations with these same inputs,
    keeping the callbacks separate can allow them to run in parallel.
    """),
    Header("Dash App With Chained Callbacks"),
    dcc_markdown("""
    You can also chain outputs and inputs together: the output of one callback function could
    be the input of another callback function.

    This pattern can be used to create dynamic UIs where one input component updates the available
    options of the next input component. Here's a simple example:
    """),
    callback_chain.source_code,

    Example(callback_chain.layout),

    dcc_markdown("""

    The first callback updates the available options in the second
    `RadioItems` component.

    The second callback sets an initial value when the `options` property
    changes: it sets it to the first value in that `options` array.

    The final callback displays the selected `value` of each component. If you
    change the `value` of the countries `RadioItems` component, Dash will wait until
    the `value` of the cities component is updated before calling the final callback.
    This prevents your callbacks from being called with inconsistent state like with
    "America" and "Montreal".
    """),
    Header("Dash Apps With State"),
    dcc_markdown("""
    In some cases, you might have a "form"-type pattern in your appliication.
    In such a situation, you might want to read the value of the input component,
    but only when the user is finished entering all of his or her information in
    the form.

    Attaching a callback to the input values directly can look like this:
    """),
    basic_input.source_code,

    Example(basic_input.layout),
    dcc_markdown("""

    In this example, the callback function is fired whenever any of the attributes
    described by the `Input` change. Try it for yourself by entering data in the inputs
    above.

    `State` allows you pass along extra values without firing the callbacks. Here's the same
    example as above but with the `Input` as `State` and a button as an `Input`.
    """),
    basic_state.source_code,

    Example(basic_state.layout),
    dcc_markdown("""

    In this example, changing text in the `Input` boxes won't fire the callback but clicking on the button will.
    The current values of the `Input` values are still passed into the callback even though they don't trigger the
    callback function itself.

    Note that we're triggering the callback by listening to the `n_clicks` property of the `Button` component. `n_clicks`
    is a property that gets incremented every time the component has been clicked on. It is available in
    every component in the `DashHtmlComponents` components library.
    """),
    Header("Summary"),
    dcc_markdown("""
    We've covered the fundamentals of callback functions in Dash. Dash apps are built on a set of
    simple but powerful principles: declarative UIs that are customizable through
    reaction and functional Julia callbacks. Every element attribute of the
    declarative components can be updated through a callback and a subset of the
    attributes, like the `value` property of a `Dropdown`, are editable by the user
    in the interface.
    """),
    dcc_markdown("""
    The next part of the Dash  tutorial covers interactive graphing.

    [Dash Tutorial Part 4: Interactive Graphing](/interactive-graphing)
    """),
    html_hr(),
    dcc_markdown("[Back to the table of contents](/)")

end

end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs