text:
#= Simple classifiers
This is where we stand using the default parameters for 5 basic classifiers
Start the Julia REPL, drag to highlight a section, and use ctrl-enter to run it
This only needs one file, audio_recoded_combined_binary.csv, in the base folder
-Matei
=#

using Statistics

# Logistic Regression ##########################################################
# https://coinslab.github.io/CognitiveModelingLabWebsite/page/LogisticRegression/
#
# Results (test data):
#               precision    recall  f1-score   support
#            0       0.84      1.00      0.91     22560
#            1       0.00      0.00      0.00      4332
#     accuracy                           0.84     26892
#    macro avg       0.42      0.50      0.46     26892
# weighted avg       0.70      0.84      0.77     26892
#
# Note that training data looks the same, it's not overfitting/memorizing
# it's just failing completely and guessing false on all inputs
#
# Cross-validation:
#   0.8362637362637363
#   0.8364468864468865
#   0.8362637362637363
#   0.8361571572488323
#   0.8363403242055133
#
#   83.6% accuracy, because 83.6% of the data set is non-REM and it's guessing
#   non-REM for every single input, so it's guessing 83.6% of them right and
#   all of the rest wrong.
#
# Confusion Matrix:
#   0    45661    16
#   1     8921     0
#            1     0
#
# This is fantastically useless: it's predicting 0 (false) on every single item
# Unless something is fixable here, solving this problem is probably impossible
# with logistic regression and it's best to move on
# TODO: Someone better with the theory can probably figure out why
################################################################################

using DataFrames, CSV, ScikitLearn, PyPlot, Statistics

data = CSV.File("audio_recoded_combined_binary.csv") |> DataFrame

X = convert(Array, data[!,Not(:Class)])
y = convert(Array, data[!,:Class]) # :Class is our target variable

@sk_import model_selection: train_test_split
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.33, random_state=42)

@sk_import linear_model: LogisticRegression
simplelogistic =LogisticRegression(penalty=:none)

fit!(simplelogistic,X_train,y_train)

@sk_import metrics: classification_report
# Training performance
y_pred = predict(simplelogistic,X_train)
print(classification_report(y_train,y_pred))
# Testing performance
y_pred = predict(simplelogistic,X_test)
print(classification_report(y_test,y_pred))

# Cross validation
@sk_import model_selection: cross_val_score
@time cva = (cross_val_score(LogisticRegression(penalty=:none ), X_train, y_train, cv=10))
mean(cva)
@time cvp = cross_val_score(LogisticRegression(penalty=:none ), X_train, y_train, cv=10, scoring="precision")
mean(cvp)
@time cvr = cross_val_score(LogisticRegression(penalty=:none ), X_train, y_train, cv=10, scoring="recall")
mean(cvr)

# Confusion matrix
@sk_import metrics: plot_confusion_matrix
plot_confusion_matrix(simplelogistic,X_train,y_train)
PyPlot.gcf()
plot_confusion_matrix(simplelogistic,X_test,y_test)
PyPlot.gcf()

# Neural network ###############################################################
# https://coinslab.github.io/CognitiveModelingLabWebsite/page/NeuralNetwork/
#
# Classification report (test set) 1:
#
#           precision    recall  f1-score   support
#    0           0.99      0.98      0.98     22560
#    1           0.91      0.93      0.92      4332
#
# accuracy                           0.97     26892
# macro avg      0.95      0.95      0.95     26892
# weighted avg   0.97      0.97      0.97     26892
#
# Cross Validation 1:
#   0.9662087912087912
#   0.9668498168498169
#   0.8365384615384616    -- This set failed to train completely! Why?
#   0.9733492078029123
#   0.968220533015844
#
# Cross Validation 2:
#   0.9686813186813187
#   0.9574175824175825
#   0.9741758241758242
#   0.9707848704093781
#   0.9638245260554996
#
# Cross Validation 3 (10 k-sets this time):
#   0.9681318681318681
#   0.9686813186813187
#   0.9467032967032967
#   0.9701465201465201
#   0.9635531135531136
#   0.9622710622710623
#   0.9672161172161172
#   0.9683150183150183
#   0.9690419490749221
#   0.9736215424070342
#
# Confusion Matrix 1:
#   0    44353   1324
#   1      807   8114
#            1      0
#
# Looks great except for the bad k-set in the first cross-validation.
# Other than that accuracy is a solid 97% and recall is from 83% to 93% so it's
# not at the price of discarding too many of the REM inputs.
#
# Cross-validation k=10: 182.194503 seconds, 166.427203 seconds, etc
#
################################################################################

using DataFrames, CSV, ScikitLearn, PyPlot, Statistics

data = CSV.File("audio_recoded_combined_binary.csv") |> DataFrame

X = convert(Array, data[!,Not(:Class)])
y = convert(Array, data[!,:Class]) # :Class is our target variable

@sk_import model_selection: train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

@sk_import neural_network: MLPClassifier
mlp_6layer = MLPClassifier(hidden_layer_sizes=(30, 50, 60, 10, 10, 10))

fit!(mlp_6layer,X_train,y_train)

@sk_import metrics: classification_report
# Training performance
y_pred = predict(mlp_6layer,X_train)
print(classification_report(y_train,y_pred))
# Testing performance
y_pred = predict(mlp_6layer,X_test)
print(classification_report(y_test,y_pred))

# Cross validation
@sk_import model_selection: cross_val_score
@time cva = cross_val_score(MLPClassifier(hidden_layer_sizes=(30, 50, 60, 10, 10, 10)), X_train, y_train, cv=10)
mean(cva)
@time cvp = cross_val_score(MLPClassifier(hidden_layer_sizes=(30, 50, 60, 10, 10, 10)), X_train, y_train, cv=10, scoring="precision")
mean(cvp)
@time cvr = cross_val_score(MLPClassifier(hidden_layer_sizes=(30, 50, 60, 10, 10, 10)), X_train, y_train, cv=10, scoring="recall")
mean(cvr)

# Confusion matrix
@sk_import metrics: plot_confusion_matrix
plot_confusion_matrix(mlp_6layer,X_train,y_train)
PyPlot.gcf()

@sk_import metrics: plot_confusion_matrix
plot_confusion_matrix(mlp_6layer,X_test,y_test)
PyPlot.gcf()

# Decision Tree ################################################################
# https://coinslab.github.io/CognitiveModelingLabWebsite/page/tree/
#
# Classification report (test set):
#               precision    recall  f1-score   support
#            0       0.98      0.98      0.98     22560
#            1       0.90      0.90      0.90      4332
#     accuracy                           0.97     26892
#    macro avg       0.94      0.94      0.94     26892
# weighted avg       0.97      0.97      0.97     26892
#
# Cross-Validation:
#   0.9652014652014652
#   0.9666666666666667
#   0.967948717948718
#   0.9664835164835165
#   0.969047619047619
#   0.9659340659340659
#   0.9661172161172161
#   0.9677655677655678
#   0.9673932954753618
#   0.9629968858765342
#
# Results are great, comparable with the neural network
# However this runs MUCH faster (2 orders of magnitude vs neural network)
# k=10 cross validation runs in  1.579159 seconds
#
################################################################################

using DataFrames, CSV, ScikitLearn, PyPlot, Statistics

data = CSV.File("audio_recoded_combined_binary.csv") |> DataFrame
X = convert(Array, data[!,Not(:Class)])
y = convert(Array, data[!,:Class]) # :Class is our target variable

@sk_import model_selection: train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

@sk_import tree: DecisionTreeClassifier
tree = DecisionTreeClassifier()
fit!(tree,X_train,y_train)

@sk_import metrics: classification_report
y_pred = predict(tree,X_train)
print(classification_report(y_train,y_pred))
y_pred = predict(tree,X_test)
print(classification_report(y_test,y_pred))

@sk_import metrics: plot_confusion_matrix
plot_confusion_matrix(tree,X_train,y_train)
PyPlot.gcf()

@sk_import metrics: plot_confusion_matrix
plot_confusion_matrix(tree,X_test,y_test)
PyPlot.gcf()

@sk_import model_selection: cross_val_score
@time cva = cross_val_score(DecisionTreeClassifier(), X_train, y_train, cv=10)
mean(cva)
@time cvp = cross_val_score(DecisionTreeClassifier(), X_train, y_train, cv=10, scoring="precision")
mean(cvp)
@time cvr = cross_val_score(DecisionTreeClassifier(), X_train, y_train, cv=10, scoring="recall")
mean(cvr)

# Naive Bayes ##################################################################
# https://coinslab.github.io/CognitiveModelingLabWebsite/page/nbclassifier/
#
#               precision    recall  f1-score   support
#            0       0.87      1.00      0.93     22560
#            1       0.99      0.22      0.35      4332
#     accuracy                           0.87     26892
#    macro avg       0.93      0.61      0.64     26892
# weighted avg       0.89      0.87      0.84     26892
#
# Cross-Validation:
#   0.8697802197802198
#   0.8703296703296703
#   0.8725274725274725
#   0.8716117216117216
#   0.8745421245421245
#   0.8732600732600733
#   0.8710622710622711
#   0.8708791208791209
#   0.8728704891005679
#   0.8754350613665507
#
# This doesn't work well. It's just guessing "no" on almost all the inputs and
# getting absolutely abysmal recall (22%).
# TODO: Find out why.
# On the upside, it runs fast.  k=10 cross validation runs in 0.104042 seconds
#
################################################################################

using DataFrames, CSV, ScikitLearn, PyPlot, Statistics

data = CSV.File("audio_recoded_combined_binary.csv") |> DataFrame
X = convert(Array, data[!,Not(:Class)])
y = convert(Array, data[!,:Class]) # :Class is our target variable

@sk_import model_selection: train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

@sk_import naive_bayes: GaussianNB
gnb = GaussianNB()
fit!(gnb,X_train,y_train)

@sk_import metrics: classification_report

y_pred = predict(gnb,X_train)
print(classification_report(y_train,y_pred))
y_pred = predict(gnb,X_test)
print(classification_report(y_test,y_pred))

@sk_import metrics: plot_confusion_matrix
plot_confusion_matrix(gnb,X_train,y_train)
PyPlot.gcf()
plot_confusion_matrix(gnb,X_test,y_test)
PyPlot.gcf()

@sk_import model_selection: cross_val_score
@time cva = cross_val_score(GaussianNB(), X_train, y_train, cv=10)
mean(cva)
@time cvp = cross_val_score(GaussianNB(), X_train, y_train, cv=10, scoring="precision")
mean(cvp)
@time cvr = cross_val_score(GaussianNB(), X_train, y_train, cv=10, scoring="recall")
mean(cvr)

# Random Forest ################################################################
# https://coinslab.github.io/CognitiveModelingLabWebsite/page/nbclassifier/
#
# Classification report (test set):
#               precision    recall  f1-score   support
#            0       0.98      0.99      0.99     22560
#            1       0.96      0.92      0.94      4332
#     accuracy                           0.98     26892
#    macro avg       0.97      0.96      0.96     26892
# weighted avg       0.98      0.98      0.98     26892
#
# Cross-Validation:
#   0.9782051282051282
#   0.9782051282051282
#   0.9793040293040293
#   0.9774725274725274
#   0.978021978021978
#   0.9798534798534798
#   0.9787545787545787
#   0.9805860805860805
#   0.9803993405385601
#   0.9778347682725774
#
# This is the best so far.
# Not as fast as Naive Bayes or Tree but 3 times as fast as the neural net.
# k=10 cross validation runs in 54.345652 seconds
#
################################################################################

using DataFrames, CSV, ScikitLearn, PyPlot, Statistics

data = CSV.File("audio_recoded_combined_binary.csv") |> DataFrame
X = convert(Array, data[!,Not(:Class)])
y = convert(Array, data[!,:Class]) # :Class is our target variable

@sk_import model_selection: train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

@sk_import ensemble: RandomForestClassifier
rfc = RandomForestClassifier()
fit!(rfc,X_train,y_train)

@sk_import metrics: classification_report
y_pred = predict(rfc,X_train)
print(classification_report(y_train,y_pred))
y_pred = predict(rfc,X_test)
print(classification_report(y_test,y_pred))

@sk_import metrics: plot_confusion_matrix
plot_confusion_matrix(rfc,X_train,y_train)
PyPlot.gcf()
plot_confusion_matrix(rfc,X_test,y_test)
PyPlot.gcf()

@sk_import model_selection: cross_val_score
@time cva = cross_val_score(RandomForestClassifier(), X_train, y_train, cv=10)
mean(cva)
@time cvp = cross_val_score(RandomForestClassifier(), X_train, y_train, cv=10, scoring="precision")
mean(cvp)
@time cvr = cross_val_score(RandomForestClassifier(), X_train, y_train, cv=10, scoring="recall")
mean(cvr)


text:
@safetestset "AppServer functionality" begin

  @safetestset "Start/stop servers" begin
    using Genie
    using Genie.AppServer

    servers = Genie.AppServer.startup(; open_browser = false)
    @test servers.webserver.state == :runnable
    @test Genie.AppServer.SERVERS.webserver.state == :runnable

    servers = Genie.AppServer.down()
    sleep(2)
    @test servers.webserver.state == :done
    @test Genie.AppServer.SERVERS.webserver.state == :done

    servers = Genie.AppServer.startup(; open_browser = false)
    Genie.AppServer.down(; webserver = false)
    sleep(2)
    @test servers.webserver.state == :runnable
    @test Genie.AppServer.SERVERS.webserver.state == :runnable

    servers = Genie.AppServer.down(; webserver = true)
    sleep(2)
    @test servers.webserver.state == :done
    @test Genie.AppServer.SERVERS.webserver.state == :done
  end;

  @safetestset "Update config when custom startup args" begin
    using Genie
    using Genie.AppServer

    port = Genie.config.server_port
    ws_port = Genie.config.websockets_port

    Genie.AppServer.up(port+1_000; ws_port = ws_port+1_000, open_browser = false)

    @test Genie.config.server_port == port+1_000
    @test Genie.config.websockets_port == ws_port+1_000

    Genie.config.server_port = port
    Genie.config.websockets_port = ws_port

    Genie.AppServer.down()
  end;

end;

ext: jl


text:
include(joinpath(@__DIR__, "extension.jl"))
include(joinpath(@__DIR__, "node_bb.jl"))
include(joinpath(@__DIR__, "variable_info.jl"))
include(joinpath(@__DIR__, "functions/functions.jl"))
include(joinpath(@__DIR__, "problem_parsed.jl"))
include(joinpath(@__DIR__, "global_options.jl"))
include(joinpath(@__DIR__, "global_optimizer.jl"))
include(joinpath(@__DIR__, "optimizer.jl"))


text:
################################################################################
### FUNCTIONS FOR SOLVING COORDINATED SUPPLY CHAIN MODELS WITH JuMP
"""
    Output, Stats = OptimizeSSCase(A,N,P,D,S,T,L,Sets,Pars; CaseDataDirectory=pwd(), PrintModel=false, WriteReport=true, PrintOutput=true, ModelOutputFileName="_Model.txt", SolutionOutputFileName="_SolutionData.txt", PrintSpacer="*"^50)

# Arguments

A - arc structure from ImportFunctions
N - node structure from ImportFunctions
P - product structure from ImportFunctions
D - demand structure from ImportFunctions
S - supply structure from ImportFunctions
T - technology structure from ImportFunctions
L - technology site structure from ImportFunctions
Sets - structure containing additional sets from ImportFunctions
Pars - structure constaning additional parameters from ImportFunctions
CaseDataDirectory=pwd(): (optional keyword arg) file directory for supply chain case study data; defaults to current Julia directory if not specified
PrintModel=false: (optional keyword arg) Print the model to the REPL
WriteReport=true: (optional keyword arg) Write report to text
PrintOutput=true: (optional keyword arg) Print outputs to Julia REPL
ModelOutputFileName="_Model.txt": (optional keyword arg) Change the output text file name for the model
SolutionOutputFileName="_SolutionData.txt": (optional keyword arg) Change the output text file name for the solution report
PrintSpacer="*"^50: (optional keyword arg): Change the characters used to space out the report

# Returns

Output - Data strucuture with outputs from JuMP model

Stats - Data structure with statistics from JuMP model
```
Optimizes a steady state supply chain problem using data loaded by ImportFunctions()

Optional arguments (see arguments list) allow you to customize output behaviour
"""
function OptimizeSSCase(A,N,P,D,S,T,L,Sets,Pars; CaseDataDirectory=pwd(), PrintModel=false, WriteReport=true, PrintOutput=true, ModelOutputFileName="_Model.txt", SolutionOutputFileName="_SolutionData.txt", PrintSpacer="*"^50)
    ################################################################################
    ### MODEL STATEMENT
    MOD = JuMP.Model(optimizer_with_attributes(Clp.Optimizer))

    ################################################################################
    ### VARIABLES
    @variable(MOD, 0 <= s[n=N.ID,p=P.ID] <= Pars.sMAX[n,p]) # nodal supply allocation
    @variable(MOD, 0 <= d[n=N.ID,p=P.ID] <= Pars.dMAX[n,p]) # nodal demand allocation
    @variable(MOD, 0 <= sl[i=S.ID] <= S.cap[i]) # individual supplier allocation
    @variable(MOD, 0 <= dl[j=D.ID] <= D.cap[j]) # individual consumer allocation
    @variable(MOD, 0 <= f[a=A.ID,p=P.ID] <= Pars.fMAX[a,p]) # transport allocation (arc-indexed)
    @variable(MOD, 0 <= ξcon[t=T.ID,n=N.ID,q=P.ID] <= Pars.ξconMAX[t,n,q]) # tech. consumption, standard: products q in P committed
    @variable(MOD, 0 <= ξgen[t=T.ID,n=N.ID,p=P.ID] <= Pars.ξgenMAX[t,n,p]) # tech. generation, standard: products p in P produced

    ################################################################################
    ### EQUATIONS
    # nodal supply and nodal demand total equal total of individual supplies and individual demands
    @constraint(MOD, SupplyBalance[n=N.ID,p=P.ID], s[n,p] == sum(sl[i] for i in S.ID if S.node[i] == n && S.prod[i] == p))
    @constraint(MOD, DemandBalance[n=N.ID,p=P.ID], d[n,p] == sum(dl[j] for j in D.ID if D.node[j] == n && D.prod[j] == p))

    # System mass balance
    @constraint(MOD, Balance[n=N.ID,p=P.ID],
        s[n,p] + sum(f[a,p] for a in Sets.Ain[n]) + sum(ξgen[t,n,p] for t in Sets.NPt[n,p]) ==
        d[n,p] + sum(f[a,p] for a in Sets.Aout[n]) + sum(ξcon[t,n,p] for t in Sets.NQt[n,p]))

    # Conversion relationships (yield-based)
    @constraint(MOD, Conversion[n=N.ID,(t,p,q) in Sets.TPQ], ξgen[t,n,p] == Pars.α[t,p,q]*ξcon[t,n,q])

    ################################################################################
    ### OBJECTIVE
    # Cost expressions
    demand_revenue = @expression(MOD, sum(D.bid[j]*dl[j] for j in D.ID))
    supply_revenue = @expression(MOD, sum(S.bid[i]*sl[i] for i in S.ID))
    transport_cost = @expression(MOD, sum(P.transport_cost[p]*A.len[a]*f[a,p] for a in A.ID,p in P.ID))
    operating_cost = @expression(MOD, sum(ξcon[t,n,q]*T.bid[t] for t in T.ID, n in N.ID, q in P.ID if (Sets.NQT[n,q,t] && T.InputRef[t] == q)))

    # Full objective
    @objective(MOD, Max, demand_revenue - supply_revenue - transport_cost - operating_cost)

    ################################################################################
    ### DISPLAY MODEL FORMULATION
    filename = open(CaseDataDirectory*"/"*ModelOutputFileName,"w")
    print(filename, MOD)
    close(filename)
    if PrintModel
        print(MOD)
    end

    ################################################################################
    ### SOLVE AND DATA RETRIEVAL
    # Display statistics
    NumVars = length(all_variables(MOD))
    TotalIneqCons = num_constraints(MOD,AffExpr, MOI.LessThan{Float64})+num_constraints(MOD,AffExpr, MOI.GreaterThan{Float64})+num_constraints(MOD,VariableRef, MOI.LessThan{Float64})+num_constraints(MOD,VariableRef, MOI.GreaterThan{Float64})
    TotalEqCons = num_constraints(MOD,VariableRef, MOI.EqualTo{Float64})+num_constraints(MOD,AffExpr, MOI.EqualTo{Float64})
    NumVarBounds = num_constraints(MOD,VariableRef, MOI.LessThan{Float64})+num_constraints(MOD,VariableRef, MOI.GreaterThan{Float64})
    ModelIneqCons = num_constraints(MOD,AffExpr, MOI.LessThan{Float64})+num_constraints(MOD,AffExpr, MOI.GreaterThan{Float64})
    ModelEqCons = num_constraints(MOD,AffExpr, MOI.EqualTo{Float64})

    println(PrintSpacer*"\nModel statistics:")
    println("Variables: "*string(NumVars))
    println("Total inequality constraints: "*string(TotalIneqCons))
    println("Total equality constraints: "*string(TotalEqCons))
    println("Variable bounds: "*string(NumVarBounds))
    println("Model inequality constraints: "*string(ModelIneqCons))
    println("Model equality constraints: "*string(ModelEqCons))
    println(PrintSpacer)

    # Solve
    println("Solving original problem...")
    JuMP.optimize!(MOD)
    println("Primal status: "*string(primal_status(MOD)))
    println("Dual status: "*string(dual_status(MOD)))
    z_out = JuMP.objective_value(MOD)
    s_out = JuMP.value.(s)
    d_out = JuMP.value.(d)
    sl_out = JuMP.value.(sl)
    dl_out = JuMP.value.(dl)
    f_out = JuMP.value.(f)
    ξcon_out = JuMP.value.(ξcon)
    ξgen_out = JuMP.value.(ξgen)
    π_out = JuMP.dual.(Balance)

    ################################################################################
    ### PROFIT CALCULATIONS

    PhiSupply = Dict()
    for i in S.ID
        PhiSupply[i] = (π_out[S.node[i],S.prod[i]] - S.bid[i])*sl_out[i]
    end

    PhiDemand = Dict()
    for j in D.ID
        PhiDemand[j] = (D.bid[j] - π_out[D.node[j],D.prod[j]])*dl_out[j]
    end

    PhiTransport = Dict()
    π_transport = Dict()
    for a in A.ID
        for p in P.ID
            # transport prices (notation: s -> r for p)
            π_transport[a,p] = π_out[A.node_r[a],p] - π_out[A.node_s[a],p]
            # transport profits
            PhiTransport[a,p] = (π_transport[a,p] - P.transport_cost[p]*A.len[a])*f_out[a,p]
        end
    end

    PhiTech = DictInit([T.ID,N.ID],0)
    π_tech = DictInit([T.ID,N.ID],0)
    for t in T.ID
        for n in N.ID
            # technology prices
            π_tech[t,n] = sum(π_out[n,p]*T.OutputYields[t,p] for p in T.Outputs[t]) - sum(π_out[n,q]*T.InputYields[t,q] for q in T.Inputs[t])
        end
    end
    for ts in L.ID
        n = L.node[ts]
        t = L.tech[ts]
        p = T.InputRef[t]
        PhiTech[t,n] = (π_tech[t,n] - T.InputYields[t,p]*T.bid[t])*ξcon_out[t,n,p]
    end

    ################################################################################
    ### WRITE OUTPUT TO FILE AND DISPLAY TO REPL
    if WriteReport
        filename = open(CaseDataDirectory*"/"*SolutionOutputFileName,"w")
        print(filename, PrintSpacer*"\nModel statistics:")
        print(filename, "\nVariables: "*string(NumVars))
        print(filename, "\nTotal inequality constraints: "*string(TotalIneqCons))
        print(filename, "\nTotal equality constraints: "*string(TotalEqCons))
        print(filename, "\nVariable bounds: "*string(NumVarBounds))
        print(filename, "\nModel inequality constraints: "*string(ModelIneqCons))
        print(filename, "\nModel equality constraints: "*string(ModelEqCons))
        print(filename, "\n"*PrintSpacer*"\nObjective value: "*string(z_out))
        FilePrint(s_out,[N.ID,P.ID],filename,DataName="Supply values:",VarName="s")
        FilePrint(d_out,[N.ID,P.ID],filename,DataName="Demand values:",VarName="d")
        FilePrint(f_out,[A.ID,P.ID],filename,DataName="Transport values:",VarName="f")
        FilePrint(ξcon_out,[T.ID,N.ID,P.ID],filename,DataName="Consumption values:",VarName="ξcon")
        FilePrint(ξgen_out,[T.ID,N.ID,P.ID],filename,DataName="Generation values:",VarName="ξgen")
        FilePrint(π_out,[N.ID,P.ID],filename,DataName="Nodal clearing prices:",VarName="π")
        FilePrint(π_transport,[A.ID,P.ID],filename,DataName="Transport clearing prices:",VarName="π_f")
        FilePrint(π_tech,[T.ID,N.ID],filename,DataName="Technology clearing prices:",VarName="π_t")
        FilePrint(PhiDemand,[D.ID],filename,DataName="Demand profits:",VarName="Φ_d")
        FilePrint(PhiSupply,[S.ID],filename,DataName="Supply profits:",VarName="Φ_s")
        FilePrint(PhiTransport,[A.ID,P.ID],filename,DataName="Transport profits:",VarName="Φ_f")
        FilePrint(PhiTech,[T.ID,N.ID],filename,DataName="Technology profits:",VarName="Φ_ξ")
        close(filename)
    end

    # And print to REPL
    if PrintOutput
        println(PrintSpacer)
        println("Objective value: ", z_out)
        # NOTE to display all values, replace Nonzeros() with DictInit([A,B,...Z],true), or (for specific values) with a TruthTable
        PrettyPrint(s_out, [N.ID,P.ID], Nonzeros(s_out, [N.ID,P.ID]), DataName="Supply values:", VarName="s")
        PrettyPrint(d_out, [N.ID,P.ID], Nonzeros(d_out, [N.ID,P.ID]), DataName = "Demand values:", VarName="d")
        PrettyPrint(f_out, [A.ID,P.ID], Nonzeros(f_out, [A.ID,P.ID]), DataName = "Transport values:", VarName="f")
        PrettyPrint(ξcon_out, [T.ID,N.ID,P.ID], Nonzeros(ξcon_out, [T.ID,N.ID,P.ID]), DataName = "Consumption values:", VarName="ξcon")
        PrettyPrint(ξgen_out, [T.ID,N.ID,P.ID], Nonzeros(ξgen_out, [T.ID,N.ID,P.ID]), DataName = "Generation values:", VarName="ξgen")
        PrettyPrint(π_out, [N.ID,P.ID], Nonzeros(π_out, [N.ID,P.ID]), DataName = "Nodal clearing prices:", VarName="π")
        PrettyPrint(π_transport, [A.ID,P.ID], Nonzeros(π_transport, [A.ID,P.ID]), DataName = "Transport clearing prices:", VarName="π_f")
        PrettyPrint(π_tech, [T.ID,N.ID], Nonzeros(π_tech, [T.ID,N.ID]), DataName = "Technology clearing prices:", VarName="π_t")
        PrettyPrint(PhiDemand, [D.ID], Nonzeros(PhiDemand, [D.ID]), DataName = "Demand profits:", VarName="Φ_d")
        PrettyPrint(PhiSupply, [S.ID], Nonzeros(PhiSupply, [S.ID]), DataName = "Supply profits:", VarName="Φ_s")
        PrettyPrint(PhiTransport, [A.ID,P.ID], Nonzeros(PhiTransport, [A.ID,P.ID]), DataName = "Transport profits:", VarName="Φ_f")
        PrettyPrint(PhiTech, [T.ID,N.ID], Nonzeros(PhiTech, [T.ID,N.ID]), DataName = "Technology profits:", VarName="Φ_ξ")
    end

    ################################################################################
    ### PACKAGE INTO RETURN STRUCTURE
    Output = OutputStruct(z_out,sl_out,dl_out,s_out,d_out,f_out,ξcon_out,ξgen_out,π_out,π_transport,π_tech,PhiDemand,PhiSupply,PhiTransport,PhiTech);
    Stats = StatStruct(NumVars,TotalIneqCons,TotalEqCons,NumVarBounds,ModelIneqCons,ModelEqCons)

    ################################################################################
    ### UPDATE USER
    println(PrintSpacer,"\nModel Complete\n",PrintSpacer)

    ################################################################################
    ### RETURN
    return(Output,Stats);
    # e.g.: Output = OptimizeSSCase(A,N,P,D,S,T,L,Sets,Pars,CaseDataDirectory="/Users/ptominac/Documents/Code/CoordinatedSupplyChains/TestCases/TestV36");
end


text:
using Bonobo
using Documenter

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

makedocs(;
    modules=[Bonobo],
    authors="Ole Kroeger <o.kroeger@opensourc.es> and contributors",
    repo="https://github.com/Wikunia/Bonobo.jl/blob/{commit}{path}#{line}",
    sitename="Bonobo.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Wikunia.github.io/Bonobo.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "How-To" => "how_to.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/Wikunia/Bonobo.jl",
    devbranch="main",
)


text:
module libcimgui
	using CBinding
	using CImGui_jll

	c`$([
		Main.SYSROOT...,

		"-DCIMGUI_DEFINE_ENUMS_AND_STRUCTS=1",
		"-I$(CImGui_jll.artifact_dir)/include",
		"-L$(dirname(CImGui_jll.libcimgui_path))", "-lcimgui",
	])`

	const c"int8_t" = Int8
	const c"int16_t" = Int16
	const c"int32_t" = Int32
	const c"int64_t" = Int64
	const c"uint8_t" = UInt8
	const c"uint16_t" = UInt16
	const c"uint32_t" = UInt32
	const c"uint64_t" = UInt64
	const c"size_t" = Csize_t

	# define them as Cvoid since they are usually used as opaque
	const c"FILE"    = Cvoid
	const c"va_list" = Cvoid

	c"""
		#include <cimgui.h>
	"""ji
end


text:
module DemoPackageZW

print_greeting(io::IO =stdout) = printstyled(io, "Hello, Pingchuan!"; color=:cyan)

end


text:
function dijkstra(net::AbstractNetwork,
                  src::T,
                  linkcosts::AbstractVector{U}) where {T<:Integer, U<:Real}

    nnodes = numnodes(net)
    costs = fill(typemax(U), nnodes)
    parentvec = zeros(T, nnodes)

    P = PriorityQueue{T,U}()
    P[src] = costs[src] = 0
    parentvec[src] = src

    while !isempty(P)
        u = dequeue!(P)
        for v in outneighbors(net, u)
            lidx = idx(net, u, v)
            alt = costs[u] + linkcosts[lidx]
            if (costs[v] > alt) && (throughflowallowed(net, u) || (u == src))
                P[v] = costs[v] = alt
                parentvec[v] = u
            end
        end
    end

    return (costs = costs, parentvec = parentvec)
end

function dijkstra(net::AbstractNetwork,
                  srcs::AbstractVector{T},
                  linkcosts::AbstractVector{U}) where {T<:Integer, U<:Real}

    nsrcs = length(srcs)
    nnodes = numnodes(net)
    costs, parentmx = zeros(U, nsrcs, nnodes), zeros(T, nsrcs, nnodes)
    @sync @distributed for i in 1:nsrcs
        res = dijkstra(net, srcs[i], linkcosts)
        costs[i,:] = res[1]
        parentmx[i,:] = res[2]
    end

    return (costs = costs, parentmx = parentmx)
end


text:
"""
    struct AffineOperator{A<:AbstractMatrix,B<:AbstractVector} <: NonlinearOperator
      matrix::A
      vector::B
    end
"""
struct AffineOperator{A<:AbstractMatrix,B<:AbstractVector} <: NonlinearOperator
  matrix::A
  vector::B
end

"""
    get_matrix(operator)

Return the matrix corresponding to the assembled left hand side of the operator.
This matrix incorporates all boundary conditions and constraints.
"""
get_matrix(op::AffineOperator) = op.matrix

"""
    get_vector(operator)

Return the vector corresponding to the assembled right hand side of the operator.
This vector includes all boundary conditions and constraints.
"""
get_vector(op::AffineOperator) = op.vector

function residual!(b::AbstractVector,op::AffineOperator,x::AbstractVector)
  mul!(b,op.matrix,x)
  b .-= op.vector
  b
end

function jacobian!(A::AbstractMatrix,op::AffineOperator,x::AbstractVector)
  copy_entries!(A,op.matrix)
  A
end

function jacobian(op::AffineOperator,x::AbstractVector)
  op.matrix
end

function zero_initial_guess(op::AffineOperator)
  x = allocate_in_domain(typeof(op.vector),op.matrix)
  fill!(x,zero(eltype(x)))
  x
end

function allocate_residual(op::AffineOperator,x::AbstractVector)
  allocate_in_range(typeof(op.vector),op.matrix)
end

function allocate_jacobian(op::AffineOperator,x::AbstractVector)
  op.matrix
end

"""
    abstract type LinearSolver <: NonlinearSolver end

- [`symbolic_setup(::LinearSolver,mat::AbstractMatrix)`](@ref)
- [`test_linear_solver`](@ref)

"""
abstract type LinearSolver <: NonlinearSolver end

"""
    symbolic_setup(::LinearSolver,mat::AbstractMatrix) -> SymbolicSetup
"""
function symbolic_setup(::LinearSolver,mat::AbstractMatrix)
  @abstractmethod
end

"""
    solve(ls::LinearSolver,A::AbstractMatrix,b::AbstractVector)
"""
function solve(ls::LinearSolver,A::AbstractMatrix,b::AbstractVector)
  ss = symbolic_setup(ls,A)
  ns = numerical_setup(ss,A)
  x = similar(b)
  solve!(x,ns,b)
  x
end

"""
    solve!(x::AbstractVector,ls::LinearSolver,A::AbstractMatrix,b::AbstractVector)
"""
function solve!(x::AbstractVector,ls::LinearSolver,A::AbstractMatrix,b::AbstractVector)
  ss = symbolic_setup(ls,A)
  ns = numerical_setup(ss,A)
  solve!(x,ns,b)
  x
end

"""
    abstract type SymbolicSetup <: GridapType end

- [`numerical_setup(::SymbolicSetup,mat::AbstractMatrix)`](@ref)
"""
abstract type SymbolicSetup <: GridapType end

"""
    numerical_setup(::SymbolicSetup,mat::AbstractMatrix) -> NumericalSetup
"""
function numerical_setup(::SymbolicSetup,mat::AbstractMatrix)
  @abstractmethod
end

"""
    abstract type NumericalSetup <: GridapType end

- [`numerical_setup!(::NumericalSetup,mat::AbstractMatrix)`](@ref)
- [`solve!(x::AbstractVector,::NumericalSetup,b::AbstractVector)`](@ref)
"""
abstract type NumericalSetup <: GridapType end

"""
    numerical_setup!(::NumericalSetup,mat::AbstractMatrix)
"""
function numerical_setup!(::NumericalSetup,mat::AbstractMatrix)
  @abstractmethod
end

"""
    solve!(x::AbstractVector,::NumericalSetup,b::AbstractVector)
"""
function solve!(x::AbstractVector,::NumericalSetup,b::AbstractVector)
  @abstractmethod
end

"""
    test_linear_solver(
      ls::LinearSolver,
      A::AbstractMatrix,
      b::AbstractVector,
      x::AbstractVector)
"""
function test_linear_solver(
  ls::LinearSolver,
  A::AbstractMatrix,
  b::AbstractVector,
  x::AbstractVector)

  y = solve(ls,A,b)
  @test x ≈ y

  ss = symbolic_setup(ls,A)
  ns = numerical_setup(ss,A)
  numerical_setup!(ns,A)
  solve!(y,ns,b)
  @test x ≈ y

  y .= zero(eltype(y))
  solve!(y,ls,A,b)
  @test x ≈ y

  op = AffineOperator(A,b)
  x0 = copy(x)
  test_nonlinear_solver(ls,op,x0,x)

end

# function solve!(x::AbstractVector,ls::LinearSolver,op::NonlinearOperator,cache)
#   @unreachable "A LinearSolver can only solve an AffineOperator"
# end

struct LinearSolverCache <: GridapType
  A::AbstractMatrix
  b::AbstractVector
  ns::NumericalSetup
end

function solve!(x::AbstractVector,
                ls::LinearSolver,
                op::NonlinearOperator,
                cache::Nothing)
  fill!(x,zero(eltype(x)))
  b = residual(op, x)
  A = jacobian(op, x)
  ss = symbolic_setup(ls, A)
  ns = numerical_setup(ss,A)
  rmul!(b,-1)
  solve!(x,ns,b)
  LinearSolverCache(A,b,ns)
end

function solve!(x::AbstractVector,
                ls::LinearSolver,
                op::NonlinearOperator,
                cache)
  fill!(x,zero(eltype(x)))
  b = cache.b
  A = cache.A
  ns = cache.ns
  residual!(b, op, x)
  numerical_setup!(ns,A)
  rmul!(b,-1)
  solve!(x,ns,b)
  cache
end

function solve!(x::AbstractVector,ls::LinearSolver,op::AffineOperator,cache::Nothing)
  A = op.matrix
  b = op.vector
  ss = symbolic_setup(ls,A)
  ns = numerical_setup(ss,A)
  solve!(x,ns,b)
  ns
end

function solve!(x::AbstractVector,ls::LinearSolver,op::AffineOperator,cache)
  newmatrix = true
  solve!(x,ls,op,cache,newmatrix)
  cache
end

"""
    solve!(
      x::AbstractVector,
      ls::LinearSolver,
      op::AffineOperator,
      cache,
      newmatrix::Bool)
"""
function solve!(
  x::AbstractVector,ls::LinearSolver,op::AffineOperator,cache,newmatrix::Bool)
  A = op.matrix
  b = op.vector
  ns = cache
  if newmatrix
    numerical_setup!(ns,A)
  end
  solve!(x,ns,b)
  cache
end

function solve!(
  x::AbstractVector,ls::LinearSolver,op::AffineOperator,cache::Nothing,newmatrix::Bool)
  solve!(x,ls,op,cache)
end

# Concrete implementations

"""
    struct LUSolver <: LinearSolver end

Wrapper of the LU solver available in julia
"""
struct LUSolver <: LinearSolver end

struct LUSymbolicSetup <: SymbolicSetup end

mutable struct LUNumericalSetup{F} <: NumericalSetup
  factors::F
end

symbolic_setup(::LUSolver,mat::AbstractMatrix) = LUSymbolicSetup()

numerical_setup(::LUSymbolicSetup,mat::AbstractMatrix) = LUNumericalSetup(lu(mat))

function numerical_setup!(ns::LUNumericalSetup, mat::AbstractMatrix)
  fac = lu(mat)
  ns.factors = fac
  ns
end

function numerical_setup!(ns::LUNumericalSetup, mat::SparseMatrixCSC)
  lu!(ns.factors,mat)
  ns
end

function solve!(
  x::AbstractVector,ns::LUNumericalSetup,b::AbstractVector)
  ldiv!(x,ns.factors,b)
  x
end

"""
    struct BackslashSolver <: LinearSolver end

Wrapper of the backslash solver available in julia
This is typically faster than LU for a single solve
"""
struct BackslashSolver <: LinearSolver end

struct BackslashSymbolicSetup <: SymbolicSetup end

mutable struct BackslashNumericalSetup{T<:AbstractMatrix} <: NumericalSetup
  A::T
end

symbolic_setup(::BackslashSolver,mat::AbstractMatrix) = BackslashSymbolicSetup()

numerical_setup(::BackslashSymbolicSetup,mat::AbstractMatrix) = BackslashNumericalSetup(mat)

function numerical_setup!(ns::BackslashNumericalSetup, mat::AbstractMatrix)
  ns.A = mat
end

function solve!(
  x::AbstractVector,ns::BackslashNumericalSetup,b::AbstractVector)
  copy_entries!(x, ns.A\b)
end


text:
const c=299792458.
const mu0=4*pi*1e-7
const mur=1
const eps0=1./(mu0*c^2)
const epsr=1

using PyPlot

function urc(Q,f,f0)
    w=f/f0
    I=1/sqrt(1+Q^2*(w-1/w)^2)
    return I
end

const a=2.9#3.11
const b=3.7#3.62
const d=8.7#4.84
const f0=40e6
const f1=1000e6

#rf=.01
#nf=int(1/log(1+rf)*log(f1/f0))
#f=zeros(nf)
#for i=1:nf
#  f[i]=f0*(1+rf)^i
#end

const df=1e6
const nf=int((f1-f0)/df)+1
const f=linspace(f0,f1,nf)#array([20e6,40e6,45e6,56e6,65e6])#linspace(20e6,.3e9,281)#linspace(20e6,3e8,1121)
const V=a*b*d
const Q4=16*pi^2*V*f.^3/(1*c^3)
const Q1=3*V./(2*1*(2./sqrt(2*pi*f.*mur*mu0*1*1e5)))/200
const Q=1./(1./Q1+1./Q4)

const Npt=150
const x=rand(Npt)*(a-1)+0.5
const y=rand(Npt)*(b-1)+0.5
const z=rand(Npt)*(d-1)+0.5

#n_pos=10
#x_p=linspace(0.5,a-0.5,n_pos)#rand(n_pos)*(a-2)+1
#y_p=linspace(0.5,b-0.5,n_pos)#rand(n_pos)*(b-2)+1
#z_p=linspace(0.5,d-0.5,n_pos)#rand(n_pos)*(d-2)+1

#using Grid
#x_pi=InterpGrid(x_p,BCnil,InterpQuadratic)
#y_pi=InterpGrid(y_p,BCnil,InterpQuadratic)
#z_pi=InterpGrid(z_p,BCnil,InterpQuadratic)

n_steps=360
phi=linspace(2*pi/n_steps,2*pi,n_steps)
n_sources = 6

x_b=zeros(n_sources,n_steps)
y_b=zeros(n_sources,n_steps)
z_b=zeros(n_sources,n_steps)

for i=1:n_sources
  phi0=rand()*2pi
  r0=rand()
  x_b[i,:]=ones(n_steps)*(i*(a-1)/6+0.5)
  y_b[i,:]=1.5+r0*sin(phi+phi0)
  z_b[i,:]=1.5+r0*cos(phi+phi0)
end

#plot3D(x_pi[1:0.1:n_pos],y_pi[1:0.1:n_pos],z_pi[1:0.1:n_pos])

E0=1
H0=1/(c*mu0*mur)

Nmodesx=zeros(nf)
Nx=zeros(nf)

Ex=zeros(Complex128,nf,n_steps,Npt)
for u=1:nf
  n_x=int(f[u]*(1+100/Q[u])*2*a/c)
  n_y=int(f[u]*(1+100/Q[u])*2*b/c)
  n_z=int(f[u]*(1+100/Q[u])*2*d/c)
  M=[]
  N=[]
  P=[]
  F=[]
  AMP=[]
  for mm=0:n_x
    for nn=0:n_y
      for pp=0:n_z
        fmnp=c/2*sqrt((mm/a)^2+(nn/b)^2+(pp/d)^2)
        #if abs(fmnp-f[u])<100*f[u]/Q[u]
        Amp=urc(Q[u],f[u],fmnp)
        if Amp>0.01
          M=vcat(M,mm)
          N=vcat(N,nn)
          P=vcat(P,pp)
          F=vcat(F,fmnp)
          AMP=vcat(AMP,Amp)
        end
      end
    end
  end
  nm=length(M)
  println("$(f[u]), $nm resonant freq.")
  f_coeffTMx_x=zeros(n_steps,nm)
  f_coeffTMx_y=zeros(n_steps,nm)
  f_coeffTMx_z=zeros(n_steps,nm)
  f_coeffTEx_y=zeros(Complex128,n_steps,nm)
  f_coeffTEx_z=zeros(Complex128,n_steps,nm)
  for pos=1:n_steps
    ETMx_x=zeros(Npt)
    ETMx_y=zeros(Npt)
    ETMx_z=zeros(Npt)
    ETEx_y=zeros(Complex128,Npt)
    ETEx_z=zeros(Complex128,Npt)
    xe=x_b[:,pos]
    ye=y_b[:,pos]
    ze=z_b[:,pos]
    #println(" Step $pos/$n_steps")
    for k=1:nm
      fmnp=F[k]
      Amp=AMP[k]
      kx=M[k]*pi/a
      ky=N[k]*pi/b
      kz=P[k]*pi/d
      k2mnp=kx^2+ky^2+kz^2
      ssc_ze=sum(sin(kx*xe).*sin(ky*ye).*cos(kz*ze))
      scs_ze=sum(sin(kx*xe).*cos(ky*ye).*sin(kz*ze))
      css_ze=sum(cos(kx*xe).*sin(ky*ye).*sin(kz*ze))
      ssc_xe=sum(sin(ky*ye).*sin(kz*ze).*cos(kx*xe))
      scs_xe=sum(sin(ky*ye).*cos(kz*ze).*sin(kx*xe))
      css_xe=sum(cos(ky*ye).*sin(kz*ze).*sin(kx*xe))
      ssc_ye=sum(sin(kz*ze).*sin(kx*xe).*cos(ky*ye))
      scs_ye=sum(sin(kz*ze).*cos(kx*xe).*sin(ky*ye))
      css_ye=sum(cos(kz*ze).*sin(kx*xe).*sin(ky*ye))

      if M[k]>0 && N[k]>0 && P[k]>0
        #TM modes z
        f_coeffTMx_z[pos,k] = -E0*Amp*kx*kz/(k2mnp-kz^2)*css_ze
        #f_coeffTMy_z[pos,k] = E0*Amp*ky*kz/(k2mnp-kz^2)*scs_ze
        #f_coeffTMz_z[pos,k] = E0*Amp*ssc_ze
        #TE modes z
        f_coeffTEx_z[pos,k] = -H0*Amp*1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kz^2)*css_ze
        #f_coeffTEy_z[pos,k] = H0*Amp*1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-kz^2)*scs_ze

        #TM modes x
        #f_coeffTMy_x[pos,k] = -E0*Amp*ky*kx/(k2mnp-kx^2)*css_xe
        #f_coeffTMz_x[pos,k] = E0*Amp*kz*kx/(k2mnp-kx^2)*scs_xe
        f_coeffTMx_x[pos,k] = E0*Amp*ssc_xe
        #TE modes x
        #f_coeffTEy_x[pos,k] = -H0*Amp*1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-kx^2)*css_xe
        #f_coeffTEz_x[pos,k] = H0*Amp*1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kx^2)*scs_xe

        ##TM modes y
        #f_coeffTMz_y[pos,k] += -H0*Amp*E0*Amp*kz*ky/(k2mnp-ky^2)*css_ye
        f_coeffTMx_y[pos,k] = H0*Amp*E0*Amp*kx*ky/(k2mnp-ky^2)*scs_ye
        #f_coeffTMy_y[pos,k] = H0*Amp*E0*Amp*ssc_ye
        ##TE modes y
        #f_coeffTEz_y[pos,k] = -H0*Amp*1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-ky^2)*css_ye
        f_coeffTEx_y[pos,k] = H0*Amp*1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-ky^2)*scs_ye
      end

      #if M[k]>0 && N[k]>0 && P[k]==0
        #f_coeffTMz_z[pos,k] = E0*Amp*ssc_ze #TM mode z
        #f_coeffTEz_y[pos,k] = -H0*Amp*1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-ky^2)*css_ye #TE mode y
        #f_coeffTEz_x[pos,k] = H0*Amp*1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kx^2)*scs_xe #TE mode x
      #end

      if M[k]==0 && N[k]>0 && P[k]>0
        f_coeffTEx_z[pos,k] += -H0*Amp*1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kz^2)*css_ze #TE mode z
        f_coeffTMx_x[pos,k] += E0*Amp*ssc_xe #TM mode x
        f_coeffTEx_y[pos,k] += H0*Amp*1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-ky^2)*scs_ye #TE mode y
      end

      #if M[k]>0 && N[k]==0 && P[k]>0
        #f_coeffTEy_z[pos,k] = H0*Amp*1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-kz^2)*scs_ze #TE mode z
        #f_coeffTMy_y[pos,k] = E0*Amp*ssc_ye #TM mode y
        #f_coeffTEy_x[pos,k] = -H0*Amp*1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-kx^2)*css_xe #TE mode x
      #end

      for i=1:Npt
        ssc_z=sin(kx*x[i])*sin(ky*y[i])*cos(kz*z[i])
        scs_z=sin(kx*x[i])*cos(ky*y[i])*sin(kz*z[i])
        css_z=cos(kx*x[i])*sin(ky*y[i])*sin(kz*z[i])
        ssc_x=sin(ky*y[i])*sin(kz*z[i])*cos(kx*x[i])
        scs_x=sin(ky*y[i])*cos(kz*z[i])*sin(kx*x[i])
        css_x=cos(ky*y[i])*sin(kz*z[i])*sin(kx*x[i])
        ssc_y=sin(kz*z[i])*sin(kx*x[i])*cos(ky*y[i])
        scs_y=sin(kz*z[i])*cos(kx*x[i])*sin(ky*y[i])
        css_y=cos(kz*z[i])*sin(kx*x[i])*sin(ky*y[i])

        if M[k]>0 && N[k]>0 && P[k]>0#TM modes z
          #TM modes z
          ETMx_z[i] += f_coeffTMx_z[pos,k]*css_z
          #ETMy_z[i] += f_coeffTMy_z[pos,k]*scs_z
          #ETMz_z[i] += f_coeffTMz_z[pos,k]*ssc_z

          #TE modes z
          ETEx_z[i] += f_coeffTEx_z[pos,k]*css_z
          #ETEy_z[i] += f_coeffTEy_z[pos,k]*scs_z

          #TM modes x
          #ETMy_x[i] += f_coeffTMy_x[pos,k]*css_x
          #ETMz_x[i] += f_coeffTMz_x[pos,k]*scs_x
          ETMx_x[i] += f_coeffTMx_x[pos,k]*ssc_x

          #TE modes x
          #ETEy_x[i] += f_coeffTEy_x[pos,k]*css_x
          #ETEz_x[i] += f_coeffTEz_x[pos,k]*scs_x

          ##TM modes y
          #ETMz_y[i] += f_coeffTMz_y[pos,k]*css_y
          ETMx_y[i] += f_coeffTMx_y[pos,k]*scs_y
          #ETMy_y[i] += f_coeffTMy_y[pos,k]*ssc_y

          ##TE modes y
          #ETEz_y[i] += f_coeffTEz_y[pos,k]*css_y
          ETEx_y[i] += f_coeffTEx_y[pos,k]*scs_y
        end

        #if M[k]>0 && N[k]>0 && P[k]==0
          #ETMz_z[i] += f_coeffTMz_z[pos,k]*ssc_z #TM mode z
          #ETEz_y[i]+= f_coeffTEz_y[pos,k]*css_y #TE mode y
          #ETEz_x[i]+= f_coeffTEz_x[pos,k]*scs_x #TE mode x
        #end

        if M[k]==0 && N[k]>0 && P[k]>0
          ETEx_z[i] += f_coeffTEx_z[pos,k]*css_z #TE mode z
          ETMx_x[i] += f_coeffTMx_x[pos,k]*ssc_x #TM mode x
          ETEx_y[i] += f_coeffTEx_y[pos,k]*scs_y #TE mode y
        end

        #if M[k]>0 && N[k]==0 && P[k]>0
          #ETEy_z[i] += f_coeffTEy_z[pos,k]*scs_z #TE mode z
          #ETMy_y[i] += f_coeffTMy_y[pos,k]*ssc_y #TM mode y
          #ETEy_x[i] += f_coeffTEy_x[pos,k]*css_x #TE mode x
        #end
      end
    end
    Ex[u,pos,:]=ETMx_x+ETMx_y+ETMx_z+ETEx_y+ETEx_z
    #Ey[u,pos,:]=ETMy_x+ETMy_y+ETMy_z+ETEy_x+ETEy_z
    #Ez[u,pos,:]=ETMz_x+ETMz_y+ETMz_z+ETEz_y+ETEz_x
  end

  f_coeffTMx=f_coeffTMx_x+f_coeffTMx_y+f_coeffTMx_z
  f_coeffTEx=f_coeffTEx_y+f_coeffTEx_z
  modesTE=zeros(n_z)
  modesTM=zeros(n_z)
  maxmax=maximum([maximum(abs(f_coeffTMx)),maximum(abs(f_coeffTEx))])
  for k=1:nm
    if maximum(abs(f_coeffTEx[:,k]))/maxmax>0.3
      modesTE[P[k]] += 1
    end
    if maximum(abs(f_coeffTMx[:,k]))/maxmax>0.3
      modesTM[P[k]] += 1
    end
  end

  ax=cor(squeeze(Ex[u,:,:],1)')
  Nx[u]=n_steps^2/sum(abs(ax).^2)
  Nmodesx[u]=sum((modesTE+modesTM).>0)
  println("Nx=$(Nx[u]), Nmodesx=$(Nmodesx[u])")

  figure(num=1,figsize=(14,10))
  ftext=round(f[u]*1000/1e6)/1000
  suptitle("\$ f=$ftext\$ MHz",fontsize=15)
  subplot(211)
  Nxtxt=round(Nx[u]*100)/100
  title("\$ \\vert E_x\\vert\$, \$ N_x=$Nxtxt \$ pos.")
  #imshow(abs(ax))
  #xlabel("\$z\$")
  #ylabel("\$z\$/cm")
  #colorbar(shrink=0.7)
  #clim(0,1)
  plot(phi,abs(squeeze(Ex[u,:,:],1)),"k",lw=0.3)
  xlim(0,2*pi)
  #xlabel("\$z\$/m")
  ylabel("V/m")
  grid()

  subplot(212)
  if Nmodesx[u]>1
    title("Main modes, \$N_m=$(Nmodesx[u])\$ uncorrelated modes")
  else
    title("Main modes, \$N_m=$(Nmodesx[u])\$ uncorrelated mode")
  end

  for i=1:nm
    if maximum(abs(f_coeffTEx[:,i]))/maxmax.>0.3
      Mtxt=M[i]
      Ntxt=N[i]
      Ptxt=P[i]
      label="$Mtxt \\ $Ntxt \\ $Ptxt"
      plot(phi,imag(f_coeffTEx[:,i])/maxmax,label="\$ TE^{yz}_{$label}\$")
    else
      plot(phi,imag(f_coeffTEx[:,i])/maxmax,"k",lw=0.3)
    end

    if maximum(abs(f_coeffTMx[:,i]))/maxmax.>0.3
      Mtxt=M[i]
      Ntxt=N[i]
      Ptxt=P[i]
      label="$Mtxt \\ $Ntxt \\ $Ptxt"
      plot(phi,f_coeffTMx[:,i]/maxmax,"--",label="\$ TM^{xyz}_{$label}\$")
    else
      plot(phi,f_coeffTMx[:,i]/maxmax,"--k",lw=0.3)
    end
  end
  xlim(0,2*pi)
  ylim(-1.1,1.1)
  grid()
  xlabel("\$ \\phi \$/rad")
  ylabel("A.U.")
  legend(fontsize="small")
  savefig("N_1stirrer$u.png")
  clf()
  M=[]
  N=[]
  P=[]
end

plot(f,Nx)
plot(f,Nmodesx)
grid()

using NPZ
npzwrite("Exstirrer.npz", ["Ex" => Ex, "f" => f, "Q" => Q, "x" => x,"y" => y,"z" => z,"Nx"=>Nx,"Nmodesx"=>Nmodesx])


text:
function maximumsimultlicenses(io::IO)
    out, maxout, maxtimes = 0, -1, String[]
    for job in readlines(io)
        out += ifelse(occursin("OUT", job), 1, -1)
        if out > maxout
            maxout = out
            empty!(maxtimes)
        end
        if out == maxout
            push!(maxtimes, split(job)[4])
        end
    end
    return maxout, maxtimes
end

let (maxout, maxtimes) = open(maximumsimultlicenses, "data/mlijobs.txt")
    println("Maximum simultaneous license use is $maxout at the following times: \n - ", join(maxtimes, "\n - "))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Adapted from Example 3.99 of [BPT12].
#
# [BPT12] Blekherman, G.; Parrilo, P. & Thomas, R.
# Semidefinite Optimization and Convex Algebraic Geometry
# Society for Industrial and Applied Mathematics, 2012

using Test
using SumOfSquares
using DynamicPolynomials

function BPT12e399_test(optimizer, config::MOIT.TestConfig)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)

    @polyvar x y
    cref = @constraint(model, 10 - (x^2 + α*y) in SOSCone(),
                       domain = @set x^2 + y^2 == 1)

    @objective(model, Max, α)
    JuMP.optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 6.0 atol=atol rtol=rtol

    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    @test JuMP.value(α) ≈ 6 atol=1e-6

    @test_throws SumOfSquares.ValueNotSupported value(cref)
    p = gram_matrix(cref)
    @test getmat(p) ≈ [1 -3; -3 9] atol=atol rtol=rtol
    @test p.x == [y, 1]

    @test dual_status(model) == MOI.FEASIBLE_POINT
    μ = dual(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 3
    @test moment_value(moments(μ)[1]) ≈ -8/3 atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == x^2
    @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
    @test monomial(moments(μ)[2]) == y
    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
    @test monomial(moments(μ)[3]) == 1

    μ = moments(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 3
    @test moment_value(moments(μ)[1]) ≈ 3 atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == y^2
    @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
    @test monomial(moments(μ)[2]) == y
    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
    @test monomial(moments(μ)[3]) == 1

    @objective(model, Min, α)
    JuMP.optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ -6.0 atol=atol rtol=rtol

    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    @test JuMP.value(α) ≈ -6 atol=1e-6

    @test_throws SumOfSquares.ValueNotSupported value(cref)
    p = gram_matrix(cref)
    @test getmat(p) ≈ [1 3; 3 9] atol=atol rtol=rtol
    @test p.x == [y, 1]

    @test dual_status(model) == MOI.FEASIBLE_POINT
    μ = dual(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 3
    @test moment_value(moments(μ)[1]) ≈ -8/3 atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == x^2
    @test moment_value(moments(μ)[2]) ≈ -1 atol=atol rtol=rtol
    @test monomial(moments(μ)[2]) == y
    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
    @test monomial(moments(μ)[3]) == 1

    μ = moments(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 3
    @test moment_value(moments(μ)[1]) ≈ 3 atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == y^2
    @test moment_value(moments(μ)[2]) ≈ -1 atol=atol rtol=rtol
    @test monomial(moments(μ)[2]) == y
    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
    @test monomial(moments(μ)[3]) == 1
end

sd_tests["BPT12e399"] = BPT12e399_test


text:
#################################################################################
#
# sed_1d_v5.jl - a Julia script for modeling 1-D sediment transport
#
# concept is to exchange sediments in water column with sediment bed layers
# based on differences in river sediment-carrying capacity
#
#################################################################################

### constants

const rho_w = 1000.		# density of suspending medium (~water)
const rho_s = 2650.     # density of sediment particles
const rho_b = 1600.     # bulk sediment density
const g = 9.807         # gravitational acceleration
const nu = 0.9e-6       # kinematic viscosity of water (m^2/sec)
const s = rho_s/rho_w  	# sediment particle specific gravity

### sediment physical characteristics

type Sediment
	name::AbstractString
    d50::Float64
    d90::Float64
    floc::Float64 					# multiplier for increasing d50 to account for floccing
    upgradient::Float64 			# fractional abundance of this sediment class included in the upgradient boundary condition
	dnd::Float64 					# dimensionless particle size
    vs::Float64 					# settling velocity
end

function AssignSeds()
	# read sediment characteristics file and assign individual sediments to respective type
	sediment = Sediment[]
	data = readdlm("sediments.txt", '\t', header=true)
	for i = 1:size(data[1], 1)
		name = data[1][i, 1]
		d50 = Float64(data[1][i, 2])
		d90 = Float64(data[1][i, 3])
		floc = Float64(data[1][i, 4])
		upgradient = Float64(data[1][i, 5])
		dnd = Dnd(d50, 1.0)
		vs = SettlingVel(d50, floc)
		push!(sediment, Sediment(name, d50, d90, floc, upgradient, dnd, vs))
	end
	return sediment
end

function Dnd(d50::Float64, floc::Float64)
	# non-dimensional particle diameter
    return ((s-1.)*g/nu^2)^(1./3.) * (d50 * floc)
end

function SettlingVel(d50::Float64, floc::Float64)
	# particle settling velocity (Cheng, 1997 approximation)
    return nu/(d50*floc) * ((25. + 1.2*Dnd(d50,floc)^2)^0.5 - 5.)^1.5
end

### dredging history

type Dredge
	t_ex::Array{Float64, 1}			# time of dredging event
	dz::Array{Float64, 2} 			# dredging depth; rows = time events, columns = cells
end

function AssignDredge()
	# read dredge event file and assign all to Dredge type
	data = readdlm("dredge.txt", '\t', header=true)
	t_ex = data[1][:,1]
	dz = data[1][:,2:end]
	dredge = Dredge(t_ex, dz)
	return dredge
end

### model parameters and supporting functions

type Model
    num_cells::Int64 			# number of river cells included in the model
	w_up::Float64 				# wdith of river section upgradient of first model cell
    hw::Float64 				# height (thickness) of water column
    elev_w::Float64 			# elevation of top of water column above datum
    dt::Float64 				# time step size (fixed)
    t_max::Float64 				# simulation end time
    t_profile::Float64 			# time step size for writing profile output files
    t_series::Float64 			# time step size for monitoring output
    num_slices::Int64 			# (initial) number of slices for sediment bed
    mixing_depth::Float64 		# scour/mixing depth (fixed)
    mode_Q::Int64 				# 0 = fixed discharge rate; 1 = variable (probabilistic) discharge rate
    log_min_Q::Float64 			# minimum discharge rate (log scale)
    log_max_Q::Float64 			# maximum discharge rate (log scale)
    log_avg_Q::Float64 			# average discharge rate (log scale)
    log_stdev_Q::Float64 		# standard deviation of the log discharge rate
    dt_Q::Float64 				# time interval after which to select a new discharge rate if mode_Q == 1
end

function ModelParams()
	# read model parameters file and assign values to Model type
	data = readdlm("model_params.txt", '\t', header=false)
	num_cells = Int64(data[1, 2])
	w_up = Float64(data[2, 2])
	hw = Float64(data[3, 2])
	elev_w = Float64(data[4, 2])
	dt = Float64(data[5, 2])
	t_max = Float64(data[6, 2])
	t_profile = Float64(data[7, 2])
	t_timeseries = Float64(data[8, 2])
	num_slices = Int64(data[9, 2])
	mixing_depth = Float64(data[10, 2])
	mode_Q = Int64(data[11, 2])
	log_min_Q = Float64(data[12, 2])
	log_max_Q = Float64(data[13, 2])
	log_avg_Q = Float64(data[14, 2])
	log_stdev_Q = Float64(data[15, 2])
	dt_Q = Float64(data[16, 2])
	model = Model(num_cells, w_up, hw, elev_w, dt, t_max, t_profile, t_timeseries, num_slices, mixing_depth, mode_Q, log_min_Q, log_max_Q, log_avg_Q, log_stdev_Q, dt_Q)
	return model
end

function Velocity(Q::Float64, h::Float64, w::Float64)
    # mean flow velocity
    return Q/(w * h)
end

function Me(q::Float64, ucr::Float64, sed::Sediment)
    # mobility parameter for simplified load models (for all class sizes)
    if q > ucr
		me = (q-ucr)/((s-1)*g*sed.d50)^0.5
    else
		me = 0.
	end
    return me
end

function Ucr(h::Float64, sed::Sediment)
    # critical velocity for currents (for all class sizes)
    if sed.d50 < 0.0005
		ucr = 0.19*sed.d50^0.1 * log10(12*h/(3.*sed.d90))
    else
		ucr = 8.5*sed.d50^0.6 * log10(12*h/(3.*sed.d90))
	end
    return ucr
end

function TransBL(q::Float64, h::Float64, sed::Sediment)
    # simplified bed load transport formula for steady-state current flow (for all class sizes)
    ucr = Ucr(h, sed)
    me = Me(q, ucr, sed)
    return 0.015 * rho_s * q * h * (sed.d50/h)^1.2 * me^1.5      # dimensions = mass length^-1 time^-1
end

function TransSL(q::Float64, h::Float64, sed::Sediment)
	# simplified suspended load transport formula for steady-state current flow
    ucr = Ucr(h, sed)
    M_e = Me(q, ucr, sed)
    return 0.012 * rho_s * q * sed.d50 * M_e^2.4 * sed.dnd^(-0.6)    # dimensions = mass length^-1 time^-1
end

### stratigraphy: data and methods

type Strat
    dz::Float64 									# model individual slice thickness (fixed at the start of the simulation)
    active_dz::Float64 								# thickness of the top (active) slice, the only slice which can have a thickness != dz
    lith_frac::Array{Array{Float64,1},1}			# array of volume proportions of different sediment classes represented in each slice
end

function AssignStrat(model::Model)

	# read stratigraphy file (per model cell), which specifies sediment layer characteristics
	data = readdlm("stratigraphy.txt", '\t', header=true)
	strat = Strat[]
	num_seds = size(data[1], 2)-3 					# note: sediment ordering is assumed to be preserved in input file; no check is made ...
	start_line = zeros(Int64, model.num_cells) 		# starting point for cell[i] in the table
	num_layers = zeros(Int64, model.num_cells) 		# number of layers corresponding to cell[i] in the table
	cell_indices = data[1][:,1]
	for i = 1:model.num_cells
		start_line[i] = minimum(findin(cell_indices,i))
		num_layers[i] = length(findin(cell_indices,i))
	end

	for icell = 1:model.num_cells 					# for each cell
		bottom = Float64[]
		top = Float64[]
		lith_frac_matrix = Array[] 					# this is just a container for sets of Float64 arrays
		for i = start_line[icell]:(start_line[icell] + num_layers[icell]) - 1		# for each defined layer
			push!(bottom, Float64(data[1][i, 2]))
			push!(top, Float64(data[1][i, 3]))
			lith_frac_layer = Float64[]
			for j = 1:num_seds
				push!(lith_frac_layer, Float64(data[1][i, 3+j]))
			end
			push!(lith_frac_matrix, copy(lith_frac_layer))
		end

		# distribute layer characteristics across appropriate model slices
		dz = (maximum(top) - minimum(bottom))/model.num_slices     			# model slice thickness will always remain fixed, even for new added layers
		active_dz = dz                            							# initial active slice thickness
		lith_frac = Array[] 												# just a container for sets of Float64 arrays
		for i = 1:model.num_slices
			midpoint = (i+0.5) * dz
			for j = 1:num_layers[icell]
				if midpoint >= bottom[j] && midpoint < top[j]
					push!(lith_frac, copy(lith_frac_matrix[j]))  			# this 2D array is all that ever needs to be tracked subsequently, since dz is already specified
					break
				end
			end
		end
		push!(strat, Strat(dz, active_dz, lith_frac))

	end

	return strat

end

function Dig(strat::Strat, water_depth::Float64, dredge_depth::Float64)
    # remove sediment to match dredging event authorized depth
	dig = max(dredge_depth - water_depth, 0.0)
    num_dig_layers = Int64(floor(dig/strat.dz))
	for i = 1:num_dig_layers
		deleteat!(strat.lith_frac, 1)
	end
    strat.active_dz = strat.dz
	return dig
end

function Blend(strat::Strat, d::Float64)
	# determine size class composition of mixing zone at top of sediment bed
	blend_frac = zeros(strat.lith_frac[1]) 			# composite sediment class composition of mixed zone (initialize array)
	blend_w = Float64[] 							# weighting fractions associated with active slices
	d_1 = max(d - strat.active_dz,0.)        		# extent blending zone beneath active layer
	d_2 = d_1 % strat.dz 							# division remainder
	num_mix_layer = Int64(floor(d_1/strat.dz)) 		# division quotient
	# weight the influences of blending layers on average sediment composition
	push!(blend_w,min(strat.active_dz/d, 1.))       # active layer
	for i = 1:num_mix_layer
		push!(blend_w, strat.dz/d) 					# complete layers underneath active layer
	end
	if d_2 > 0.
		push!(blend_w, d_2/d)                  		# partial layer on bottom
	end
	for j = 1:length(blend_frac), i = 1:length(blend_w)
		blend_frac[j] += blend_w[i] * strat.lith_frac[i][j]
	end
	return blend_frac
end

function Deposit(strat::Strat, d::Float64, blended_lith::Array{Float64,1}, dh::Array{Float64,1})

    # mix in deposited material (positive) or remove (negative) from blended zone
	for j = 1:length(dh)
		dh[j] = max(dh[j],-d*blended_lith[j])          # net change in sediment thickness
	end
    delta_d = sum(dh)

    # (1) determine updated average composition of sediment mixture down to the blending depth
	for j = 1:length(dh)
		blended_lith[j] = ((d-delta_d)*blended_lith[j] + dh[j])/d
	end

	# (2) determine any layering changes
	if delta_d + strat.active_dz > strat.dz
		# need to add new layer(s)

		global num_added_layers = Int64(floor((delta_d + strat.active_dz)/strat.dz))		# division quotient

		strat.active_dz = (delta_d + strat.active_dz) % strat.dz 					# division remainder
		for i = 1:num_added_layers
			insert!(strat.lith_frac, 2, copy(blended_lith)) 						# insert just beneath active layer
		end
	elseif (delta_d < 0.) && (delta_d + strat.active_dz <= 0.)
		# need to remove layer(s)
		del_active_dz = abs(delta_d + strat.active_dz) % strat.dz 						# division remainder
		num_del_layers = Int64(floor(abs(delta_d + strat.active_dz)/strat.dz)) + 1		# division quotient (plus one layer)
		strat.active_dz = strat.dz - del_active_dz
		for i = 1:num_del_layers
			deleteat!(strat.lith_frac, 1)
		end
	else
		# active layer thickness changes, but no change in layer numbering
		strat.active_dz += delta_d
	end

	# (3) apply compositional changes to top layers from surface down to depth not compensated by layering changes
    strat.lith_frac[1] = copy(blended_lith) 										# active layer
    comp_depth = max(d - strat.active_dz, 0.)         								# depth into sediment pile, beyond active layer
	bottom_blend_d = comp_depth % strat.dz 											# division remainder
	num_comp_layers = Int64(floor(comp_depth/strat.dz))								# division quotient
    for i = 1:num_comp_layers
		strat.lith_frac[i+1] = copy(blended_lith) 									# layers beneath active layer
	end
	# bottom boundary of mixing zone
	for j = 1:length(dh)
		strat.lith_frac[num_comp_layers+2][j] = (bottom_blend_d/strat.dz)*blended_lith[j] + (1.-bottom_blend_d/strat.dz)*strat.lith_frac[num_comp_layers+2][j]
	end

end

### external sources

type Source
	t_src::Array{Array{Float64,1},1} 					# list of source time arrays (one array per active-source cell-sediment combo)
	J::Array{Array{Float64,1},1} 						# matching list of flux arrays
	index::Array{Tuple{Int64, Int64}, 1} 				# list of tuples (cell, material index)
    load::Array{Float64, 1} 							# cumulative loading, per sediment size class (i.e., 1-D array)
	sed_indices::Array{Int64,1} 						# set of sediment list indices employed as sources
	sed_names::Array{AbstractString,1} 					# list of sediment names accompanying sed_indices
end

function AssignSources(sediment::Array{Sediment,1})

	# assign t_src and J arrays to source arrays; create corresponding list of tuples (cell, material index) to track source arrays

	t_src  = Array[] 									# initial multi-source time and flux array container
	J = Array[]
	index_list = Tuple{Int64, Int64}[] 					# initialize list of tuples to hold cell and sediment indices for source(s)
	load = zeros(Float64,length(sediment)) 				# initialize cumulative loading array (size = number of sediment particle classes)
	sed_name = AbstractString[] 						# names of all sediment classes in sediment
	sed_names = AbstractString[] 						# subset of names of those classes with external sources
	sed_indices = Int64[] 								# sediment class indices associated with external sources
	sed_names = AbstractString[] 						# subset of names of those classes with external sources
	time_array = Float64[] 								# temporary holders for housing time and flux arrays
	flux_array = Float64[]

	for i = 1:length(sediment)
		push!(sed_name, sediment[i].name)
	end
	data = readdlm("sources.txt", '\t', header=true)

    for i = 1:1:size(data[1], 1)

		cell_index = data[1][i, 4]
		sed_index = findin(sed_name,[data[1][i, 1]])[1]
		push!(sed_indices, sed_index)

		if i == 1
			# this is a new cell and/or sediment type (in source definition)
			push!(index_list, (cell_index, sed_index))
			push!(time_array, data[1][i, 2])
			push!(flux_array, data[1][i, 3])

		elseif (i > 1) && ((cell_index, sed_index) != index_list[end])
			# this is a new cell and/or sediment type (in source definition); data exists for a prior source term
			push!(t_src, time_array)
			push!(J, flux_array)
			push!(index_list, (cell_index, sed_index))
			time_array = Float64[]
			flux_array = Float64[]
			push!(time_array, data[1][i, 2])
			push!(flux_array, data[1][i, 3])

		else
			# this is more of the same source term; add to time and mass flux arrays
			push!(time_array, data[1][i, 2])
			push!(flux_array, data[1][i, 3])

		end

	end

	# capture final read
	push!(t_src, time_array) #####
	push!(J, flux_array)

	# match names on sed_name list to sed_indices
	sed_indices = sort(collect(Set(sed_indices)))
	for i in sed_indices
		push!(sed_names, sed_name[i])
	end

	# populate source type
	source = Source(t_src, J, index_list, load, sed_indices, sed_names)

	#println(source.t_src[1])

	return source

end

function SourceTerm(icell::Int64, ised::Int64, t::Float64, source::Source)
	# return a matching source flux for ised in icell at t
	loc = (icell, ised)
	find_flag = false
	source_index = 0
	for (i, value) in enumerate(source.index)
		if loc == value
			find_flag = true 							# note that i is the index to use in the t and J lists
			source_index = i
			break
		end
	end
	if find_flag == true
		# compute flux from arrays
		J_index = 1 + round(Int, (t - source.t_src[source_index][1]) / (source.t_src[source_index][end] - source.t_src[source_index][1]) * length(source.t_src[source_index]))
		if (J_index < 1) || (J_index > length(source.t_src[source_index]))
			flux = 0.
		else
			flux = source.J[source_index][J_index]
		end
	else
		flux = 0.
	end
	return flux
end

### cell geometry

type Cell
    w::Float64 					# river cell width
    area::Float64 				# areal footprint of the river cell
end

function CellGeom(num_cells::Int64)
	# read cell geometry specs file
	cell = Cell[]
	data = readdlm("cells.txt", '\t', header=true)
	for i = 1:num_cells
		w = Float64(data[1][i, 2])
		L = Float64(data[1][i, 3])
		area = w * L
		push!(cell, Cell(w, area))
	end
	return cell
end

### utility functions

function SettlingDep(J0_SL::Array{Float64}, J0_BL::Array{Float64}, J_SL::Array{Float64,2}, J_BL::Array{Float64,2}, sediment::Array{Sediment,1}, source::Source, model::Model, Q::Float64, t::Float64)
	# compute net deposition per time step using the settling velocity method (for comparison to model algorithm)
	dS_2 = zeros(Float64, model.num_cells) 							# net deposition per time step, settling velocity method
	num_seds = length(sediment)
	dep = zeros(Float64, num_seds) 							# sediment height change, per sediment type
	for icell = 1:model.num_cells
		if icell == 1
			for i = 1:num_seds
				C = ((0.5*J0_SL[i] + 0.5*J_SL[icell, i]) + (0.5*J0_BL[i] + 0.5*J_BL[icell, i])  +  SourceTerm(icell, i, t, source)) / Q
				flux = C * sediment[i].vs
				dep[i] = model.dt * flux / rho_b
			end
		else
			for i = 1:num_seds
				C = ((0.5*J_SL[icell-1, i] + 0.5*J_SL[icell, i]) + (0.5*J_BL[icell-1, i] + 0.5*J_BL[icell, i])  +  SourceTerm(icell, i, t, source)) / Q
				flux = C * sediment[i].vs
				dep[i] = model.dt * flux / rho_b
			end
		end
		dS_2[icell] = sum(dep)
	end
	return dS_2
end

function RandNorm(avg::Float64, stdev::Float64, lower_bound::Float64, upper_bound::Float64)
	# clipped,normally-distributed random number generation
	r = 0.
	accept = 0
	while accept == 0
		r = avg + stdev*randn(1)[1]
		if (r >= lower_bound) && (r <= upper_bound)
			accept = 1
		end
	end
	return r
end

function WaterDepth(model::Model, strat::Array{Strat, 1})
	# return updated array of water depths (i.e., water-column thicknesses) across all model cells
	w_depth = zeros(Float64, model.num_cells)
	for icell = 1:model.num_cells
		w_depth[icell] = model.elev_w - strat[icell].dz*(length(strat[icell].lith_frac)-1) + strat[icell].active_dz
	end
	return w_depth
end

### output functions

function WriteProfile(fname::AbstractString, strat::Strat, sediment::Array{Sediment, 1})
	# write sediment class volume proportions, by layer, to output file
	csvfile = open(fname, "w")
	line_out = "z"
	for i = 1:length(sediment)
		line_out = line_out * "," * sediment[i].name
	end
	println(csvfile,line_out)
	num_layers = length(strat.lith_frac)
	for i = 1:num_layers
		z = (num_layers - (i+0.5))*strat.dz
		line_out = string(z)
		for j = 1:length(sediment)
			line_out = line_out * "," * string(strat.lith_frac[i][j])
		end
		println(csvfile,line_out)
	end
	close(csvfile)
end

function WriteThickSeries(t::Float64, Q::Float64, strat::Array{Strat, 1}, init_flag::Bool)
    # write select variables to time series file
	fname = "sed_thick.csv"
	num_cells = length(strat)
	if init_flag == true
		# create time series output file and header row
		csvfile = open(fname,"w")
		line_out = "t" * "," * "Q"
		for i = 1:num_cells
			line_out = line_out * "," * "cell_" * string(i)
		end
		println(csvfile,line_out)
	else
		csvfile = open(fname,"a")
	end
	# update ...
	line_out = string(t) * "," * string(Q)
	for i = 1:num_cells
		thick = (length(strat[i].lith_frac)-1)*strat[i].dz + strat[i].active_dz
		line_out = line_out * "," * string(thick)
	end
	println(csvfile,line_out)
	close(csvfile)
end

function WriteDepSeries(t::Float64, Q::Float64, dS_1::Array{Float64}, dS_2::Array{Float64}, init_flag::Bool)
	# write deposition thicknesses, per time step, to output file
	fname = "dep_compare.csv"
	num_cells = length(dS_1)
	if init_flag == true
		# create header row
		csvfile = open(fname,"w")
		line_out = "t" * "," * "Q"
		for j = 1:2
			for i = 1:num_cells
				line_out = line_out * "," * "dS_" * string(j) * "_cell_" * string(i)
			end
		end
		println(csvfile,line_out)
	else
		csvfile = open(fname,"a")
	end
	# update ...
	line_out = string(t) * "," * string(Q)
	# depositional calc method #1 (employed): differences in sediment carrying capacity
	for i = 1:num_cells
		line_out = line_out * "," * string(dS_1[i])
	end
	# depositional calc method #2: settling velocity X concentration
	for i = 1:num_cells
		line_out = line_out * "," * string(dS_2[i])
	end
	println(csvfile,line_out)
	close(csvfile)
end

function WriteSourceSeries(t::Float64, strat::Array{Strat, 1}, source::Source, cell::Array{Cell,1}, init_flag::Bool)
    # write select variables to time series file
	fname = "external_residual.csv"
	num_cells = length(strat)
	if init_flag == true
		# create time series output file and header row
		csvfile = open(fname,"w")
		line_out = "t"
		for sed in source.sed_names
			line_out = line_out * "," * sed * "_input"
			for icell = 1:num_cells
				line_out = line_out * "," * sed * "_cell_" * string(icell)
			end
		end
		println(csvfile,line_out)
	else
		csvfile = open(fname,"a")
	end

	# update ...
	line_out = string(t)
	for ised in source.sed_indices
		line_out = line_out * "," * string(source.load[ised])
		for icell = 1:num_cells
			sum_mass = strat[icell].lith_frac[1][ised] * strat[icell].active_dz
			for islice = 2:length(strat[icell].lith_frac)
				sum_mass += strat[icell].lith_frac[islice][ised] * strat[icell].dz
			end
			sum_mass *= cell[icell].area * rho_b
			line_out = line_out * "," * string(sum_mass)
		end
	end

	println(csvfile,line_out)
	close(csvfile)
end

function WriteSpoils(t::Float64, tot_spoils::Float64, dredge_count::Int64)
	# write dredging spoils uncorrected volume estimates, per dredge evenet, to output file
	fname = "spoils.csv"
	if dredge_count == 1
		# create header row
		csvfile = open(fname,"w")
		line_out = "t" * "," * "spoils_vol"
		println(csvfile, line_out)
	else
		csvfile = open(fname,"a")
	end
	# update ...
	line_out = string(t) * "," * "," * string(tot_spoils)
	println(csvfile,line_out)
	close(csvfile)
end

function WriteSedFluxes(t::Float64, J_influx::Array{Float64, 1}, J_outflux::Array{Float64, 1}, sediment::Array{Sediment, 1}, init_flag::Bool)
	# update cumulative sediment influxes and outflows, by sediment class
	fname = "sed_fluxes.csv"

	if init_flag == true
		# create time series output file and header row
		csvfile = open(fname,"w")
		line_out = "t"
		for i = 1:length(sediment)
			line_out = line_out * "," * sediment[i].name * "_in"
		end
		for i = 1:length(sediment)
			line_out = line_out * "," * sediment[i].name* "_out"
		end
		println(csvfile,line_out)
	else
		csvfile = open(fname,"a")
	end

	# update ...
	line_out = string(t)
	for i = 1:length(sediment)
		line_out = line_out * "," * string(J_influx[i])
	end
	for i = 1:length(sediment)
		line_out = line_out * "," * string(J_outflux[i])
	end
	println(csvfile,line_out)
	close(csvfile)

end

### main script

function SedTran()

	println("Initializing sediment-1D model ...")

	sediment = AssignSeds()
	println("... read sediment characteristics")

	model = ModelParams()
	println("... read model parameters")

	cell = CellGeom(model.num_cells)
	println("... read cell geometry")

	strat = AssignStrat(model)
	println("... read initial stratigraphy")

	source = AssignSources(sediment)
	println("... read external source terms.")

	dredge = AssignDredge()
	dredge_count = 0
	println("... read dredge history")

	# initialize flux containers
	t = 0.
	num_seds = length(sediment)
	J0_SL = zeros(Float64, num_seds) 								# incoming (upstream) suspended load fluxes
	J0_BL = zeros(Float64, num_seds)								# incoming (upstream) bed load fluxes
	dh = zeros(Float64,num_seds) 									# change in sediment thickness, per size class
	J_SL = zeros(Float64,(model.num_cells, num_seds)) 				# outgoing (current-cell) suspended load fluxes; 1st = cell index, 2nd = sediment size class
	J_BL = zeros(Float64,(model.num_cells, num_seds)) 				# outgoing (current-cell) bed load fluxes; 1st = cell index, 2nd = sediment size class
	dS_1 = zeros(Float64, model.num_cells) 							# net deposition per time step, sediment-carrying capacity method
	J_influx = zeros(Float64, num_seds)								# cumulative sediment fluxes, by size class, across upstream boundary
	J_outflux = zeros(Float64, num_seds)							# cumulative sediment fluxes, by size class, out of last cell

	# initialize discharge
	if model.mode_Q == 1
		Q = 10.^RandNorm(model.log_avg_Q, model.log_stdev_Q, model.log_min_Q, model.log_max_Q)
	else
		Q = 10.^model.log_avg_Q
	end

	# time-step loop ...

	while t < model.t_max

		t += model.dt

		# remove sediment thickness for a simplified dredging event
		dredge_match = findin(dredge.t_ex,[t]) 						# check if current time corresponds to a dredging event
		if length(dredge_match) == 1
			dredge_count += 1
			spoils = zeros(Float64, model.num_cells) 						# depth excavated (per cell)
			w_depth = WaterDepth(model, strat)   							# compute/update water depth array
			println("Modeling dredge event # ", dredge_count, ", time = ", t)
			for i = 1:model.num_cells
				dig = Dig(strat[i], w_depth[i], dredge.dz[dredge_count, i])
				spoils[i] = dig * cell[i].area
			end
			WriteSpoils(t, sum(spoils), dredge_count) 						# append to dredging spoils volume file
		end

		# boundary condition for first cell
		q = Velocity(Q, model.hw, model.w_up)
		for i = 1:num_seds
			J0_SL[i] = sediment[i].upgradient * TransSL(q, model.hw, sediment[i]) * model.w_up
			J0_BL[i] = sediment[i].upgradient * TransBL(q, model.hw, sediment[i]) * model.w_up
		end

		# sediment fluxes per model cell
		blended_lith = Array[]
		w_depth = WaterDepth(model, strat)   							# compute/update water depth array
		for icell = 1:model.num_cells

			# average current sediment composition across blending depth for icell
			push!(blended_lith, Blend(strat[icell], model.mixing_depth))

			# outgoing water-column fluxes
			if icell == 1
				q = 0.5*Velocity(Q, model.hw, model.w_up) + 0.5*Velocity(Q, w_depth[icell], cell[icell].w)
			else
				q = 0.5*Velocity(Q, w_depth[icell-1], cell[icell-1].w) + 0.5*Velocity(Q, w_depth[icell], cell[icell].w)
			end

			for i = 1:num_seds
				J_SL[icell,i] = blended_lith[icell][i] * TransSL(q, w_depth[icell], sediment[i]) * cell[icell].w
				J_BL[icell,i] = blended_lith[icell][i] * TransBL(q, w_depth[icell], sediment[i]) * cell[icell].w
			end

		end

		# track cumulative sediment fluxes across boundaries
		for i = 1:num_seds
			J_influx[i] += model.dt * (J0_SL[i] + J0_BL[i])
			J_outflux[i] += model.dt * (J_SL[end, i] + J_BL[end, i])
		end

		# calculate implied net deposition from differences in sediment carrying capacity + source term(s)
		for icell = 1:model.num_cells
			if icell == 1 				# first cell in river sees incoming sediment fluxed from upgradient boundary condition
				for i = 1:num_seds
					dh[i] = model.dt * ((J0_SL[i]-J_SL[icell, i]) + (J0_BL[i]-J_BL[icell, i])  +  SourceTerm(icell, i, t, source)) / (rho_b * cell[icell].area)
				end
			else 						# cell is downstream from boundary condition
				for i = 1:num_seds
					dh[i] = model.dt * ((J_SL[icell-1, i]-J_SL[icell, i]) + (J_BL[icell-1, i]-J_BL[icell, i])  +  SourceTerm(icell, i, t, source)) / (rho_b * cell[icell].area)
				end
			end
			dS_1[icell] = sum(dh)
			if dS_1[icell] > model.mixing_depth
				println("Time step is too large; dS_1 = ", dS_1[icell])
				assert(1==0) 		# throw an exception
			end

			# modify active layer (by mixing), or change active layer by strat layer addition or subtraction
			Deposit(strat[icell], model.mixing_depth, blended_lith[icell], dh)
		end

		# summarize source loading
		for loc in source.index
			source.load[loc[2]] += SourceTerm(loc[1], loc[2], t, source) * model.dt
		end

		# write profiles to output files, as warranted
		if (t % model.t_profile == 0) || (dredge_match == 1)
			println("Writing output at t = " * string(t) * ", Q = " * string(Q))
			for icell = 1:model.num_cells
				fname = string(icell) * "_" * string(t) * "_out.csv"
				WriteProfile(fname, strat[icell], sediment)
			end
		end

		# write/append to time series to output files, as warranted
		if t % model.t_series == 0
			init_flag = (t == model.t_series)
			dS_2 = SettlingDep(J0_SL, J0_BL, J_SL, J_BL, sediment, source, model, Q, t)
			WriteDepSeries(t, Q, dS_1, dS_2, init_flag)
			WriteThickSeries(t, Q, strat, init_flag)
			WriteSourceSeries(t, strat, source, cell, init_flag)
			WriteSedFluxes(t, J_influx, J_outflux, sediment, init_flag)
		end

		# select a new value for Q, if warranted by t (and mode_Q selection)
		if (model.mode_Q == 1) && (t % model.dt_Q == 0)
			Q = 10.^RandNorm(model.log_avg_Q, model.log_stdev_Q, model.log_min_Q, model.log_max_Q)
		end

	end

	println("Done.")

end

### run script

SedTran()


text:
using KernelAbstractions
using Test

import KernelAbstractions.NDIteration: NDRange, StaticSize, DynamicSize

@kernel function index(A)
    I = @index(Global, NTuple)
    @show A[I...]
end

@kernel function literal_pow(A)
    A[1] = 2^11
end

@kernel function square(A, B)
    A[1] = B[1]^2
end

@kernel function pow(A, B)
    A[1] = A[1]^B[1]
end

@kernel function checked(A, a, b)
    A[1] = Base.Checked.checked_add(a, b)
end

function check_for_overdub(stmt)
    if stmt isa Expr
        if stmt.head == :invoke
            mi = first(stmt.args)::Core.MethodInstance
            if mi.def.name === :overdub
                @show stmt
                return true
            end
        end
    end
    return false
end

function compiler_testsuite(backend, ArrayT)
    kernel = index(CPU(), DynamicSize(), DynamicSize())
    iterspace = NDRange{1, StaticSize{(128,)}, StaticSize{(8,)}}();
    ctx = KernelAbstractions.mkcontext(kernel, 1, nothing, iterspace, Val(KernelAbstractions.NoDynamicCheck()))
    @test KernelAbstractions.__index_Global_NTuple(ctx, CartesianIndex(1)) == (1,)

    A = ArrayT{Int}(undef, 1)
    let (CI, rt) = @ka_code_typed literal_pow(backend())(A, ndrange=1)
        # test that there is no invoke of overdub
        @test !any(check_for_overdub, CI.code)
    end

    A = ArrayT{Float64}(undef, 1)
    let (CI, rt) = @ka_code_typed square(backend())(A, A, ndrange=1)
        # test that there is no invoke of overdub
        @test !any(check_for_overdub, CI.code)
    end

    A = ArrayT{Float64}(undef, 1)
    B = ArrayT{Float64}(undef, 1)
    let (CI, rt) = @ka_code_typed pow(backend())(A, B, ndrange=1)
        # test that there is no invoke of overdub
        @test !any(check_for_overdub, CI.code)
    end

    A = ArrayT{Float64}(undef, 1)
    B = ArrayT{Int32}(undef, 1)
    let (CI, rt) = @ka_code_typed pow(backend())(A, B, ndrange=1)
        # test that there is no invoke of overdub
        @test !any(check_for_overdub, CI.code)
    end

    A = ArrayT{Int}(undef, 1)
    let (CI, rt) = @ka_code_typed checked(backend())(A, 1, 2, ndrange=1)
        # test that there is no invoke of overdub
        @test !any(check_for_overdub, CI.code)
    end
end


text:
export get_active_force
export update_active_force

"""
    get_active_force(active_brownian)

Returns the components of the propulsion force `γv` for active Brownian motion.
This function is used to avoid some type stability issues.
"""
@inline function get_active_force(active_brownian::ActiveBrownian)
    return active_brownian.γv_x, active_brownian.γv_y
end

"""
    update_active_force!(active_brownian, particle; dt)

Advances the active Brownian force by a timestep `dt`.  Currently this
requires `particle` as an input, which shouldn't be necessary.  This will be
removed in future versions.
"""
function update_active_force(active_brownian::ActiveBrownian; particle_orientation::Float64, dt::Float64)
    if active_brownian.align
        active_brownian.θ = particle_orientation
    else
        active_brownian.θ += active_brownian.amplitude * sqrt(dt) * randn()
    end
    sine, cosine = sincos(active_brownian.θ)
    active_brownian.γv_x = active_brownian.γv * cosine
    active_brownian.γv_y = active_brownian.γv * sine

    return particle_orientation
end

"""
    get_active_force(run_and_tumble)

Returns the components of the propulsion force `γv` for run-and-tumble
motion.  This function is used to avoid some type stability issues.
"""
@inline function get_active_force(run_and_tumble::RunAndTumble)
    return run_and_tumble.γv_x, run_and_tumble.γv_y
end

"""
    update_active_force!(active_brownian, particle; dt)

Advances the run-and-tumble force by a timestep `dt`.  Currently this
requires `particle` as an input, which shouldn't be necessary.  This will be
removed in future versions.
"""
function update_active_force(run_and_tumble::RunAndTumble; particle_orientation::Float64, dt::Float64)
    run_and_tumble.τ_run -= dt
    if run_and_tumble.τ_run < 0.0
        run_and_tumble.θ = 2 * pi * rand()
        run_and_tumble.τ_run = -log(rand()) / run_and_tumble.α

        if run_and_tumble.align
            particle_orientation = run_and_tumble.θ
        end

        sine, cosine = sincos(run_and_tumble.θ)
        run_and_tumble.γv_x = run_and_tumble.γv * cosine
        run_and_tumble.γv_y = run_and_tumble.γv * sine
    else
        if run_and_tumble.align
             run_and_tumble.θ = particle_orientation

             sine, cosine = sincos(run_and_tumble.θ)
             run_and_tumble.γv_x = run_and_tumble.γv * cosine
             run_and_tumble.γv_y = run_and_tumble.γv * sine
        end
    end

    return particle_orientation
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/05/2014

module CollisionAvoidanceSystem

export
    AbstractCollisionAvoidanceSystem,

    initialize,
    update,

    testThreat,
    selectRA,

    SimpleTCAS,
    SimpleTCASInput,
    SimpleTCASResolutionAdvisory,

    CoordSimpleTCAS,

    ACASX_CCAS,
    ACASX_ADD,
    ACASXInput,
    ACASXOutput

import CommonInterfaces: initialize, update

using AbstractCollisionAvoidanceSystemImpl

using SimpleTCASImpl

using CoordSimpleTCASImpl

using ACASXCommonImpl

pkgs = keys(Pkg.installed())
if "CCAS" in pkgs && "CASInterface" in pkgs
    @eval using ACASX_CCAS_Impl
end

using ACASX_ADD_Impl

end


text:
# GPU test environment
using Pkg
Pkg.add("Test")
Pkg.add("Pkg")
Pkg.add("CUDAdrv")
Pkg.add("CUDAnative")
Pkg.add("CuArrays")
Pkg.add("Adapt")

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
@reset
@set printstyle none

Econv = [-75.4196061687054,  -39.5638172309568,  -38.9378598408943]
Edf =   [-75.41960080317435, -39.56380446089254, -38.93785583198393]
@testset "UHF" begin
    @testset "Conventional Chonky" begin
        Fermi.Options.set("df", false)
        for i = eachindex(uhf_boys)
            # Read molecule
            path = joinpath(@__DIR__, "xyz/"*uhf_boys[i][1]*".xyz")
            mol = open(f->read(f,String), path)

            # Define options
            Fermi.Options.set("molstring", mol)
            Fermi.Options.set("basis", uhf_boys[i][2])
            Fermi.Options.set("reference", "uhf")
            Fermi.Options.set("charge", uhf_boys[i][3])
            Fermi.Options.set("multiplicity", uhf_boys[i][4])
            I = Fermi.Integrals.IntegralHelper(eri_type=Fermi.Integrals.Chonky())

            wf = @energy I => uhf
            @test isapprox(wf.energy, Econv[i], rtol=tol) # Energy from Psi4
        end
    end

    @testset "Conventional SparseERI" begin
        Fermi.Options.set("df", false)
        for i = eachindex(uhf_boys)
            # Read molecule
            path = joinpath(@__DIR__, "xyz/"*uhf_boys[i][1]*".xyz")
            mol = open(f->read(f,String), path)

            # Define options
            Fermi.Options.set("molstring", mol)
            Fermi.Options.set("basis", uhf_boys[i][2])
            Fermi.Options.set("reference", "uhf")
            Fermi.Options.set("charge", uhf_boys[i][3])
            Fermi.Options.set("multiplicity", uhf_boys[i][4])
            I = Fermi.Integrals.IntegralHelper(eri_type=Fermi.Integrals.SparseERI())

            wf = @energy I => uhf
            @test isapprox(wf.energy, Econv[i], rtol=tol) # Energy from Psi4
        end
    end

    @testset "Density Fitting" begin
        Fermi.Options.set("df", true)
        Fermi.Options.set("jkfit", "cc-pvtz-jkfit")
        for i in [1,3]
            # Read molecule
            path = joinpath(@__DIR__, "xyz/"*uhf_boys[i][1]*".xyz")
            mol = open(f->read(f,String), path)

            # Define options
            Fermi.Options.set("molstring", mol)
            Fermi.Options.set("basis", uhf_boys[i][2])
            Fermi.Options.set("reference", "uhf")
            Fermi.Options.set("charge", uhf_boys[i][3])
            Fermi.Options.set("multiplicity", uhf_boys[i][4])

            wf = @energy uhf
            @test isapprox(wf.energy, Edf[i], rtol=tol) # Energy from Psi4
        end
    end
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
#===============================================================================

A row of five black square tiles is to have a number of its tiles replaced with
coloured oblong tiles chosen from red (length two), green (length three), or
blue (length four).

If red tiles are chosen there are exactly seven ways this can be done. If green
tiles are chosen there are three ways. And if blue tiles are chosen there are
two ways.Assuming that colours cannot be mixed there are 7 + 3 + 2 = 12 ways of
replacing the black tiles in a row measuring five units in length.

How many different ways can the black tiles in a row measuring fifty units in
length be replaced if colours cannot be mixed and at least one coloured tile
must be used

===============================================================================#

using Memoize
using Base.Test

@memoize function ncases_swb(blocklen, rowlen)
    # Return the number of cases _s_tarting _w_ith a _b_lock

    if rowlen < blocklen
        return 0
    end

    result = 1
    # case for which there is only one block

    rowlen -= blocklen
    # remove room taken by the first block

    for nspace = 0:rowlen-blocklen
        # nbr spaces allowed between first block and second one

        result += ncases_swb(blocklen, rowlen-nspace)
    end

    result
end

function ncases(blocklen, rowlen)
    result = 0
    for nspace = 0:rowlen-blocklen
        # nbr spaces before first block
        result += ncases_swb(blocklen, rowlen-nspace)
    end

    result
end

function euler(rowlen)
    ncases(2, rowlen) + ncases(3, rowlen) + ncases(4, rowlen)
end

@test ncases(2, 5) == 7
@test ncases(3, 5) == 3
@test ncases(4, 5) == 2
@test euler(5) == 12

@printf("Euler 0116: %d", euler(50))

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
"""
    to_file(model::ModelChains, outpath::AbstractString)

This function writes the results of the MCMC runs into files. The destination of
the files is specified using `outpath`. It will create a files for each chain. A
`params_x.log` file storing each parameter sample. In this case `x` specifies the
index of the chain. The file is compatible with MCMC analysis tools like `Tracer`
(http://tree.bio.ed.ac.uk/software/tracer/). If in addition trees are sampled,
they are stored in newick format in a file called `trees_x.nwk`, where `x` again
specifies the index of the respective chain.
"""
function to_file(model::ModelChains, outpath::AbstractString)

    for run in 1:size(model.value,3)

        df = DataFrame(model.value[:,:,run])
        rename!(df, Symbol.(model.names))
        thin = model.range.step
        if isassigned(model.trees, 1)
            tdf = DataFrame(model.trees[:,:,run])
            to_file(df, tdf, outpath, string(run), thin, model.tree_names)
        else
            to_file(df, outpath, string(run), thin)
        end
    end

end

function to_file(df::DataFrame, outpath::AbstractString, run::AbstractString, thin::Int64)
    insertcols!(df,1, :it=>1:nrow(df))
    df[!, 1] .*= thin
    CSV.write(string(outpath, "params_"*run*".log"), df, writeheader=true, delim="\t")

end

function to_file(df::DataFrame, tdf::DataFrame, outpath::AbstractString,
                run::AbstractString, thin::Int64, tree_names::Array{A}) where A <: AbstractString

    insertcols!(df,1, :it=>1:nrow(df))
    df[!, 1] .*= thin
    CSV.write(string(outpath, "params_"*run*".log"), df, header=true, delim="\t")
    for (ind, n_name) in enumerate(tree_names)
        io = open(string(outpath, "trees_"*n_name*"_"*run*".nwk"), "w")
        for x = 1:length(tdf[:,ind])
            write(io, tdf[x,ind])
            write(io, "\n")
        end
        close(io)
    end
end


text:
# Do not share a stream between processes
# The token would be shared so putting would give InvalidSequenceTokenException a lot
struct CloudWatchLogHandler{F<:Formatter} <: Handler{F, Union{}}
    stream::CloudWatchLogStream
    channel::Channel{LogEvent}  # only one task should read from this channel
    fmt::F
end

"""
    CloudWatchLogHandler(
        config::AWSConfig,
        log_group_name,
        log_stream_name,
        formatter::Memento.Formatter,
    )

Construct a Memento Handler for logging to a CloudWatch Log Stream.
This constructor creates a task which asynchronously submits logs to the stream.

A CloudWatch Log Event has only two properties: `timestamp` and `message`.

If a `Record` has a `date` property it will be used as the `timestamp`, otherwise the
current time will be captured when `Memento.emit` is called.
All `DateTime`s will be assumed to be in UTC.

The `message` will be generated by calling `Memento.format` on the `Record` with this
handler's `formatter`.
"""
function CloudWatchLogHandler(
    config::AWSConfig,
    log_group_name::AbstractString,
    log_stream_name::AbstractString,
    formatter::F=DefaultFormatter(),
) where F<:Formatter
    ch = Channel{LogEvent}(Inf)
    handler = CloudWatchLogHandler(
        CloudWatchLogStream(config, log_group_name, log_stream_name),
        ch,
        formatter,
    )

    tsk = @async process_logs!(handler)
    # channel will be closed if task fails, to avoid unknowingly discarding logs
    bind(ch, tsk)

    return handler
end

function process_available_logs!(handler::CloudWatchLogHandler)
    events = Vector{LogEvent}()
    batch_size = 0

    while isready(handler.channel) && length(events) < MAX_BATCH_LENGTH
        event = fetch(handler.channel)
        batch_size += aws_size(event)
        if batch_size <= MAX_BATCH_SIZE
            take!(handler.channel)
            push!(events, event)
        else
            break
        end
    end

    if isempty(events)
        warn(LOGGER, string(
            "Channel was ready but no events were found. ",
            "Is there another task pulling logs from this handler?",
        ))
    end

    try
        @mock submit_logs(handler.stream, events)
    catch e
        warn(LOGGER, CapturedException(e, catch_backtrace()))
    end
end

"""
    process_logs!(handler::CloudWatchLogHandler)

Continually pulls logs from the handler's channel and submits them to AWS.
This function should terminate silently when the channel is closed.
"""
function process_logs!(handler::CloudWatchLogHandler)
    group = handler.stream.log_group_name
    stream = handler.stream.log_stream_name

    debug(LOGGER, "Handler for group '$group' stream '$stream' initiated")

    try
        while isopen(handler.channel)  # might be able to avoid the error in this case
            wait(handler.channel)
            process_available_logs!(handler)
            sleep(PUTLOGEVENTS_RATE_LIMIT)  # wait at least this long due to AWS rate limits
        end
    catch err
        if !(err isa InvalidStateException && err.state === :closed)
            log(
                LOGGER,
                :error,
                "Handler for group '$group' stream '$stream' terminated unexpectedly",
            )
            error(LOGGER, CapturedException(err, catch_backtrace()))
        end
    end

    debug(LOGGER, "Handler for group '$group' stream '$stream' terminated normally")

    return nothing
end

function Memento.emit(handler::CloudWatchLogHandler, record::Record)
    dt = isdefined(record, :date) ? record.date : Dates.now(tz"UTC")
    message = format(handler.fmt, record)
    event = LogEvent(message, dt)
    put!(handler.channel, event)
end


text:
using SEAL
using Printf

function print_example_banner(title)
  if isempty(title)
    return nothing
  end

  title_length = length(title)
  banner_length = title_length + 2 * 10
  banner_border = "+" * "-"^(banner_length - 2) * "+"
  banner_center = "|" * " "^9 * title * " "^9 * "|"

  println()
  println(banner_border)
  println(banner_center)
  println(banner_border)
end

function print_parameters(context::SEALContext)
  context_data = key_context_data(context)
  encryption_parms = parms(context_data)
  scheme_type = scheme(encryption_parms)
  if scheme_type == SchemeType.bfv
    scheme_name = "BFV"
  elseif scheme_type == SchemeType.ckks
    scheme_name = "CKKS"
  else
    error("unsupported scheme")
  end

  println("/")
  println("| Encryption parameters :")
  println("|   scheme: ", scheme_name)
  println("|   poly_modulus_degree: ", poly_modulus_degree(encryption_parms))

  print("|   coeff_modulus size: ", total_coeff_modulus_bit_count(context_data), " (")
  bit_counts = [bit_count(modulus) for modulus in coeff_modulus(encryption_parms)]
  print(join(bit_counts, " + "))
  println(") bits")

  if scheme_type == SchemeType.bfv
    println("|   plain_modulus: ", value(plain_modulus(encryption_parms)))
  end

  println("\\")
end

function print_line(line_number)
  if line_number < 10
    padding = "  "
  elseif line_number < 100
    padding = " "
  else
    padding = ""
  end
  print("Line ", padding, line_number, " --> ")
end

function print_vector(vector, print_size=4, prec=3)
  slot_count = length(vector)
  println()
  if slot_count <= 2 * print_size
    print("    [ ")
    if prec == 7
      @printf "%.7f" first(vector)
    else
      @printf "%.3f" first(vector)
    end
    for i in 2:slot_count
      if prec == 7
        @printf ", %.7f" vector[i]
      else
        @printf ", %.3f" vector[i]
      end
    end
    println(" ]")
  else
    print("    [ ")
    if prec == 7
      @printf "%.7f" first(vector)
    else
      @printf "%.3f" first(vector)
    end
    for i in 2:print_size
      if prec == 7
        @printf ", %.7f" vector[i]
      else
        @printf ", %.3f" vector[i]
      end
    end
    if slot_count > 2 * print_size
      print(", ...")
    end
    for i in (slot_count - print_size + 1):slot_count
      if prec == 7
        @printf ", %.7f" vector[i]
      else
        @printf ", %.3f" vector[i]
      end
    end
    println(" ]")
  end
  println()
end

function print_matrix(matrix, row_size)
  print_size = 5

  println()

  # First row
  print("    [")
  for i in 1:print_size
    @printf("%3d,", matrix[i])
  end
  print(" ...,")
  for i in (row_size - print_size + 1):row_size
    if i != row_size
      @printf("%3d,", matrix[i])
    else
      @printf("%3d ]\n", matrix[i])
    end
  end

  # Second row
  print("    [")
  for i in (row_size + 1):(row_size + print_size)
    @printf("%3d,", matrix[i])
  end
  print(" ...,")
  for i in (2 * row_size - print_size + 1):(2 * row_size)
    if i != 2 * row_size
      @printf("%3d,", matrix[i])
    else
      @printf("%3d ]\n", matrix[i])
    end
  end

  println()
end

"""
    @elapsedus(ex)

Return integer number of elapsed microseconds required for executing `ex`. Modified from
`Base.@elapsed`, which returns the number of seconds as a floating point number.

See also: [`@elapsed`](@ref)
"""
macro elapsedus(ex)
  quote
    while false; end # compiler heuristic: compile this block (alter this if the heuristic changes)
    local t0 = time_ns()
    $(esc(ex))
    Int(div(time_ns() - t0, 1000))
  end
end


text:
### accessing libc functions (XOPEN_SOURCE >= 700, POSIX_C_SOURCE >= 200809L glibc>=2.24)
using ResourceBundles.CLocales
using .LC

import .CLocales: newlocale_c, strcoll_c, nl_langinfo_c

const P0 = Ptr{Nothing}(0)

if Sys.islinux() && get(Base.ENV, "NO_CLOCALE", "") != "1"

@test newlocale_c(LC._MASK_ALL, "invalidxxx", P0)  == P0
@test newlocale_c(LC._MASK_ALL, "en_US.utf8", P0)  != P0

test_locale_C = newlocale_c(LC._MASK_ALL, "C", P0)
test_locale_ca = newlocale_c(LC._MASK_ALL, "en_US.utf8", P0)
@test duplocale(test_locale_ca) != P0

@test unsafe_string(nl_langinfo_c(Cint(0xffff), test_locale_ca)) == "en_US.utf8"

COLL_TESTS_C = [
    ( "a", "b", -1 ),
    ( "A", "b", -1 ),
    ( "a", "B", +1 ),
]

COLL_TESTS_FR = [
    ( "a", "b", -1 ),
    ( "A", "b", -1 ),
    ( "a", "B", -1 ),
]

@testset "string comparisons '$a' ~ '$b'" for (a, b, r) in COLL_TESTS_C
    @test sign(strcoll_c(a, b, test_locale_C)) == r
end

@testset "string comparisons '$a' ~ '$b'" for (a, b, r) in COLL_TESTS_FR
    @test sign(strcoll_c(a, b, test_locale_ca)) == r
end

@test freelocale(test_locale_ca) == nothing
freelocale(test_locale_C)
end
#################################


text:
abstract type Loss end

"""
    get_loss(l, y, ŷ)

Return a loss.
"""
get_loss(l::Loss, y, ŷ) = throw("unimplemented")

struct LogLoss <: Loss
end

get_loss(l::LogLoss, y, ŷ) = logloss(y, ŷ)

struct MSE <: Loss
end

get_loss(l::MSE, y, ŷ) = mse(y, ŷ)

max_stars_repo_name: taraspiotr/DeepJulia

max_issues_repo_name: taraspiotr/DeepJulia

max_forks_repo_name: taraspiotr/DeepJulia


text:
# ---
# title: 827. Making A Large Island
# id: problem827
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search
# link: <https://leetcode.com/problems/making-a-large-island/description/>
# hidden: true
# ---
#
# In a 2D grid of `0`s and `1`s, we change at most one `0` to a `1`.
#
# After, what is the size of the largest island? (An island is a 4-directionally
# connected group of `1`s).
#
# **Example 1:**
#
#
#
#     Input: [[1, 0], [0, 1]]
#     Output: 3
#     Explanation: Change one 0 to 1 and connect two 1s, then we get an island with area = 3.
#
#
# **Example 2:**
#
#
#
#     Input: [[1, 1], [1, 0]]
#     Output: 4
#     Explanation: Change the 0 to 1 and make the island bigger, only one island with area = 4.
#
# **Example 3:**
#
#
#
#     Input: [[1, 1], [1, 1]]
#     Output: 4
#     Explanation: Can't change any 0 to 1, only one island with area = 4.
#
#
#
# Notes:
#
#   * `1 <= grid.length = grid[0].length <= 50`.
#   * `0 <= grid[i][j] <= 1`.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end