text:
# io = open("C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/train-images-idx3-ubyte/train-images.idx3-ubyte", n=28)
#write(io)

#open("C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/train-images-idx3-ubyte/train-images.idx3-ubyte") do file
#    for ln in eachline(file)
#        println("$(length(ln)), $(ln)")
#    end
#end

# Back propogation neural network
# 'Hello World' of neural networks
# Andrew Bannerman 4.21.2019

# Resources
# Cheat Sheet = https://ml-cheatsheet.readthedocs.io/en/latest/forwardpropagation.html
# Matrix multiplication Rules - https://www.mathsisfun.com/algebra/matrix-multiplying.html
# How to get MNIST data as a CSV https://pjreddie.com/projects/mnist-in-csv/
# Kaggle - https://www.kaggle.com/ngbolin/mnist-dataset-digit-recognizer

# Back propogation tutorial: https://dev.to/shamdasani/build-a-flexible-neural-network-with-backpropagation-in-python

# Needed packages
using DataFrames
using CSV
using Distributions
using Plots
using Random

#=
Random.seed!(123) # Setting the seed
d = Normal()
x = rand(d, 100)
x = rand(Normal(100, 2), 100)
vars = var(x)
sqrt(vars)
histogram(x)
=#

# Train and test set data preparation

# Load MNIST Data set from .csv format
train_set = CSV.read("C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/mnist_train.csv",header=false,types=fill(Float64,785))
test_set = CSV.read("C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/mnist_test.csv",header=false,types=fill(Float64,785))

# Plot Example image
example = Float64.(test_set[2,2:size(test_set,2)])
example = rotl90(reshape(example,28,28))
heatmap(example,c = :greys)
print(example)

# Extract labels
train_label = train_set[:,1]
test_label = test_set[:,1]

# Extract Data and exclude the label in the first column
train_set = train_set[:,2:size(train_set,2)]
test_set = test_set[:,2:size(test_set,2)]

# Pre process - Scale the inputs to a range of 0.01 to 1.00
m_train = Array{Float64}(zeros(nrow(train_set),size(train_set,2))) # empty matrix
@inbounds for i in 1:nrow(train_set)
    @inbounds for j in 1:size(train_set,2)
    m_train[i,j] = (train_set[i,j] / 255.0 * 0.99) + 0.01
    end
end

m_test = Array{Float64}(zeros(nrow(test_set),size(test_set,2))) # empty matrix
@inbounds for i in 1:nrow(test_set)
    @inbounds for j in 1:size(test_set,2)
    m_test[i,j] = (test_set[i,j] / 255.0 * 0.99) + 0.01
    end
end

# Convert matrix to DataFrame and append the labels back to the train and test set
scaled_orig_train_inputs = DataFrame(hcat(train_label,m_train))
scaled_test_inputs = DataFrame(hcat(test_label,m_test))

# Plot the new scale example
example = Float64.(scaled_orig_train_inputs[1,2:size(scaled_orig_train_inputs,2)])
example = rotl90(reshape(example,28,28))
heatmap(example, c = :greys)

#---------
# Data Augmentation
#---------

# Rotate images
training_set_aug_minus = zeros(size(scaled_orig_train_inputs,1),size(scaled_orig_train_inputs,2))
training_set_aug_plus = zeros(size(scaled_orig_train_inputs,1),size(scaled_orig_train_inputs,2))
degree_rot_minus = -.3
degree_rot_plus = .3
i=1
for i in 1:size(scaled_orig_train_inputs,1)
        label = scaled_orig_train_inputs[i,1]
        example = Float64.(scaled_orig_train_inputs[i,2:size(scaled_orig_train_inputs,2)])
        example = reshape(example,28,28)
        #rot_minus = LinearMap(RotMatrix(degree_rot_minus))
        #rot_plus = LinearMap(RotMatrix(degree_rot_plus))
        #imgw_minus = warp(example_minus, rot_minus, axes(example_plus)) #  axes(example_minus) = crops
        #imgw_plus = warp(example_plus, rot_plus, axes(example_plus)) # axes() = crops
        imgw_minus = imrotate(example, degree_rot_minus, axes(example))
        imgw_plus = imrotate(example, degree_rot_plus, axes(example))
        imgw_minus = reshape(imgw_minus,1,784)
        imgw_plus = reshape(imgw_plus,1,784)
            for k in 1:size(imgw_minus,2)
                if isnan(imgw_minus[k]) # Convert NaN to 0.0
                    imgw_minus[k] = 0.01
                end
                if isnan(imgw_plus[k]) # Convert NaN to 0.0
                    imgw_plus[k] = 0.01
                end
            end
            #imgw_plus .= (imgw_plus ./ 255.0 .* 0.99) .+ 0.01 # Scale the outputs
            #imgw_minus .= (imgw_minus ./ 255.0 .* 0.99) .+ 0.01 # Scale the outputs
            training_set_aug_minus[i,1] = label
            training_set_aug_minus[i,2:size(scaled_orig_train_inputs,2)] = imgw_minus
            training_set_aug_plus[i,1] = label
            training_set_aug_plus[i,2:size(scaled_orig_train_inputs,2)] = imgw_plus
        print("This is iteration i",i,"\n")
end

    maximum(imgw_plus)

# Combine augmented data with the original
# note may want to add shuffling?
scaled_train_inputs = DataFrame(vcat(scaled_orig_train_inputs,training_set_aug_minus,training_set_aug_plus))
#scaled_train_inputs = DataFrame(training_set_aug_plus)

# Shuffle the augmented + original data
scaled_train_inputs = scaled_train_inputs[shuffle(1:end), :]

#---------
# Activation Functions
#---------

function sigmoid(x)
    return 1 / (1.0 + exp(-x))
end

function sigmoid_derivative(x)
    s = 1 / (1 + exp(-x))
    return s * (1 - s)
end

function ReLU(x)
    #return x*(x>0)
    return maximum([0.0,x])
end

function ReLU_derivative(x)
    if x <= 0
    end
    return 0
end

#---------
# Layers
# https://medium.freecodecamp.org/building-a-3-layer-neural-network-from-scratch-99239c4af5d3
#---------

function dropout(d_in::Array{Float64,2}; p::Float64=.8)
    out = (d_in .* rand(Binomial(1,p), size(d_in)[1])) ./ p
return out
end

activation = "sigmoid"
drop_out = false
function forward_prop(W1::Array{Float64,2}, b1::Array{Float64,2}, W2::Array{Float64,2}, b2::Array{Float64,2}, X::Array{Float64,2}, Y::Array{Float64,2}; drop_out::Bool=true, p::Float64=.5, activation::String="sigmoid")
    Z1 = W1 * X + b1
    if drop_out == false
        if activation == "sigmoid"
    A1 = sigmoid.(Z1)
elseif activation == "relu"
    A1 = ReLU.(Z1)
    end
    elseif drop_out == true
        if activation == "sigmoid"
    A1 = sigmoid.(Z1)
    A1 = dropout(A1,p=.8)
    elseif activation == "relu"
    A1 = ReLU.(Z1)
    A1 = dropout(A1,p=.8) # Apply drop out by hidden output * 0,1's / p (1/p - scaled so we dont need to apply at test time)(inverted dropout)
        end
    end
    Z2 = (W2 * A1) + b2
    if activation == "sigmoid"
    A2 = sigmoid.(Z2)
elseif activation == "relu"
    A2 = ReLU.(Z2)
    end
    return Z1, A1, Z2, A2
end

#----
# back prop math: http://deeplizard.com/learn/video/G5b4jRBKNxw
# https://github.com/mbadry1/DeepLearning.ai-Summary/tree/master/1-%20Neural%20Networks%20and%20Deep%20Learning#gradient-descent-on-m-examples
#----
function backward_prop(W1::Array{Float64,2}, b1::Array{Float64,2}, W2::Array{Float64,2}, b2::Array{Float64,2}, X::Array{Float64,2}, Y::Array{Float64,2}, Z1::Array{Float64,2}, A1::Array{Float64,2}, Z2::Array{Float64,2}, A2::Array{Float64,2}; learning_rate::Float64=.01,activation::String="relu")
    dZ2 = A2 - Y # output errors
    dW2 = (dZ2 * A1')
    db2 =  sum(dZ2, dims = 2)
    if activation == "sigmoid"
        dZ1 = (W2' * dZ2) .* sigmoid_derivative.(A1) #(1 .- A1.^2)
    elseif activation == "relu"
        dZ1 = (W2' * dZ2) .* ReLU_derivative.(A1) #(1 .- A1.^2)
    end
    dW1 = (dZ1 * X')
    db1 =  sum(dZ1,dims = 2)
   # Update params
    W1 = W1 .- (learning_rate .* dW1)
     b1 = b1 .- (learning_rate .* db1)
     W2 = W2 .- (learning_rate * dW2)
     b2 = b2 .- (learning_rate .* db2)
    # print("check w1 is updating",W2[1])

    return  dW1, db1, dW2, db2, W1, b1, W2, b2
end

#---------
# Cost Function
#---------

function cost_function(A2::Array{Float64,2}, Y::Array{Float64,2}; activation::String="sigmoid")
    if activation == "sigmoid"
    m = size(Y)[1]
    #cost = (- 1 / m) * np.sum(Y * np.log(A) + (1 - Y) * (np.log(1 - A)))  # compute cost
    cost = (-1/m) * sum((Y' * log.(A2)) .+ ((1 .- Y')* log.(1 .- A2)))
        #cost = (-1/m) * sum((Y * log.(A2')) .+ ((1 .- Y)* log.(1 .- A2')))
elseif activation == "relu"
    return Inf
end
return cost
end

#---------
# Initialize network
#---------

# Define neural network
# Number of input, hidden and output nodes
input_nodes = size(train_set,2) # equal to the size of the 28 * 28 image (784 elements)
hidden_nodes = 256
output_nodes = 10 # Equal to the number of labels in this case numbers 0 to 9 (10 elements total)

# Create instance of neural network
inodes = input_nodes
hnodes = hidden_nodes
onodes = output_nodes
# Weight initialization
# w = np.random.randn(n) * sqrt(2.0/n) for relu
# w = np.random.randn(n) / sqrt(n) # other
#W1 = randn(hnodes,inodes) * inodes^(-0.5)
#W2 = randn(onodes,hnodes) * hnodes^(-0.5)

W1 = randn(hnodes,inodes) * sqrt(2.0/inodes)
W2 = randn(onodes,hnodes) * sqrt(2.0/hnodes)
# Initializing the biases as 0's
b1 = zeros(hnodes,1)
b2 = zeros(onodes,1)

epochs = 5
# Check that the neural network is converging
# Initialize training progress
learning_curve = fill(0,nrow(scaled_train_inputs))
training_accuracy = fill(0.0,nrow(scaled_train_inputs))
i = 1
j = 1
@inbounds for i in 1:epochs
    # go through all records in the training data set
    @inbounds for j in 1:nrow(scaled_train_inputs)
        # Forward propogation
        label = Int64.(scaled_train_inputs[j,1])
        X = reshape(Float64.(scaled_train_inputs[j,2:size(scaled_train_inputs,2)]),784,1)
        # create the target output values (all 0.01, except the desired label which is 0.99)
        Y = zeros(10,1)
        Y .= 0.01 .+ Y
        # Set the target .99 to the correct index position within the target array
        Y[label+1] = 0.99 # We +1 because the number convention starts at 0 to 9. Thus array position 6 is actually target number 5.
global Z1, A1, Z2, A2 = forward_prop(W1, b1, W2, b2, X, Y; drop_out=false, p=.5,activation="relu")
cost_rate = cost_function(A2,Y;activation="relu")
global dW1, db1, dW2, db2, W1, b1, W2, b2 = backward_prop(W1, b1, W2, b2, X, Y, Z1, A1, Z2, A2; learning_rate=.01,activation="relu")
#print("check training weight being update", dW1[1],"\n")
network_train_output = argmax(A2)[1]-1 # adjust index position to account from 0 start
if (label == network_train_output)
    # network's answer matches correct answer, add 1 to learning curve
    learning_curve[j] = 1
else
    # network's answer doesn't match correct answer, add 0 to learning curve
    learning_curve[j] = 0
end
if (j == 1)
    global learning_curve = fill(0,nrow(scaled_train_inputs))
end
print("Cost after iteration ",j," ",cost_rate,"\n")
print("Training Accuracy ",(sum(learning_curve) / j),"\n")
print("This is epoch ",i,"\nThis is iteration ",j,"\n")

    end
end

sum(learning_curve) / 60000

print(learning_curve)

# Plot the learning curve
plot(cumsum(learning_curve),title="Learning Curve")

#---------
# Test Function
#---------

# Function for applying the network to unseen or test data not seen (trained) by the network before
function test_nn(W1::Array{Float64,2}, W2::Array{Float64,2},b1::Array{Float64,2}, b2::Array{Float64,2},X::Array{Float64,2})
    Z1 = W1 * X + b1
        if activation == "sigmoid"
    A1 = sigmoid.(Z1)
elseif activation == "relu"
    A1 = ReLU.(Z1)
    end
    Z2 = (W2 * A1) + b2
    if activation == "sigmoid"
    A2 = sigmoid.(Z2)
elseif activation == "relu"
    A2 = ReLU.(Z2)
    end
    return A2
end

#---------
# Test Network on Individual image
#---------

# Call the network on an indivdual entry
example_call = Float64.(scaled_test_inputs[500,2:size(scaled_test_inputs,2)])
example_call = reshape(example_call,784,1)
example_call_image = rotl90(reshape(example_call,28,28))
# Plot Number
heatmap(example_call_image,c = :greys)
# Grab the known correct label
correct_label = scaled_test_inputs[500,1]
# Apply the trained network to the input example
outputs = test_nn(W1, W2, b1,b2,example_call)
# Check what the label thought the number was
network_label = argmax(outputs)[1]-1
if network_label ==  correct_label
    print("The neural network correctly identified the hand written number - ", network_label)
else
    print("The neural network incorrectly identified the hand written number")
end

print(example_call_image)

#---------
# Test Neural Network on Test Data
#---------

# Initialize the scorecard output the same size as the number of test set labels
scorecard = fill(0,nrow(scaled_test_inputs))
test_accuracy = fill(0.0,nrow(scaled_test_inputs))

# Input the test set data over the trained neural network to see how well it does
i=1
@inbounds for i in 1:nrow(scaled_test_inputs)
    let wih = wih, who = who
    correct_label = Int64.(scaled_test_inputs[i,1])
    inputs = reshape(Float64.(scaled_test_inputs[i,2:size(scaled_test_inputs,1)]),784,1) # rotr90() for changing dim to [784,1] from [1,784] to meet matrix multiplication rules
    # query the network
    outputs = test_nn(W1, W2, b1,b2,inputs) # Test the network over the test set data (wih and who already set during training - no weight updating happens this time round)
    # Find element position of what the network thinks the output is ie if output is a .96 in element position 6 then the network thinks the label was 5.
    # This is because we have hand written numbers 0 to 9, total 10 outputs. Starting from index position 1 = label 0 hence index position 6 then = label 5
    label = argmax(outputs)[1]-1 # adjust index position to account from 0 start
    # append correct or incorrect to list
    if (label == correct_label)
        # network's answer matches correct answer, add 1 to scorecard
        scorecard[i] = 1
    else
        # network's answer doesn't match correct answer, add 0 to scorecard
        scorecard[i] = 0
        end
        print("Accuracy ",sum(scorecard) / i,"\n")
        test_accuracy[i] = sum(scorecard) / i
    end
end

# Percentage Correct
percentage = sum(scorecard) / size(scorecard,1)

# Change log
# 5.12.2019 - Changed the back propogation derivatives to allow different activation functions
# 5.12.2019 - Added bias
# 4.23.2019 - Added first implemntation of drop out
# 4.22.2019 - added a % print function to check train / test accuracy


text:
 __precompile__()

module Plots

using Compat
using Reexport
@reexport using Colors

export
  Plot,
  Subplot,
  SubplotLayout,
  GridLayout,
  FlexLayout,
  AVec,
  AMat,

  plot,
  plot!,
  # plot_display,
  # plot_display!,
  subplot,
  subplot!,

  current,
  default,
  with,

  scatter,
  scatter!,
  bar,
  bar!,
  histogram,
  histogram!,
  density,
  density!,
  heatmap,
  heatmap!,
  hexbin,
  hexbin!,
  sticks,
  sticks!,
  hline,
  hline!,
  vline,
  vline!,
  ohlc,
  ohlc!,
  pie,
  pie!,
  contour,
  contour!,
  surface,
  surface!,
  wireframe,
  wireframe!,
  path3d,
  path3d!,
  scatter3d,
  scatter3d!,

  title!,
  xlabel!,
  ylabel!,
  xlims!,
  ylims!,
  xticks!,
  yticks!,
  annotate!,
  xflip!,
  yflip!,
  xaxis!,
  yaxis!,

  savefig,
  png,
  gui,

  backend,
  backends,
  aliases,
  dataframes,

  Shape,
  text,
  font,
  stroke,
  brush,
  Surface,
  OHLC,

  colorscheme,
  ColorScheme,
  ColorGradient,
  ColorVector,
  ColorWrapper,
  ColorFunction,
  ColorZFunction,
  getColor,
  getColorZ,

  debugplots,

  supportedArgs,
  supportedAxes,
  supportedTypes,
  supportedStyles,
  supportedMarkers,
  subplotSupported,

  Animation,
  frame,
  gif,

  # recipes
  PlotRecipe,
  # EllipseRecipe,
  # spy,
  corrplot

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

const IMG_DIR = Pkg.dir("Plots") * "/img/"

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

include("types.jl")
include("utils.jl")
include("colors.jl")
include("components.jl")
include("plotter.jl")
include("args.jl")
include("plot.jl")
include("subplot.jl")
include("recipes.jl")
include("animation.jl")
include("output.jl")

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

scatter(args...; kw...)    = plot(args...; kw...,  linetype = :scatter)
scatter!(args...; kw...)   = plot!(args...; kw..., linetype = :scatter)
bar(args...; kw...)        = plot(args...; kw...,  linetype = :bar)
bar!(args...; kw...)       = plot!(args...; kw..., linetype = :bar)
histogram(args...; kw...)  = plot(args...; kw...,  linetype = :hist)
histogram!(args...; kw...) = plot!(args...; kw..., linetype = :hist)
density(args...; kw...)    = plot(args...; kw...,  linetype = :density)
density!(args...; kw...)   = plot!(args...; kw..., linetype = :density)
heatmap(args...; kw...)    = plot(args...; kw...,  linetype = :heatmap)
heatmap!(args...; kw...)   = plot!(args...; kw..., linetype = :heatmap)
hexbin(args...; kw...)     = plot(args...; kw...,  linetype = :hexbin)
hexbin!(args...; kw...)    = plot!(args...; kw..., linetype = :hexbin)
sticks(args...; kw...)     = plot(args...; kw...,  linetype = :sticks, marker = :ellipse)
sticks!(args...; kw...)    = plot!(args...; kw..., linetype = :sticks, marker = :ellipse)
hline(args...; kw...)      = plot(args...; kw...,  linetype = :hline)
hline!(args...; kw...)     = plot!(args...; kw..., linetype = :hline)
vline(args...; kw...)      = plot(args...; kw...,  linetype = :vline)
vline!(args...; kw...)     = plot!(args...; kw..., linetype = :vline)
ohlc(args...; kw...)       = plot(args...; kw...,  linetype = :ohlc)
ohlc!(args...; kw...)      = plot!(args...; kw..., linetype = :ohlc)
pie(args...; kw...)        = plot(args...; kw...,  linetype = :pie)
pie!(args...; kw...)       = plot!(args...; kw..., linetype = :pie)
contour(args...; kw...)    = plot(args...; kw...,  linetype = :contour)
contour!(args...; kw...)   = plot!(args...; kw..., linetype = :contour)
surface(args...; kw...)    = plot(args...; kw...,  linetype = :surface)
surface!(args...; kw...)   = plot!(args...; kw..., linetype = :surface)
wireframe(args...; kw...)  = plot(args...; kw...,  linetype = :wireframe)
wireframe!(args...; kw...) = plot!(args...; kw..., linetype = :wireframe)
path3d(args...; kw...)     = plot(args...; kw...,  linetype = :path3d)
path3d!(args...; kw...)    = plot!(args...; kw..., linetype = :path3d)
scatter3d(args...; kw...)  = plot(args...; kw...,  linetype = :scatter3d)
scatter3d!(args...; kw...) = plot!(args...; kw..., linetype = :scatter3d)

title!(s::@compat(AbstractString); kw...)                 = plot!(; title = s, kw...)
xlabel!(s::@compat(AbstractString); kw...)                = plot!(; xlabel = s, kw...)
ylabel!(s::@compat(AbstractString); kw...)                = plot!(; ylabel = s, kw...)
xlims!{T<:Real,S<:Real}(lims::@compat(Tuple{T,S}); kw...) = plot!(; xlims = lims, kw...)
ylims!{T<:Real,S<:Real}(lims::@compat(Tuple{T,S}); kw...) = plot!(; ylims = lims, kw...)
xlims!(xmin::Real, xmax::Real; kw...)                     = plot!(; xlims = (xmin,xmax), kw...)
ylims!(ymin::Real, ymax::Real; kw...)                     = plot!(; ylims = (ymin,ymax), kw...)
xticks!{T<:Real}(v::AVec{T}; kw...)                       = plot!(; xticks = v, kw...)
yticks!{T<:Real}(v::AVec{T}; kw...)                       = plot!(; yticks = v, kw...)
xticks!{T<:Real,S<:@compat(AbstractString)}(
              ticks::AVec{T}, labels::AVec{S}; kw...)     = plot!(; xticks = (ticks,labels), kw...)
yticks!{T<:Real,S<:@compat(AbstractString)}(
              ticks::AVec{T}, labels::AVec{S}; kw...)     = plot!(; yticks = (ticks,labels), kw...)
annotate!(anns...; kw...)                                 = plot!(; annotation = anns, kw...)
annotate!{T<:Tuple}(anns::AVec{T}; kw...)                 = plot!(; annotation = anns, kw...)
xflip!(flip::Bool = true; kw...)                          = plot!(; xflip = flip, kw...)
yflip!(flip::Bool = true; kw...)                          = plot!(; yflip = flip, kw...)
xaxis!(args...; kw...)                                    = plot!(; xaxis = args, kw...)
yaxis!(args...; kw...)                                    = plot!(; yaxis = args, kw...)

title!(plt::Plot, s::@compat(AbstractString); kw...)                  = plot!(plt; title = s, kw...)
xlabel!(plt::Plot, s::@compat(AbstractString); kw...)                 = plot!(plt; xlabel = s, kw...)
ylabel!(plt::Plot, s::@compat(AbstractString); kw...)                 = plot!(plt; ylabel = s, kw...)
xlims!{T<:Real,S<:Real}(plt::Plot, lims::@compat(Tuple{T,S}); kw...)  = plot!(plt; xlims = lims, kw...)
ylims!{T<:Real,S<:Real}(plt::Plot, lims::@compat(Tuple{T,S}); kw...)  = plot!(plt; ylims = lims, kw...)
xlims!(plt::Plot, xmin::Real, xmax::Real; kw...)                      = plot!(plt; xlims = (xmin,xmax), kw...)
ylims!(plt::Plot, ymin::Real, ymax::Real; kw...)                      = plot!(plt; ylims = (ymin,ymax), kw...)
xticks!{T<:Real}(plt::Plot, ticks::AVec{T}; kw...)                    = plot!(plt; xticks = ticks, kw...)
yticks!{T<:Real}(plt::Plot, ticks::AVec{T}; kw...)                    = plot!(plt; yticks = ticks, kw...)
xticks!{T<:Real,S<:@compat(AbstractString)}(plt::Plot,
                          ticks::AVec{T}, labels::AVec{S}; kw...)     = plot!(plt; xticks = (ticks,labels), kw...)
yticks!{T<:Real,S<:@compat(AbstractString)}(plt::Plot,
                          ticks::AVec{T}, labels::AVec{S}; kw...)     = plot!(plt; yticks = (ticks,labels), kw...)
annotate!(plt::Plot, anns...; kw...)                                  = plot!(plt; annotation = anns, kw...)
annotate!{T<:Tuple}(plt::Plot, anns::AVec{T}; kw...)                  = plot!(plt; annotation = anns, kw...)
xflip!(plt::Plot, flip::Bool = true; kw...)                           = plot!(plt; xflip = flip, kw...)
yflip!(plt::Plot, flip::Bool = true; kw...)                           = plot!(plt; yflip = flip, kw...)
xaxis!(plt::Plot, args...; kw...)                                     = plot!(plt; xaxis = args, kw...)
yaxis!(plt::Plot, args...; kw...)                                     = plot!(plt; yaxis = args, kw...)

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

# try
#   import DataFrames
#   dataframes()
# end

# const CURRENT_BACKEND = pickDefaultBackend()

# for be in backends()
#   try
#     backend(be)
#     backend()
#   catch err
#     @show err
#   end
# end

const CURRENT_BACKEND = CurrentBackend(:none)

# function __init__()
#   # global const CURRENT_BACKEND = pickDefaultBackend()
#   # global const CURRENT_BACKEND = CurrentBackend(:none)

#   # global CURRENT_BACKEND
#   # println("[Plots.jl] Default backend: ", CURRENT_BACKEND.sym)

#   # # auto init dataframes if the import statement doesn't error out
#   # try
#   #   @eval import DataFrames
#   #   dataframes()
#   # end
# end

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

end # module

ext: jl


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# knetFile *knet_open(const char *fn, const char *mode);
@defun knet_open(fn::Cstring, mode::Cstring)::Ptr{knetFile}

# knetFile *knet_dopen(int fd, const char *mode);
@defun knet_dopen(fd::Cint, mode::Cstring)::Ptr{knetFile}

# ssize_t knet_read(knetFile *fp, void *buf, size_t len);
@defun knet_read(fp::Ptr{knetFile}, buf::Ptr{Cvoid}, len::Csize_t)::Cssize_t

# off_t knet_seek(knetFile *fp, off_t off, int whence);
@defun knet_seek(fp::Ptr{knetFile}, off::Coff_t, whence::Cint)::Coff_t

# int knet_close(knetFile *fp);
@defun knet_close(fp::Ptr{knetFile})::Cint


text:
module VersionCheck

using Pkg, Logging, Dates, Random
using JSON3, UrlDownload
using Scratch

const version_info_bounds = r"start-versions-->(.*)<!--end-versions"s
const usersettings_filename = "usersettings.json"

# changelog_url = "https://genieframework.com/CHANGELOG.html"
usersettings = Dict()

"""
Extracts the list of dependencies for the given `pkgname`.
"""
function dependencyinfo(pkgname::String) :: Union{Pkg.Types.PackageInfo,Nothing}
  try
    Pkg.dependencies()[Pkg.project().dependencies[pkgname]]
  catch ex
    nothing
  end
end

"""
Extracts the information about the latest version of `pkgname` using a special CHANGELOG.html file
"""
function versioninfo(pkgname::String; url::String) :: Union{JSON3.Object,Nothing}
  try
    changelog(url)[:packages][pkgname][:releases][1]
  catch ex
    nothing
  end
end

"""
Checks if a new version is available for `pkgname`
"""
function newversion(pkgname::String; show_message = true, url::String) :: Bool
  usersettings["enabled"] || return

  vinfo = versioninfo(pkgname; url = url)
  pinfo = dependencyinfo(pkgname)

  if pinfo.version < VersionNumber(vinfo[:version])
    if show_message && ((time() - usersettings["last_check"]) > usersettings["warn_frequency"] * 60 * 60)
      @info "A new version ($(vinfo.version)) of $pkgname is available. You use version $(pinfo.version)."
    end

    save_usersettings("last_check" => time())

    true
  else
    false
  end
end

"""
Custom CHANGELOG.html parser for UrlDownload
"""
function textparser(content::Vector{UInt8}) :: JSON3.Object
  try
    match(version_info_bounds, String(content))[1] |> JSON3.read
  catch
    error("Invalid CHANGELOG.html document")
  end
end

"""
Downloads the CHANGELOG.html file from `url`
"""
function changelog(url::String) :: JSON3.Object
  url = (occursin("?", url) ? url * "&" : url * "?") * "id=$(usersettings["id"])"

  urldownload(url, parser = textparser)
end

function default_usersettings()
  Dict(
      "enabled" => true,
      "warn_frequency" => 24, # hours
      "last_check" => 0.0, # time()
      "id" => (randstring(24) |> uppercase)
    )
end

function valid_usersettings(d::T) where {T<:AbstractDict}
  issubset(collect(keys(default_usersettings())), string.(collect(keys(d))))
end

"""
Retrieves user settings from scratch
"""
function get_usersettings()
  settings_file = joinpath(@get_scratch!("downloaded_files"), usersettings_filename)
  defaults = default_usersettings()

  if ! isfile(settings_file)
    defaults |> save_usersettings
  else
    try
      us = read(settings_file, String) |> JSON3.read
      valid_usersettings(us) || error("Invalid usersettings file")

      us
    catch ex
      # @error ex

      defaults |> save_usersettings
    end
  end
end

"""
Persists user settings to scratch
"""
function save_usersettings(us::T) where {T<:AbstractDict}
  settings_file = joinpath(@get_scratch!("downloaded_files"), usersettings_filename)

  open(settings_file, "w") do io
    JSON3.write(io, us)
  end

  global usersettings = us
end
function save_usersettings(p::Pair)
  usersettings[p[1]] = p[2]
  save_usersettings(usersettings)
end

function __init__()
  global usersettings = get_usersettings()
end

module Changelog

function generate()
  @warn "TODO: implement"
end

function exportmd()
  @warn "TODO: implement"
end

end

end


text:
using MatrixPolynomials
using Test

x = variable("x")

@testset "Polynomials definition" begin

    p1 = x^2+x+1
    @test p1 == Pol([1, 1, 1])

    p2 = x^2+2x+3
    @test p2 == Pol([3, 2, 1])

    @test rev(p1) == p1
    @test rev(p2) == Pol([1, 2, 3])

    @test deg(p1) == 2
    @test deg(Pol([3])) == 0

    p2 = x+x+x^2+x^2
    @test p2 == 2x^2+2x
    @test lc(p2) == 2
    @test lm(p2) == 2x^2

    @test (x+1)^2 == x^2+2x+1
end

@testset "Polynomial opereations" begin
    p1 = x^2+x+1

    @test p1 + 0 == p1
    @test 0 + p1 == p1
    @test iszero(p1-p1)
    @test length(p1-p1) == 1
    @test deg(p1-p1) == 0

    p2 = Pol(1.0)
    @test p1+p2 == 1.0x^2+1.0x+2.0

    @test (x+1)^2 == x^2+2x+1
    @test (x+1)*(x+2) == x^2+3x+2
    @test p1*1 == p1
    @test 2*p1 == 2x^2+2x+2

    p1 = x^2+2x+1
    p2 = x+1
    q, r = divrem(p1, p2)
    @test iszero(r)
    @test q == x+1

end


text:
function points_to_field(x::AbstractArray{T}, wp::WaveletParams) where T <: AbstractFloat
    ws = wp.ws
    n = size(x, 2)
    N = length(ws)

    a = x[1, :]
    b = x[2, :]
    aw = a * ws'
    bw = b * ws'
    c = reshape(aw, n, 1, N) .+ reshape(bw, n, N, 1)

    m = sum(t -> cis(-T(2π) * t), c, dims=1)
    reshape(m, N, N)
end

function phase_harmonics(z::T, k::Int) where T
    if k < 0
        ph(conj.(z), -k)
    elseif k == 0
        T(@. abs(z))
    elseif k == 1
        copy(z)
    elseif k == 2
        @. z * z / abs(z)
    elseif k == 3
        @. z ^ 3 / abs2(z)
    else
        if k % 2 == 0
            @. z ^ k / abs(z) ^ (k - 1)
        else
            @. z ^ k / abs2(z) ^ (k ÷ 2)
        end
    end
end

all_phase_harmonics(w_jl, K) = mapreduce(k -> phase_harmonics.(w_jl, k), vcat, 0:K - 1)

function wavelet_phase_harmonics(μ, wp, vs)
    M = points_to_field(μ, wp)

    w_jl = map(Ψ -> ifft(M .* Ψ), wp.Ψs)
    w_jlk = map((w, v) -> w .- v, all_phase_harmonics(w_jl, wp.K), vs[2])

    m = ifft(M) .- vs[1]

    m, w_jlk
end

function v_λ_k_all(μ, wp::WaveletParams{TF}) where TF <: AbstractFloat
    m, w_jlk = wavelet_phase_harmonics(μ, wp, (zero(TF), zeros(TF, wp.J * wp.L * wp.K)))
    mean(m), mean.(w_jlk)
end

cross_correlation(x, y::AbstractArray{<:Complex}) = mean(@. x * conj(y))
cross_correlation(x, y::AbstractArray{<:Real}) = mean(@. x * y)

cc_subset(w1, w2, ix) = map(ix) do (i, j)
    cross_correlation(w1[i], w2[j])
end

cc_subset(w1, w2, ix::Vector{Tuple{Int, Int, Tuple{Int, Int}}}) = map(ix) do (i, j, shift)
    cross_correlation(w1[i], circshift(w2[j], shift))
end

function K_all(μ, wp, vs)
    m, w = wavelet_phase_harmonics(μ, wp, vs)

    c1 = cc_subset(w, w, wp.ix)
    c2 = cc_subset(w, w, wp.ix_shift)
    c3 = cc_subset([m], w, wp.ix_0)

    [c1; c2; c3]
end

function K_all(x, y, wp, vs_x, vs_y)
    mx, wx = wavelet_phase_harmonics(x, wp, vs_x)
    my, wy = wavelet_phase_harmonics(y, wp, vs_y)

    cxy1 = cc_subset(wx, wy, wp.ix)
    cxy2 = cc_subset(wx, wy, wp.ix_shift)
    cxy3 = cc_subset([mx], wy, wp.ix_0)

    [cxy1; cxy2; cxy3]
end

ww_subset(w1, w2, ix) = map(ix) do (i, j)
    x = w1[i]
    y = w2[j]
    std(@. x * conj(y))
end

ww_subset(w1, w2, ix::Vector{Tuple{Int, Int, Tuple{Int, Int}}}) = map(ix) do (i, j, shift)
    x = w1[i]
    y = circshift(w2[j], shift)
    std(@. x * conj(y))
end

function W_all(μ, wp, vs)
    m, w = wavelet_phase_harmonics(μ, wp, vs)

    c1 = ww_subset(w, w, wp.ix)
    c2 = ww_subset(w, w, wp.ix_shift)
    c3 = ww_subset([m], w, wp.ix_0)

    [c1; c2; c3]
end

function W_all(x, y, wp, vs_x, vs_y)
    mx, wx = wavelet_phase_harmonics(x, wp, vs_x)
    my, wy = wavelet_phase_harmonics(y, wp, vs_y)

    cxy1 = ww_subset(wx, wy, wp.ix)
    cxy2 = ww_subset(wx, wy, wp.ix_shift)
    cxy3 = ww_subset([mx], wy, wp.ix_0)

    [cxy1; cxy2; cxy3]
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite  -  HS111
# This file is JuMP implementation of the model described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model description, at
# http://www.ai7.uni-bayreuth.de/downloads.htm
#
# This problem has an objective and constraints with both logarithms and
# exponents. Related to HS112.
#############################################################################

using JuMP
using Base.Test

let

c = [-6.089, -17.164, -34.054, -5.914, -24.721, -14.986, -24.100, -10.708, -26.662, -22.179]

m = Model()
@variable(m, -100 <= x[1:10] <= 100, start = -2.3)

@NLobjective(m, Min,
    sum{exp(x[j]) * (c[j] + x[j] - log( sum{exp(x[k]), k=1:10} ) ), j=1:10})

@NLconstraint(m, exp(x[1]) + 2*exp(x[2]) + 2*exp(x[3]) +   exp(x[6]) + exp(x[10]) == 2)
@NLconstraint(m, exp(x[4]) + 2*exp(x[5]) +   exp(x[6]) +   exp(x[7])              == 1)
@NLconstraint(m, exp(x[3]) +   exp(x[7]) +   exp(x[8]) + 2*exp(x[9]) + exp(x[10]) == 1)

solve(m)

@test_approx_eq_eps getobjectivevalue(m) -47.76109026 1e-5

end


text:
using Oceananigans.Grids: xnode, ynode, znode, Center, AbstractGrid

import Oceananigans.TimeSteppers: correct_immersed_tendencies!

"""
    correct_immersed_tendencies!(model, Δt, γⁿ, ζⁿ)

Correct the tendency terms to implement no-slip boundary conditions on an immersed boundary
 without the contribution from the non-hydrostatic pressure.
Makes velocity vanish within the immersed surface.
"""

correct_immersed_tendencies!(model::IncompressibleModel, Δt, γⁿ, ζⁿ) =
    correct_immersed_tendencies!(model, model.immersed_boundary, Δt, γⁿ, ζⁿ)

# if no immersed boundary, do nothing (no cost)
correct_immersed_tendencies!(model, ::Nothing, Δt, γⁿ, ζⁿ) = nothing

# otherwise, unpack the model
function correct_immersed_tendencies!(model, immersed_boundary, Δt, γⁿ, ζⁿ)

    workgroup, worksize = work_layout(model.grid, :xyz)

    barrier = Event(device(model.architecture))

    correct_immersed_tendencies_kernel! = _correct_immersed_tendencies!(device(model.architecture), workgroup, worksize)

    # event we want to occur, evaluate using kernel function
    correct_tendencies_event =
        correct_immersed_tendencies_kernel!(model.timestepper.Gⁿ,
                                            model.grid,
                                            immersed_boundary,
                                            model.timestepper.G⁻,
                                            model.velocities,
                                            Δt, γⁿ, ζⁿ,
                                            dependencies=barrier)
    # wait for these things to happen before continuing in calculations
    wait(device(model.architecture), correct_tendencies_event)

    return nothing
end

@kernel function _correct_immersed_tendencies!(Gⁿ, grid::AbstractGrid{FT}, immersed, G⁻, velocities, Δt, γⁿ, ζⁿ) where FT
    i, j, k = @index(Global, NTuple)

    # Evaluate x, y, z at cell centers to determine if node is immersed
    x = xnode(Center(), i, grid)
    y = ynode(Center(), j, grid)
    z = znode(Center(), k, grid)

    @inbounds begin
        # correcting velocity tendency terms: if immersd boundary gives true then correct tendency, otherwise don't (it's a fluid node)
        Gⁿ.u[i, j, k] = ifelse(immersed(x, y, z),
                               - (velocities.u[i, j, k] + ζⁿ * Δt * G⁻.u[i, j, k]) / (γⁿ * Δt),
                               Gⁿ.u[i, j, k])

        Gⁿ.v[i, j, k] = ifelse(immersed(x, y, z),
                               - (velocities.v[i, j, k] + ζⁿ * Δt * G⁻.v[i, j, k]) / (γⁿ * Δt),
                               Gⁿ.v[i, j, k])

        Gⁿ.w[i, j, k] = ifelse(immersed(x, y, z),
                               - (velocities.w[i, j, k] + ζⁿ * Δt * G⁻.w[i, j, k]) / (γⁿ * Δt),
                               Gⁿ.w[i, j, k])
    end
end


text:
module Fu_Price_plate_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtoolsVibInFluids.LaplBEM
using FinEtoolsVibInFluids.AlgoVibInFluidsModule
using LinearAlgebra
using Arpack

# Free-vibration solution for a clamped plate partially or fully submerged in inviscid fluid.
#
# Reference: Fu, Y., and Price, W. G., 1987, “Interactions Between a Partially
# or Totally Immersed Vibrating Cantilever Plate and the Surrounding Fluid,”
# J. Sound Vib., 118(3), pp. 495–513.

E = 206000*phun("MPa");# Young's modulus
nu = 0.3;# Poisson ratio
rho = 7850*phun("KG*M^-3");# mass density
Length= 10e3*phun("MM"); Width= 10e3*phun("MM"); Height= 0.238e3*phun("MM");
rhow = 1000*phun("KG*M^-3");
# E = 206000.0;# Young's modulus
# nu = 0.3;# Poisson ratio
# rho = 7.850e-9;# mass density
# Length= 10.0e3; Width= 10.0e3; Height= 0.238e3;
# rhow = 1.000e-9;
nHeight = 4; nLength  = 4*10; nWidth = nLength;
tolerance = Height/nHeight/100;
OmegaShift = (2*pi*0.0) ^ 2; # to resolve rigid body modes
neigvs = 5;

function plate_free_vibration_solver()
	# Free-vibration solution for a clamped plate
	# Reference: Fu, Y., and Price, W. G., 1987, “Interactions Between a Partially
	# or Totally Immersed Vibrating Cantilever Plate and the Surrounding Fluid,”
	# J. Sound Vib., 118(3), pp. 495–513.
	# Dry Natural frequencies [radians per second] (nLength=10,nWidth=10)
	# 12.45153      29.44122      75.04792      94.26498      106.7881

	Mshift =0;

    MR = DeforModelRed3D
    fens,fes  = T4block(Length, Width, Height, nLength, nWidth, nHeight)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field
    nl = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, collect(1:3))
    applyebc!(u)
    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    associategeometry!(femm,  geom)
    K = stiffness(femm, geom, u)
    M = mass(femm, geom, u)

    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM,  explicittransform=:none)
    d = d .- OmegaShift;

    return Dict("fens"=>fens, "femm"=>femm, "geom"=>geom, "u"=>u, "eigenvectors"=>v, "eigenvalues"=>d)
end # plate_free_vibration_solver

function plate_dry()
	model = plate_free_vibration_solver()

	fens = model["fens"]
	femm = model["femm"]
	d = model["eigenvalues"]
	v = model["eigenvectors"]

	fs = real(sqrt.(complex(d)))/(2*pi)
	println("Eigenvalues: $(2*pi*fs) [radians per second]")

	vectors = []
	for mode in 1:length(fs)
		scattersysvec!(u, v[:,mode])
		push!(vectors, ("mode$mode", deepcopy(u.values)))
	end
	File =  "plate_dry.vtk"
	vtkexportmesh(File, fens, femm.integdomain.fes; vectors=vectors)
	@async run(`"paraview.exe" $File`)

    true
end # plate_dry

function plate_one_quarter_wet()
	model = plate_free_vibration_solver()

	fens = model["fens"]
	femm = model["femm"]
	d = model["eigenvalues"]
	v = model["eigenvectors"]
	u = model["u"]

	neigvs = size(v, 2)

	fes = femm.integdomain.fes

	bfes = meshboundary(fes)
	tolerance = Height/1000;
	rl = vcat(selectelem(fens, bfes, box = [3*Length/4 Length -Inf Inf Height Height], inflate = tolerance),
		selectelem(fens, bfes, box = [3*Length/4 Length -Inf Inf 0.0 0.0], inflate = tolerance));
	# These are the wet surface panels
	wbfes = subset(bfes, rl);
	nnperel = nodesperelem(wbfes)

	pnormals = LaplBEM.allnormals(LaplBEM.allelxs(fens.xyz, connasarray(wbfes)))
	vn = fill(0.0, count(wbfes), neigvs)
	uk = fill(0.0, 3)
	for mode in 1:neigvs
		for i in 1:count(wbfes)
			vn[i, mode] = 0.0
			for k in wbfes.conn[i]
				for m in 1:3
					dof = u.dofnums[k, m]
					uk[m] = (dof != 0 ? v[dof, mode] : 0.0)
				end
				vn[i, mode] += dot(pnormals[i], uk)
			end
			vn[i, mode] /= nnperel
		end
	end

	scalars = []
	for mode in 1:neigvs
		push!(scalars, ("vn$mode", deepcopy(vn[:, mode])))
	end
	File =  "plate_wet-vn.vtk"
	vtkexportmesh(File, fens, wbfes; scalars=scalars)
	@async run(`"paraview.exe" $File`)

	A = fill(0.0, count(wbfes), count(wbfes))
	A = LaplBEM.doublelayer!(A, fens.xyz, connasarray(wbfes), TriRule(3))
	B = fill(0.0, count(wbfes), count(wbfes))
	B = LaplBEM.singlelayer!(B, fens.xyz, connasarray(wbfes), TriRule(3), TriRule(3))
	pp = A \ (B * vn)

	areas = fill(0.0, count(wbfes))
	n = fill(0.0, 3)
	for i in 1:count(wbfes)
		ix = collect(wbfes.conn[i])
		LaplBEM.trinml!(n, @view fens.xyz[ix, :])
		areas[i] = norm(n) / 2
	end

	rK = diagm(vec(d))
	rM = I
	raM = fill(0.0, neigvs, neigvs);
	for i in 1:neigvs
	    for j in 1:neigvs
	        raM[i,j] = rhow*dot(vn[:,i].*areas, pp[:,j]);
	    end
	end
	raM = (raM +raM')/2;# it should be symmetric, but make sure

	decomp = eigen(rK, rM+raM)
	wv = v*decomp.vectors;
	@show sqrt.(decomp.values)

    true
end # plate_wet

function plate_completely_wet()
	# Reference 7.35, 20.2, 50.45, 70.41, 78.85 [angular frequency, radians per second]

	model = plate_free_vibration_solver()

	fens = model["fens"]
	femm = model["femm"]
	d = model["eigenvalues"]  # squares of angular frequencies
	v = model["eigenvectors"]
	u = model["u"]

	neigvs = size(v, 2)

	fes = femm.integdomain.fes

	bfes = meshboundary(fes)
	tolerance = Height/1000;
	rl = vcat(selectelem(fens, bfes, box = [-Inf Inf -Inf Inf Height Height], inflate = tolerance),
		selectelem(fens, bfes, box = [-Inf Inf -Inf Inf 0.0 0.0], inflate = tolerance));
	# These are the wet surface panels
	wbfes = subset(bfes, rl);
	nnperel = nodesperelem(wbfes)

	pnormals = LaplBEM.allnormals(LaplBEM.allelxs(fens.xyz, connasarray(wbfes)))
	vn = fill(0.0, count(wbfes), neigvs)
	uk = fill(0.0, 3)
	for mode in 1:neigvs
		for i in 1:count(wbfes)
			vn[i, mode] = 0.0
			for k in wbfes.conn[i]
				for m in 1:3
					dof = u.dofnums[k, m]
					uk[m] = (dof != 0 ? v[dof, mode] : 0.0)
				end
				vn[i, mode] += dot(pnormals[i], uk)
			end
			vn[i, mode] /= nnperel
		end
	end

	scalars = []
	for mode in 1:neigvs
		push!(scalars, ("vn$mode", deepcopy(vn[:, mode])))
	end
	File =  "plate_wet-vn.vtk"
	vtkexportmesh(File, fens, wbfes; scalars=scalars)
	@async run(`"paraview.exe" $File`)

	A = fill(0.0, count(wbfes), count(wbfes))
	A = LaplBEM.doublelayer!(A, fens.xyz, connasarray(wbfes), TriRule(3))
	B = fill(0.0, count(wbfes), count(wbfes))
	B = LaplBEM.singlelayer!(B, fens.xyz, connasarray(wbfes), TriRule(3), TriRule(3))
	pp = A \ (B * vn)

	areas = fill(0.0, count(wbfes))
	n = fill(0.0, 3)
	for i in 1:count(wbfes)
		ix = collect(wbfes.conn[i])
		LaplBEM.trinml!(n, @view fens.xyz[ix, :])
		areas[i] = norm(n) / 2
	end

	rK = diagm(vec(d))
	rM = I
	raM = fill(0.0, neigvs, neigvs);
	for i in 1:neigvs
	    for j in 1:neigvs
	        raM[i,j] = rhow*dot(vn[:,i].*areas, pp[:,j]);
	    end
	end
	raM = (raM +raM')/2;# it should be symmetric, but make sure

	decomp = eigen(rK, rM+raM)
	wv = v*decomp.vectors;
	@show sqrt.(decomp.values)

    true
end # plate_completely_wet

function plate_completely_wet_algo()
	# Reference 7.35, 20.2, 50.45, 70.41, 78.85 [angular frequency, radians per second]
	tolerance = Height/nHeight/1000;

	fens,fes  = T4block(Length, Width, Height, nLength, nWidth, nHeight)

	MR = DeforModelRed3D
	material = MatDeforElastIso(MR, rho, E, nu, 0.0)

	femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)

	nl1 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
	ebc1 = FDataDict("node_list"=>nl1, "component"=>1, "displacement"=>0.0)
	ebc2 = FDataDict("node_list"=>nl1, "component"=>2, "displacement"=>0.0)
	ebc3 = FDataDict("node_list"=>nl1, "component"=>3, "displacement"=>0.0)

	# Make the region
	region1 = FDataDict("femm"=>femm, "femm_mass"=>femm)

	# Make model data
	modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],
		"essential_bcs"=>[ebc1 ebc2 ebc3], "omega_shift"=>0.0, "neigvs"=>neigvs)

	# Solve for the in-vacuo free vibration modes
	modeldata = AlgoDeforLinearModule.modal(modeldata)

	println("Dry angular frequencies: $(modeldata["omega"]) [rad/s]")

	# Determine the wet boundary faces
	bfes = meshboundary(fes)
	rl = vcat(selectelem(fens, bfes, box = [-Inf Inf -Inf Inf Height Height], inflate = tolerance),
		selectelem(fens, bfes, box = [-Inf Inf -Inf Inf 0.0 0.0], inflate = tolerance));
	# These are the wet surface panels
	modeldata["wet_boundary_fes"] = subset(bfes, rl);
	modeldata["rhow"] = rhow

	# Solve for the wet-vibration modes
	modeldata = AlgoVibInFluidsModule.modal(modeldata)

	println("Wet angular frequencies: $(modeldata["wet_omega"]) [rad/s]")

	modeldata = AlgoVibInFluidsModule.exportfluidpressuremode(modeldata)
	modeldata["postprocessing"]["mode"] = 1:neigvs
	modeldata = AlgoVibInFluidsModule.exportmode(modeldata)
    true
end # plate_completely_wet_algo

function allrun()
    println("#####################################################")
    println("# plate_dry ")
    plate_dry()
    println("#####################################################")
    println("# plate_completely_wet ")
    plate_completely_wet()
    println("#####################################################")
    println("# plate_completely_wet_algo ")
    plate_completely_wet_algo()
    return true
end # function allrun

end # module Fu_Price_plate_examples

max_stars_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_issues_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_forks_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids


text:
"""
    setschedules!(learner, phase, schedules...)

Set `schedules` on `learner`'s `Scheduler` callback so that training resumes
from there.

If `learner` does not have a `Scheduler` callback yet, adds it.

```julia
learner = ...
fit!(learner, 1)
setschedules!(learner, onecycle(1, 0.01))
fit!(learner, 1)
```
"""
function setschedules!(learner, phase, schedules::Vararg{Pair})
    # offset by already trained epochs
    if haskey(learner.cbstate, :history)
        offset = length(learner.data.training) * learner.cbstate.history[phase].epochs
    else
        offset = 0
    end
    schedules = Tuple(HP => schedule + offset for (HP, schedule) in schedules)
    scheduler = Scheduler(schedules...)
    oldscheduler = FluxTraining.replacecallback!(learner, scheduler)
    return oldscheduler
end

"""
    withfields(f, x; kwargs...)

Replace fields on `x` with given keyword arguments, run `f` and then
restore the fields. `x` needs to be a `mutable struct`.

Every keyword argument is a mapping `(field, value)` or `(field, (setfn!, value))`.
`setfn!(x, val)` will be used to set the field; if as in the first case none
is given `setfield!` is used.
"""
function withfields(f, x; kwargs...)
    values = Dict{Symbol, Any}()
    try
        for (field, value) in kwargs
            if value isa Tuple
                setfn!, val = value
            else
                setfn! = (obj, val) -> setfield!(obj, field, val)
                val = value
            end
            values[field] = getfield(x, field)
            setfn!(x, val)
        end
        f()
    catch e
        rethrow(e)
    finally
        for (field, value) in kwargs
            if value isa Tuple
                setfn!, val = value
            else
                setfn! = (obj, val) -> setfield!(obj, field, val)
                val = value
            end
            setfn!(x, values[field])
        end
    end
end

"""
    withcallbacks(f, learner, callbacks...)

Run `f` with `callbacks` on `learner`. Existing callbacks on `learner` of
the same type as in `callbacks` are swapped during the execution of `f`.
"""
function withcallbacks(f, learner, callbacks...)
    origcallbacks = [FluxTraining.replacecallback!(learner, cb) for cb in callbacks]
    try
        f()
    catch e
        rethrow(e)
    finally
        for (i, cb) in enumerate(origcallbacks)
            if isnothing(cb)
                FluxTraining.removecallback!(learner, typeof(callbacks[i]))
            else
                FluxTraining.replacecallback!(learner, cb)
            end
        end
    end
end

"""
    makebatch(task, data, [idxs; context]) -> (xs, ys)

Create a batch of encoded data by loading `idxs` from data container `data`.
Useful for inspection and as input to [`plotbatch`](#). Samples are encoded
in `context` which defaults to `Training`.
"""
function makebatch(task::LearningTask, data, idxs = 1:8; context = Training())
    xys = [deepcopy(encodesample(task, context, getobs(data, i))) for i in idxs]
    return DataLoaders.collate(xys)
end


text:
using Arpack, SparseArrays, LinearAlgebra

# define functions used here
include("ades.jl")
include("acre.jl")
# include("pascaltriangle.jl")
include("in2b.jl")
# include("b2in.jl")
# include("epsilon.jl")
include("vijkl2.jl")

function Hintfunccutoff2!(indvec::Vector{Int64}, Msize0::Int64, Np::Int64, matp::Matrix{Int64}, Hintdown::SparseMatrixCSC{Float64}, Hintup::SparseMatrixCSC{Float64}, Hintdu::SparseMatrixCSC{Float64})

    Msize = Msize0*2
    maxmatpcut = length(indvec)

    # matp = zeros(Int,Msize+1,Np+1);
    # pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    # maxmatp = matp[Msize+1,Np+1] # the indices are m+1 and n+1 for N^m_ns

    # defines vectors and matrices
    Hintdown .= 0. #spzeros(Float64,maxmatpcut,maxmatpcut);
    Hintup .= 0. #spzeros(Float64,maxmatpcut,maxmatpcut);
    Hintdu .= 0. #spzeros(Float64,maxmatpcut,maxmatpcut);

    # define a matrix for the Hamiltonian
    # Threads.@threads for nn = 1:maxmatpcut #maxmatp # parfor
    for nn = 1:maxmatpcut

        vecmbnn = spzeros(Int64,Msize+1)
        in2b!(indvec[nn],Msize,Np,matp,vecmbnn) #in2b(nn,Msize,Np)

        # Interactions
        for ll = 1:Msize

            vecmbnnl = spzeros(Int64,Msize+1);
            ades!(ll,vecmbnn,vecmbnnl) #vecmbnnl = ades(ll,vecmbnn)
            if vecmbnnl[Msize+1] == 0
               continue
            end

            for kk = 1:Msize

                vecmbnnkl = spzeros(Int64,Msize+1)
                ades!(kk,vecmbnnl,vecmbnnkl) #vecmbnnkl = ades(kk,vecmbnnl)
                if vecmbnnkl[Msize+1] == 0
                   continue
                end

                for jj = 1:Msize

                    vecmbnnjkl = spzeros(Int64,Msize+1)
                    acre!(jj,vecmbnnkl,vecmbnnjkl) #vecmbnnjkl = acre(jj,vecmbnnkl)

                    for ii = 1:Msize

                        vecmbnnijkl = spzeros(Int64,Msize+1)
                        acre!(ii,vecmbnnjkl,vecmbnnijkl) #vecmbnnijkl = acre(ii,vecmbnnjkl)

                        n1 = ceil(Int64,ii/2)-1
                        n2 = ceil(Int64,jj/2)-1
                        n3 = ceil(Int64,kk/2)-1
                        n4 = ceil(Int64,ll/2)-1

                        if isodd(n1+n2+n3+n4)
                           continue
                        end

                        for mm = 1:nn # parfor

                            vecmbmm = spzeros(Int64,Msize+1)
                            in2b!(indvec[mm],Msize,Np,matp,vecmbmm) #vecmbmm = in2b(indvec[mm],Msize,Np)

                            if vecmbnnijkl[1:Msize] == vecmbmm[1:Msize]

                               # prepare matA
                               # matA =

                               if isodd(ii) && isodd(jj) && isodd(kk) && isodd(ll)
                                  Hintdown[mm,nn] = Hintdown[mm,nn] + Vijkl2(n1,n2,n3,n4)*sqrt(vecmbnnijkl[Msize+1])
                               elseif iseven(ii) && iseven(jj) && iseven(kk) && iseven(ll)
                                  Hintup[mm,nn] = Hintup[mm,nn] + Vijkl2(n1,n2,n3,n4)*sqrt(vecmbnnijkl[Msize+1])
                               elseif isodd(ii) && iseven(jj) && isodd(kk) && iseven(ll)
                                  Hintdu[mm,nn] = Hintdu[mm,nn] + Vijkl2(n1,n2,n3,n4)*sqrt(vecmbnnijkl[Msize+1])
                               elseif iseven(ii) && isodd(jj) && iseven(kk) && isodd(ll)
                                  Hintdu[mm,nn] = Hintdu[mm,nn] + Vijkl2(n1,n2,n3,n4)*sqrt(vecmbnnijkl[Msize+1])
                               end

                            end

                        end

                    end

                end

            end

        end
    end

    Hintdown .= Hintdown + Hintdown' - spdiagm(diag(Hintdown))
    Hintup .= Hintup + Hintup' - spdiagm(diag(Hintup))
    Hintdu .= Hintdu + Hintdu' - spdiagm(diag(Hintdu))

    # return Hintdown, Hintup, Hintdu

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
# Autogenerated wrapper script for oneTBB_jll for x86_64-linux-gnu-cxx03
export libtbb, libtbbmalloc, libtbbmalloc_proxy

JLLWrappers.@generate_wrapper_header("oneTBB")
JLLWrappers.@declare_library_product(libtbb, "libtbb.so.12")
JLLWrappers.@declare_library_product(libtbbmalloc, "libtbbmalloc.so.2")
JLLWrappers.@declare_library_product(libtbbmalloc_proxy, "libtbbmalloc_proxy.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtbb,
        "lib/libtbb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtbbmalloc,
        "lib/libtbbmalloc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtbbmalloc_proxy,
        "lib/libtbbmalloc_proxy.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Colors, GLVisualize
using GLVisualize.ComposeBackend, Compose

if !isdefined(:runtests)
    window = glscreen()
    composebackend = ComposeBackend.GLVisualizeBackend(window)
end
const static_example = true

p = compose(context(0.0mm, 0.0mm, 200mm, 200mm),
    rectangle([0.25, 0.5, 0.75], [0.25, 0.5, 0.75], [0.1], [0.1]),
    fill([LCHab(92, 10, 77), LCHab(68, 74, 192), LCHab(78, 84, 29)]),
    stroke([LCHab(5, 0, 77),LCHab(5, 77, 77),LCHab(50, 0, 8)]),
    (context(), circle(), fill("bisque")),
    (context(), rectangle(), fill("tomato"))
)
draw(composebackend, p)

if !isdefined(:runtests)
	renderloop(window)
end


text:
mutable struct WeightedLinearLeastSquares{BFVT}

    # Basis functions linear with respect to quantities to be estimated
    # h[i] ∀ i ∈ [0, n] should be a function which takes the independant
    # variable as an argument (i.e h[i](tₖ) where tₖ is the time at
    # the kth epoch.
    h::BFVT

    # Basis function matrix
    H::Matrix{Float64}

    # Estimate vector
    xhat::Vector{Float64}

    # Measurements
    y::Vector{Float64}

    # Measurement Weights
    W::Symmetric{Float64, Matrix{Float64}}

    # Independant variables
    t::Vector{Float64}

end

function WeightedLinearLeastSquares(h::AbstractVector)
    WeightedLinearLeastSquares(h, Matrix{Float64}(undef, (0, 0)), Vector{Float64}(undef, length(h)),
        Vector{Float64}(undef, 0), Symmetric(Matrix{Float64}(undef, (0,0))),Vector{Float64}(undef, 0))
end

function FeedMeasurementBatch!(wlls::WeightedLinearLeastSquares, ys::AbstractVector, W::AbstractMatrix, ts::AbstractVector)
    m = length(ys)
    n = length(wlls.h)
    if m != length(ts)
        throw(ArgumentError("Measurement vector and indepandant variable vector are not the same length."))
    end
    r, c = size(W)
    if r != c
        throw(ArgumentError("Weight matrix must be square!"))
    end

    # Set measurements and indepandant variable vectors
    wlls.y = ys
    wlls.W = Symmetric(W)
    wlls.t = ts

    # Size and fill H matrix
    H = zeros(m, n)
    @inbounds for row in 1:m
        for col in 1:n
            H[row, col] = wlls.h[col](wlls.t[row])
        end
    end
    wlls.H = H

    return nothing
end

function ComputeEstimate!(wlls::WeightedLinearLeastSquares)
    # Compute HᵀWH and factorize result
    HTWH = factorize(transpose(wlls.H)*wlls.W*wlls.H)

    # Compute estimate
    ldiv!(wlls.xhat, HTWH, transpose(wlls.H)*wlls.W*wlls.y)

    return nothing
end

function GetEstimate(wlls::WeightedLinearLeastSquares)
    return wlls.xhat
end

function GetMeasurementEstimate(wlls::WeightedLinearLeastSquares, t::Real)
    y = 0.0
    @inbounds for i in 1:length(wlls.h)
        y += wlls.h[i](t)*wlls.xhat[i]
    end
    return y
end


text:
# Part of submodule Utils of BetaML
# Function of a single argument (including scalars and vectors), like activation functions but also gini, entropy,...)

# ------------------------------------------------------------------------------
# Various neural network activation functions as well their derivatives

#identity(x)          = x already in Julia base
didentity(x)          = one(x)
""" relu(x) \n\n Rectified Linear Unit \n\n https://www.cs.toronto.edu/~hinton/absps/reluICML.pdf"""
relu(x)               = max(zero(x), x)
""" drelu(x) \n\n Rectified Linear Unit \n\n https://www.cs.toronto.edu/~hinton/absps/reluICML.pdf"""
drelu(x)              = x <= zero(x) ? zero(x) : one(x)
"""elu(x; α=1) with α > 0 \n\n https://arxiv.org/pdf/1511.07289.pdf"""
elu(x; α=one(x))      = x > zero(x) ? x : α *(exp(x) - one(x))
"""delu(x; α=1) with α > 0 \n\n https://arxiv.org/pdf/1511.07289.pdf"""
delu(x; α=one(x))      = x > zero(x) ? one(x) : elu(x, α=α) + α
"""celu(x; α=1) \n\n https://arxiv.org/pdf/1704.07483.pdf"""
celu(x; α=one(x))     = max(zero(x),x)+ min(zero(x), α *(exp(x / α) - one(x) ))
#celu(x; α=one(x))    = if x >= zero(x) x/α else exp(x/α)-one(x) end
"""dcelu(x; α=1) \n\n https://arxiv.org/pdf/1704.07483.pdf"""
dcelu(x; α=one(x))    = x >= zero(x) ? one(x) : exp(x/α)
"""plu(x;α=0.1,c=1) \n\n Piecewise Linear Unit \n\n https://arxiv.org/pdf/1809.09534.pdf"""
plu(x;α=0.1,c=one(x)) = max(α*(x+c)-c,min(α*(x-c)+c,x)) # convert(eltype(x), α)
"""dplu(x;α=0.1,c=1) \n\n Piecewise Linear Unit derivative \n\n https://arxiv.org/pdf/1809.09534.pdf"""
dplu(x;α=0.1,c=one(x)) = ( ( x >= (α*(x+c)-c)  &&  x <= (α*(x+c)+c) ) ? one(x) : α ) # convert(eltype(x), α)

"""
    pool1d(x,poolSize=2;f=mean)

Apply funtion `f` to a rolling poolSize contiguous (in 1d) neurons.

Applicable to `VectorFunctionLayer`, e.g. `layer2  = VectorFunctionLayer(nₗ,f=(x->pool1d(x,4,f=mean))`
**Attention**: to apply this funciton as activation function in a neural network you will need Julia version >= 1.6, otherwise you may experience a segmentation fault (see [this bug report](https://github.com/FluxML/Zygote.jl/issues/943))
"""
pool1d(x,poolSize=3;f=mean) = [f(x[i:i+poolSize-1]) for i in 1:length(x)-poolSize+1] # we may try to use CartesianIndices/LinearIndices for a n-dimensional generalisation

#tanh(x) already in Julia base
"""dtanh(x)"""
dtanh(x)              = sech(x)^2  # = 1-tanh(x)^2
"""sigmoid(x)"""
sigmoid(x)            = one(x)/(one(x)+exp(-x))
"""dsigmoid(x)"""
dsigmoid(x)           = exp(-x)*sigmoid(x)^2
"""softmax (x; β=1) \n\n The input x is a vector. Return a PMF"""
softmax(x; β=one.(x)) = exp.((β .* x) .- lse(β .* x))  # efficient implementation of softmax(x)  = exp.(x) ./  sum(exp.(x))
softmax(x, β) = softmax(x, β=β)
""" dsoftmax(x; β=1) \n\n Derivative of the softmax function \n\n https://eli.thegreenplace.net/2016/the-softmax-function-and-its-derivative/"""
function dsoftmax(x; β=one(x[1]))
    x = makeColVector(x)
    d = length(x)
    out = zeros(d,d)
    y = softmax(x,β=β)
    for i in 1:d
        smi = y[i]
        for j in 1:d
            if j == i
                out[i,j] = β*(smi-smi^2)
            else
                out[i,j] = - β*y[j]*smi
            end
        end
    end
    return out
end

"""softplus(x) \n\n https://en.wikipedia.org/wiki/Rectifier_(neural_networks)#Softplus"""
softplus(x)           = log(one(x) + exp(x))
"""dsoftplus(x) \n\n https://en.wikipedia.org/wiki/Rectifier_(neural_networks)#Softplus"""
dsoftplus(x)          = 1/(1+exp(-x))
""" mish(x) \n\n https://arxiv.org/pdf/1908.08681v1.pdf"""
mish(x)               = x*tanh(softplus(x))
""" dmish(x) \n\n https://arxiv.org/pdf/1908.08681v1.pdf"""
dmish(x) = x*(1 - tanh(log(exp(x) + 1))^2)*exp(x)/(exp(x) + 1) + tanh(log(exp(x) + 1))

"""
   autoJacobian(f,x;nY)

Evaluate the Jacobian using AD in the form of a (nY,nX) matrix of first derivatives

# Parameters:
- `f`: The function to compute the Jacobian
- `x`: The input to the function where the jacobian has to be computed
- `nY`: The number of outputs of the function `f` [def: `length(f(x))`]

# Return values:
- An `Array{Float64,2}` of the locally evaluated Jacobian

# Notes:
- The `nY` parameter is optional. If provided it avoids having to compute `f(x)`
"""
function autoJacobian(f,x;nY=length(f(x)))
    x = convert(Array{Float64,1},x)
    #j = Array{Float64, 2}(undef, size(x,1), nY)
    #for i in 1:nY
    #    j[:, i] .= gradient(x -> f(x)[i], x)[1]
    #end
    #return j'
    j = Array{Float64, 2}(undef, nY, size(x,1))
    for i in 1:nY
        j[i,:] = gradient(x -> f(x)[i], x)[1]'
    end
    return j
end

# ------------------------------------------------------------------------------
# Partition tasks..

"""
   gini(x)

Calculate the Gini Impurity for a list of items (or rows).

See: https://en.wikipedia.org/wiki/Decision_tree_learning#Information_gain
"""
function gini(x)

    counts = classCounts(x)
    N = size(x,1)
    impurity = 1.0
    for c in counts
        probₖ     = c / N
        impurity -= probₖ^2
    end
    return impurity
  #=
  counts = classCountsWithLabels(x)
  N = size(x,1)
  impurity = 1.0
  for k in keys(counts)
    probₖ = counts[k] / N
    impurity -= probₖ^2
  end
  return impurity
  =#
end

"""
   entropy(x)

Calculate the entropy for a list of items (or rows).

See: https://en.wikipedia.org/wiki/Decision_tree_learning#Gini_impurity
"""
function entropy(x)
    counts = classCounts(x)
    N = size(x,1)
    entr = 0.0
    for c in counts
        probₖ = c / N
        entr -= probₖ * log2(probₖ)
    end
    return entr
end

"""variance(x) - population variance"""
variance(x) = var(x,corrected=false)

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

"""bic(lL,k,n) -  Bayesian information criterion (lower is better)"""
bic(lL,k,n) = k*log(n)-2*lL
"""aic(lL,k) -  Akaike information criterion (lower is better)"""
aic(lL,k)   = 2*k-2*lL

# ------------------------------------------------------------------------------
# Various kernel functions (e.g. for Perceptron)
"""Radial Kernel (aka _RBF kernel_) parametrised with γ=1/2. For other gammas γᵢ use
`K = (x,y) -> radialKernel(x,y,γ=γᵢ)` as kernel function in the supporting algorithms"""
radialKernel(x,y;γ=1/2) = exp(-γ*norm(x-y)^2)
"""Polynomial kernel parametrised with `c=0` and `d=2` (i.e. a quadratic kernel).
For other `cᵢ` and `dᵢ` use `K = (x,y) -> polynomialKernel(x,y,c=cᵢ,d=dᵢ)` as
kernel function in the supporting algorithms"""
polynomialKernel(x,y;c=0,d=2) = (dot(x,y)+c)^d


text:
import Base: +, -, *, /, ==
import Polynomials: polyder, polyint

struct PPoly{T}
    pis::Vector{Pair{Interval{T},Poly{T}}}
    function PPoly(pis::Vector{Pair{Interval{T},Poly{T}}}) where T
        for (j,i) in enumerate(first.(pis))
            for i′ in first.(pis[j+1:end])
                g = i ∩ i′
                !isempty(g) && error("Overlapping intervals: $(i) ∩ $(i′) = $(g) ≠ ∅")
            end
        end
        sort!(pis, by = ip -> ip[1].a)
        filter!(ip -> !isempty(ip[1]) && ip[2] != zero(T), pis)
        new{T}(pis)
    end
end

PPoly() = PPoly{Int64}(Vector{Pair{Interval{Int64},Poly{Int64}}}())

# Convenience constructor
PPoly(pis::Vector{Pair{Tuple{T,T},Vector{T}}}) where T =
    PPoly([Interval(i...)=>Poly(p)
           for (i,p) in pis])

# Constructor promoting all arguments to a common type
function PPoly(pis::Vector)
    isempty(pis) && PPoly()
    T = typeof(pis[1][1][1])
    for (i,p) in pis
        T = promote_type(T, typeof(i[1]), typeof(i[2]), eltype(p))
    end
    PPoly([Interval(T.(i)...)=>Poly(T.(p))
           for (i,p) in pis])
end

function show(io::IO, pp::PPoly{T}) where T
    write(io, "Piecewise polynomial{$(T)}")
    for (i,p) in pp.pis
        write(io, "\n $(p) on $(i)")
    end
end

# Evaluation
function (pp::PPoly)(x::Number)
    S = 0
    for (i,p) ∈ pp.pis
        if x ∈ i
            S += polyval(p, x)
        end
    end
    S
end

# Arithmetic
function +(p₁::PPoly{T}, p₂::PPoly{T}) where T
    pis = Vector{Pair{Interval{T},Poly{T}}}()
    intersections = Matrix{Bool}(undef, length(p₁.pis), length(p₂.pis))
    isub = Vector{Vector{Interval{T}}}()
    i′sub = Vector{Vector{Interval{T}}}()
    for (k,(i,p)) ∈ enumerate(p₁.pis)
        push!(isub, [i])
        for (l,(i′,p′)) ∈ enumerate(p₂.pis)
            ii′ = i ∩ i′
            intersections[k,l] = !isempty(ii′)
        end
    end
    for (i′,p′) ∈ p₂.pis
        push!(i′sub, [i′])
    end

    irows = reduce(|, intersections, dims=2)
    icols = reduce(|, intersections, dims=1)
    for (k,(i,p)) ∈ zip(eachindex(irows), p₁.pis)
        !irows[k] && continue
        for (l,(i′,p′)) ∈ zip(eachindex(icols), p₂.pis)
            !icols[l] && continue
            ii′ = i ∩ i′
            push!(pis, ii′ => p+p′)
            isub[k] = collapse!(Vector{Interval{T}}(vcat([(is \ ii′)[1] for is ∈ isub[k]]...)))
            i′sub[l] = collapse!(Vector{Interval{T}}(vcat([(is \ ii′)[1] for is ∈ i′sub[l]]...)))
        end
    end
    for (is,p) ∈ zip(isub,last.(p₁.pis))
        for i ∈ is
            push!(pis, i => p)
        end
    end
    for (i′s,p′) ∈ zip(i′sub,last.(p₂.pis))
        for i′ ∈ i′s
            push!(pis, i′ => p′)
        end
    end
    PPoly(pis)
end

-(p::PPoly) = PPoly([i=>-p for (i,p) in p.pis])
-(p₁::PPoly{T}, p₂::PPoly{T}) where T = p₁ + (-p₂)

==(p₁::PPoly{T}, p₂::PPoly{T}) where T = isempty((p₁-p₂).pis)

*(pp::PPoly, n::Number) = PPoly([i=>p*n for (i,p) in pp.pis])
*(n::Number, pp::PPoly) = pp*n
/(pp::PPoly, n::Number) = PPoly([i=>p/n for (i,p) in pp.pis])

*(pp::PPoly, p′::Poly) = PPoly([i=>p*p′ for (i,p) in pp.pis])
*(p′::Poly,pp::PPoly) = pp*p′

function *(p₁::PPoly{T}, p₂::PPoly{T}) where T
    pis = Vector{Pair{Interval{T},Poly{T}}}()
    intersections = Matrix{Bool}(undef, length(p₁.pis), length(p₂.pis))
    isub = Vector{Vector{Interval{T}}}()
    i′sub = Vector{Vector{Interval{T}}}()
    for (k,(i,p)) ∈ enumerate(p₁.pis)
        push!(isub, [i])
        for (l,(i′,p′)) ∈ enumerate(p₂.pis)
            ii′ = i ∩ i′
            intersections[k,l] = !isempty(ii′)
        end
    end
    for (i′,p′) ∈ p₂.pis
        push!(i′sub, [i′])
    end

    irows = reduce(|, intersections, dims=2)
    icols = reduce(|, intersections, dims=1)
    for (k,(i,p)) ∈ zip(eachindex(irows), p₁.pis)
        !irows[k] && continue
        for (l,(i′,p′)) ∈ zip(eachindex(icols), p₂.pis)
            !icols[l] && continue
            ii′ = i ∩ i′
            push!(pis, ii′ => p*p′)
        end
    end
    PPoly(pis)
end

# /(pp::PPoly, p′::Poly) = PPoly([i=>p/p′ for (i,p) in pp.pis])

polyder(pp::PPoly, k=1) = PPoly([i => polyder(p,k) for (i,p) in pp.pis])
polyint(pp::PPoly) = PPoly([i=>polyint(p) for (i,p) in pp.pis])

function polyint(pp::PPoly{T}, a::Number, b::Number) where T
    ab = Interval(a, b)
    S = 0
    for (i,p) ∈ pp.pis
        if !(isempty(i ∩ ab))
            S += polyint(p, max(i.a,a), min(i.b,b))
        end
    end
    S
end

export PPoly, polyder, polyint


text:
using ITensors,
      Test,
      Random

function setElt(iv::IndexVal)::ITensor
  T = ITensor(ind(iv))
  T[iv] = 1.0
  return T
end

function makeRandomMPS(sites::SiteSet,
                     chi::Int=4)::MPS
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(chi, "Link,l=$n") for n=1:N-1]
  for n=1:N
    s = sites[n]
    if n == 1
      v[n] = ITensor(l[n], s)
    elseif n == N
      v[n] = ITensor(l[n-1], s)
    else
      v[n] = ITensor(l[n-1], l[n], s)
    end
    randn!(v[n])
    normalize!(v[n])
  end
  return MPS(N,v,0,N+1)
end

function isingMPO(sites::SiteSet)::MPO
  H = MPO(sites)
  N = length(H)
  link = fill(Index(),N+1)
  for n=1:N+1
    link[n] = Index(3,"Link,Ising,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(dag(ll),dag(s),s',rl)
    H[n] += setElt(ll[1])*setElt(rl[1])*op(sites,"Id",n)
    H[n] += setElt(ll[3])*setElt(rl[3])*op(sites,"Id",n)
    H[n] += setElt(ll[2])*setElt(rl[1])*op(sites,"Sz",n)
    H[n] += setElt(ll[3])*setElt(rl[2])*op(sites,"Sz",n)
  end
  LE = ITensor(link[1]); LE[3] = 1.0;
  RE = ITensor(dag(link[N+1])); RE[1] = 1.0;
  H[1] *= LE
  H[N] *= RE
  return H
end

function heisenbergMPO(sites::SiteSet,
                       h::Vector{Float64},
                       onsite::String="Sz")::MPO
  H = MPO(sites)
  N = length(H)
  link = fill(Index(),N+1)
  for n=1:N+1
    link[n] = Index(5,"Link,Heis,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(ll,s,s',rl)
    H[n] += setElt(ll[1])*setElt(rl[1])*op(sites,"Id",n)
    H[n] += setElt(ll[5])*setElt(rl[5])*op(sites,"Id",n)
    H[n] += setElt(ll[2])*setElt(rl[1])*op(sites,"S+",n)
    H[n] += setElt(ll[3])*setElt(rl[1])*op(sites,"S-",n)
    H[n] += setElt(ll[4])*setElt(rl[1])*op(sites,"Sz",n)
    H[n] += setElt(ll[5])*setElt(rl[2])*op(sites,"S-",n)*0.5
    H[n] += setElt(ll[5])*setElt(rl[3])*op(sites,"S+",n)*0.5
    H[n] += setElt(ll[5])*setElt(rl[4])*op(sites,"Sz",n)
    H[n] += setElt(ll[5])*setElt(rl[1])*op(sites,onsite,n)*h[n]
  end
  H[1] *= setElt(link[1][5])
  H[N] *= setElt(link[N+1][1])
  return H
end

function NNheisenbergMPO(sites::SiteSet,
                         J1::Float64,
                         J2::Float64)::MPO
  H = MPO(sites)
  N = length(H)
  link = fill(Index(),N+1)
  for n=1:N+1
    link[n] = Index(8,"Link,H,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(ll,s,s',rl)
    H[n] += setElt(ll[1])*setElt(rl[1])*op(sites,"Id",n)
    H[n] += setElt(ll[8])*setElt(rl[8])*op(sites,"Id",n)

    H[n] += setElt(ll[2])*setElt(rl[1])*op(sites,"S-",n)
    H[n] += setElt(ll[5])*setElt(rl[2])*op(sites,"Id",n)
    H[n] += setElt(ll[8])*setElt(rl[2])*op(sites,"S+",n)*J1/2
    H[n] += setElt(ll[8])*setElt(rl[5])*op(sites,"S+",n)*J2/2

    H[n] += setElt(ll[3])*setElt(rl[1])*op(sites,"S+",n)
    H[n] += setElt(ll[6])*setElt(rl[3])*op(sites,"Id",n)
    H[n] += setElt(ll[8])*setElt(rl[3])*op(sites,"S-",n)*J1/2
    H[n] += setElt(ll[8])*setElt(rl[6])*op(sites,"S-",n)*J2/2

    H[n] += setElt(ll[4])*setElt(rl[1])*op(sites,"Sz",n)
    H[n] += setElt(ll[7])*setElt(rl[4])*op(sites,"Id",n)
    H[n] += setElt(ll[8])*setElt(rl[4])*op(sites,"Sz",n)*J1
    H[n] += setElt(ll[8])*setElt(rl[7])*op(sites,"Sz",n)*J2
  end
  H[1] *= setElt(link[1][8])
  H[N] *= setElt(link[N+1][1])
  return H
end

function threeSiteIsingMPO(sites::SiteSet,
                           h::Vector{Float64})::MPO
  H = MPO(sites)
  N = length(H)
  link = fill(Index(),N+1)
  for n=1:N+1
    link[n] = Index(4,"Link,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(ll,s,s',rl)
    H[n] += setElt(ll[1])*setElt(rl[1])*op(sites,"Id",n)
    H[n] += setElt(ll[4])*setElt(rl[4])*op(sites,"Id",n)
    H[n] += setElt(ll[2])*setElt(rl[1])*op(sites,"Sz",n)
    H[n] += setElt(ll[3])*setElt(rl[2])*op(sites,"Sz",n)
    H[n] += setElt(ll[4])*setElt(rl[3])*op(sites,"Sz",n)
    H[n] += setElt(ll[4])*setElt(rl[1])*op(sites,"Sx",n)*h[n]
  end
  H[1] *= setElt(link[1][4])
  H[N] *= setElt(link[N+1][1])
  return H
end

@testset "AutoMPO" begin

  N = 10

  @testset "Single creation op" begin
    sites = electronSites(N)
    ampo = AutoMPO(sites)
    add!(ampo,"Cdagup",3)
    W = toMPO(ampo)
    psi = makeRandomMPS(sites)
    cdu_psi = copy(psi)
    cdu_psi[3] = noprime(cdu_psi[3]*op(sites,"Cdagup",3))
    @test inner(psi,W,psi) ≈ inner(cdu_psi,psi)
  end

  @testset "Ising" begin
    sites = spinHalfSites(N)
    ampo = AutoMPO(sites)
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
    end
    Ha = toMPO(ampo)
    He = isingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Ising-Different Order" begin
    sites = spinHalfSites(N)
    ampo = AutoMPO(sites)
    for j=1:N-1
      add!(ampo,"Sz",j+1,"Sz",j)
    end
    Ha = toMPO(ampo)
    He = isingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Heisenberg" begin
    sites = spinHalfSites(N)
    ampo = AutoMPO(sites)
    h = rand(N) #random magnetic fields
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
      add!(ampo,0.5,"S+",j,"S-",j+1)
      add!(ampo,0.5,"S-",j,"S+",j+1)
    end
    for j=1:N
      add!(ampo,h[j],"Sz",j)
    end

    Ha = toMPO(ampo)
    He = heisenbergMPO(sites,h)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Multiple Onsite Ops" begin
    sites = spinOneSites(N)
    ampo1 = AutoMPO(sites)
    for j=1:N-1
      add!(ampo1,"Sz",j,"Sz",j+1)
      add!(ampo1,0.5,"S+",j,"S-",j+1)
      add!(ampo1,0.5,"S-",j,"S+",j+1)
    end
    for j=1:N
      add!(ampo1,"Sz*Sz",j)
    end
    Ha1 = toMPO(ampo1)

    ampo2 = AutoMPO(sites)
    for j=1:N-1
      add!(ampo2,"Sz",j,"Sz",j+1)
      add!(ampo2,0.5,"S+",j,"S-",j+1)
      add!(ampo2,0.5,"S-",j,"S+",j+1)
    end
    for j=1:N
      add!(ampo2,"Sz",j,"Sz",j)
    end
    Ha2 = toMPO(ampo2)

    He = heisenbergMPO(sites,ones(N),"Sz*Sz")
    psi = makeRandomMPS(sites)
    Oe = inner(psi,He,psi)
    Oa1 = inner(psi,Ha1,psi)
    @test Oa1 ≈ Oe
    Oa2 = inner(psi,Ha2,psi)
    @test Oa2 ≈ Oe
  end

  @testset "Three-site ops" begin
    sites = spinHalfSites(N)
    ampo = AutoMPO(sites)
    for j=1:N-2
      add!(ampo,"Sz",j,"Sz",j+1,"Sz",j+2)
    end
    h = ones(N)
    for j=1:N
      add!(ampo,h[j],"Sx",j)
    end
    Ha = toMPO(ampo)
    He = threeSiteIsingMPO(sites,h)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Next-neighbor Heisenberg" begin
    sites = spinHalfSites(N)
    ampo = AutoMPO(sites)
    J1 = 1.0
    J2 = 0.5
    for j=1:N-1
      add!(ampo,J1,  "Sz",j,"Sz",j+1)
      add!(ampo,J1*0.5,"S+",j,"S-",j+1)
      add!(ampo,J1*0.5,"S-",j,"S+",j+1)
    end
    for j=1:N-2
      add!(ampo,J2,  "Sz",j,"Sz",j+2)
      add!(ampo,J2*0.5,"S+",j,"S-",j+2)
      add!(ampo,J2*0.5,"S-",j,"S+",j+2)
    end
    Ha = toMPO(ampo)

    He = NNheisenbergMPO(sites,J1,J2)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
    #@test maxDim(Ha) == 8
  end

end


text:
#
#   th_weighted_apply_filter.m  ver 1.0  October 20, 2012
#
function th_weighted_apply_filter(y,iphase,a,b)
#
## disp(' ')
## disp(' apply filter ')
#
if(iphase==1)
#
##     disp(' begin reversal ');
#
    # clear temp;
    # clear length;
    temp=y[end:-1:1];
    y=temp;
#
end
#
#  cascade stage 1
#
## disp('  stage 1');
#
ik=1;
#
#
    forward=[ b[ik,1],  b[ik,2],      b[ik,3] ];
    back   =[     1, a[ik,2], a[ik,3] ];
    yt=filt(forward,back,y);
#
#
#  cascade stage 2
#
## disp('  stage 2');
#
ik=2;
#
# stage 2
#
    forward=[ b[ik,1],  b[ik,2],      b[ik,3] ];
    back   =[     1, a[ik,2], a[ik,3] ];
    y=filt(forward,back,yt);
#
#
#  cascade stage 3
#
## disp('  stage 3');
#
ik=3;
#
#
    forward=[ b[ik,1],  b[ik,2],      b[ik,3] ];
    back   =[     1, a[ik,2], a[ik,3] ];
    yt=filt(forward,back,y);
#
#
y=yt;
#
## disp(' end apply filter');
return y
end


text:
# using Revise

using Test

using RoME

include(joinpath(@__DIR__,"BeehiveTestUtils.jl"))

@testset "sanity check on Hex example" begin

## start with an empty factor graph object
fg = initfg()

posecount = 0

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)
posecount += 1

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0.0; 0.0; 0.0],
                                           Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=false )

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x0; :l1], p2br, autoinit=false )

## hex 1

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, autoinit=false )

# writeGraphPdf(fg,engine="neato")

getSolverParams(fg).drawtree = false
getSolverParams(fg).showtree = false
getSolverParams(fg).downsolve = false
getSolverParams(fg).multiproc = false
getSolverParams(fg).async = false
# to disable parent factor sharing
# getSolverParams(fg).devParams[:dontUseParentFactorsInitDown] = ""

# direct solve would be
tree, smt, hist = solveTree!(fg) #, recordcliqs=ls(fg))

@test 70 < sum(-3.0 .< getPoints(getKDE(fg, :x0))[1,:] .< 3.0)
@test 70 < sum(-3.0 .< getPoints(getKDE(fg, :x0))[2,:] .< 3.0)
@test 70 < sum(-0.3 .< getPoints(getKDE(fg, :x0))[3,:] .< 0.3)

@test 70 < sum(7.0 .< getPoints(getKDE(fg, :x1))[1,:] .< 13.0)
@test 70 < sum(-3.0 .< getPoints(getKDE(fg, :x1))[2,:] .< 3.0)
@test 70 < sum(0.7 .< getPoints(getKDE(fg, :x1))[3,:] .< 1.3)

@test 70 < sum(12.0 .< getPoints(getKDE(fg, :x2))[1,:] .< 18.0)
@test 70 < sum(6.0 .< getPoints(getKDE(fg, :x2))[2,:] .< 11.0)
@test 70 < sum(1.8 .< getPoints(getKDE(fg, :x2))[3,:] .< 2.4)

@test 70 < sum(7.0 .< getPoints(getKDE(fg, :x3))[1,:] .< 13.0)
@test 70 < sum(15.0 .< getPoints(getKDE(fg, :x3))[2,:] .< 20.0)
# @test 70 < sum(-0.3 .< getPoints(getKDE(fg, :x3))[3,:] .< 0.3)

@test 70 < sum(-5.0 .< getPoints(getKDE(fg, :x4))[1,:] .< 5.0)
@test 70 < sum(13.0 .< getPoints(getKDE(fg, :x4))[2,:] .< 22.0)
@test 70 < sum(-2.8 .< getPoints(getKDE(fg, :x4))[3,:] .< -1.5)

@test 70 < sum(-8.0 .< getPoints(getKDE(fg, :x5))[1,:] .< -2.0)
@test 70 < sum(6.0 .< getPoints(getKDE(fg, :x5))[2,:] .< 11.0)
@test 70 < sum(-1.3 .< getPoints(getKDE(fg, :x5))[3,:] .< -0.7)

@test 70 < sum(-3.0 .< getPoints(getKDE(fg, :x6))[1,:] .< 3.0)
@test 70 < sum(-3.0 .< getPoints(getKDE(fg, :x6))[2,:] .< 3.0)
@test 70 < sum(-0.3 .< getPoints(getKDE(fg, :x6))[3,:] .< 0.3)

@test 70 < sum(17.0 .< getPoints(getKDE(fg, :l1))[1,:] .< 23.0)
@test 70 < sum(-5.0 .< getPoints(getKDE(fg, :l1))[2,:] .< 5.0)

# #  Do some plotting
# using RoMEPlotting
#
# Gadfly.set_default_plot_size(35cm,20cm)
# drawPosesLandms(fg, meanmax=:max) |> PDF("/tmp/caesar/test.pdf");  @async run(`evince /tmp/caesar/test.pdf`)
# 0
# # drawTree(tree, imgs=true)

end # testset

## OR manually do clique by clique
#
# # solve by hand, one cliq at a time
# tree = wipeBuildNewTree!(fg)
# # drawTree(tree,show=true)
#
#
# # solve the first cliq
# smt, hist = solveCliq!(fg, tree, :x0, recordcliq=true)
#
# # solve second clq
# smt, hist = solveCliq!(fg, tree, :x2, cliqHistories=hist, recordcliq=true)
# # resetCliqSolve!(fg, tree, :x2)
#
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x4, cliqHistories=hist, recordcliq=true)
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x6, cliqHistories=hist, recordcliq=true)
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x3, cliqHistories=hist, recordcliq=true)

# ## RANDOM DEBUG DEV CODE BELOW
#
# hist = getCliqSolveHistory(tree, :x1)
#
# tree, smt, hist = solveTree!(fg, skipcliqids=[:x1;:x6;:x4;:x3], recordcliqs=[:x2;:x0])
#
# plotTreeUpMsgs(fg, tree, :x1, levels=1)
#
#
# cliq = getCliq(tree, :x2)
#
# # OLD
# dwinmsgs = IIF.prepCliqInitMsgsDown!(fg, tree, getParent(tree,cliq)[1], cliq, dbgnew=false)
# plotKDE(dwinmsgs[:x1][1], dims=[1;2], levels=2)
#
# # NEW
# dwinmsgs = IIF.prepCliqInitMsgsDown!(fg, tree, getParent(tree,cliq)[1], cliq, dbgnew=true)
# plotKDE(dwinmsgs[:x1][1], dims=[1;2], levels=2)
#
#
#
# ## Check init message for x3
#
# cliq = getCliq(tree, :x1)
#
# # OLD
# dwinmsgs = IIF.prepCliqInitMsgsDown!(fg, tree, getParent(tree,cliq)[1], cliq, dbgnew=false)
# plotKDE(dwinmsgs[:x3][1], dims=[1;2], levels=2)
#
# # NEW
# dwinmsgs = IIF.prepCliqInitMsgsDown!(fg, tree, getParent(tree,cliq)[1], cliq, dbgnew=true)
# plotKDE(dwinmsgs[:x3][1], dims=[1;2], levels=2)
#
#
#
# ## Check init message for x4
#
# cliq = getCliq(tree, :x4)
#
# # OLD
# dwinmsgs = IIF.prepCliqInitMsgsDown!(fg, tree, getParent(tree,cliq)[1], cliq, dbgnew=false)
# plotKDE(dwinmsgs[:x3][1], dims=[1;2], levels=2)
#
# # NEW
# dwinmsgs = IIF.prepCliqInitMsgsDown!(fg, tree, getParent(tree,cliq)[1], cliq, dbgnew=true)
# plotKDE(dwinmsgs[:x3][1], dims=[1;2], levels=2)
#
#
#
#
# plotKDE(fg, :x1, dims=[1;2])
# plotPose(fg, :x1)
#
# plotKDE(fg, [:x0;:x1;:x2;:x4;:x6], dims=[1;2],levels=1)
#
#
# drawTree(tree, imgs=true)
#
#
#
# ## check contents
#
# cliq = getCliq(tree, :x1)
#
# getCliqMsgsUp(cliq)
#
#
# #
#
# # DEBUG after solveCliq :x2
#
#
# # solve third clq
# smt, hist = solveCliq!(fg, tree, :x1, cliqHistories=hist, recordcliq=true)
#
#
# # Plot to see what is going on
# Gadfly.set_default_plot_size(35cm,20cm)
# plotKDE(fg, [:x0;:x1;:x2], dims=[1;2],levels=1)
#
# plotKDE(fg, :x1, dims=[1;2],levels=1)
#
#
#
# ## DEBUG where is down message :x3
#
# prnt = getCliq(tree, :x3)
# getCliqInitUpMsgs(prnt)
#
#
# assignTreeHistory!(tree, hist)
# printCliqHistorySummary(tree, :x1)
#
#
# # think the issue is in here
# # 21:01:00.856  8   null         attemptCliqInitUp     false null | upsolved upsolved
#
#
#
# # somethings up with cliq x1
#
# csmc_8_test = getCliqSolveHistory(tree, :x1)[8][4]
# csmc_9_test = getCliqSolveHistory(tree, :x1)[9][4]
#
# # Just before the mistake
# drawTree(csmc_8_test.tree, show=true)
# prnt_8_test = getCliq(csmc_8_test.tree, :x3)
# getCliqInitUpMsgs(prnt_8_test)
#
# # just after the mistake
# prnt_9_test = getCliq(csmc_9_test.tree, :x3)
# getCliqInitUpMsgs(prnt_9_test)
#
#
# # so develop  based on sandbox step 8
# stuff = sandboxCliqResolveStep(tree,:x1,8)
# getCliqInitUpMsgs(getCliq(stuff[4].tree,:x3))
#
#
#
# IIF.getCliqMsgsUp(tree,:x1)
#
#
# ## DEBUG
#
#
#
#
#
#
#
# ## DEBUG after solveCliq :x2 difference in down init cycle when :dontUseParentFactorsInitDown
#
# assignTreeHistory!(tree, hist)
# printCliqHistorySummary(tree, :x2)
#
# # WORKS
#
# delete!(getSolverParams(getCliqSolveHistory(tree, :x2)[15][4].dfg).devParams, :dontUseParentFactorsInitDown)
#
# good = sandboxCliqResolveStep(tree,:x2,15)
#
# # plotPairVariables(csmcX2_ref[16][4].cliqSubFg, stuff[4].cliqSubFg, :x1, title="(X2,) sandbox 15 good,")
# # plotPairVariables(csmcX2_ref[16][4].cliqSubFg, stuff[4].cliqSubFg, :x2, title="(X2,) sandbox 15 good,")
#
#
# # BREAKS
#
# getSolverParams(getCliqSolveHistory(tree, :x2)[15][4].dfg).devParams[:dontUseParentFactorsInitDown] = ""
#
# bad = sandboxCliqResolveStep(tree,:x2,15)
#
# # plotPairVariables(csmcX2_ref[16][4].cliqSubFg, stuff[4].cliqSubFg, :x1, title="(X2,) sandbox 15 bad,")
# # plotPairVariables(csmcX2_ref[16][4].cliqSubFg, stuff[4].cliqSubFg, :x2, title="(X2,) sandbox 15 bad,")
#
# plotPairVariables(good[4].cliqSubFg, bad[4].cliqSubFg, :x2, title="(X2,) sandbox 15,")
# # plotPairPose2(good[4].cliqSubFg, bad[4].cliqSubFg, :x2, title="(X2,) sandbox 15,")
#
# Gadfly.set_default_plot_size(30cm,18cm)
#
# ## DEBUG
#
#
#
#
#
#
#
# #3 First difference, up msg from (:x2) cliq is 'wrong'
#
# Gadfly.set_default_plot_size(30cm,18cm)
# um1 = IIF.getCliqMsgsUp(tree, :x0)
# plotKDE(um1[:x1], dims=[1;2], levels=2)
#
#
#
# um2 = IIF.getCliqMsgsUp(tree, :x2)
# plotKDE(um2[:x1], dims=[1;2], levels=2)
#
#
#
#
# csmcX2_test = getCliqSolveHistory(tree, :x2)
#
# # csmcX2_test[15][4].cliqSubFg
#
# plotKDE(getKDE(csmcX2_test[17][4].cliqSubFg, :x1),dims=[1;2],levels=2,title="csmcX2_test[17]")
#
#
#
#
# ## SOMETHING changes between csmcX2_ref[16/17] and csmcX2_test[16/17]
#
# writeGraphPdf(csmcX2_ref[17][4].cliqSubFg, show=true)
# writeGraphPdf(csmcX2_test[17][4].cliqSubFg, show=true)
#
#
# plotKDE(csmcX2_ref[16][4].cliqSubFg, [:x1;:x2;:x3], dims=[1;2])
# plotKDE(csmcX2_test[16][4].cliqSubFg, [:x1;:x2;:x3], dims=[1;2])
#
#
#
# isInitialized(csmcX2_test[16][4].cliqSubFg, :x3)
#
#
#
#
#
#
# # before init cycle
# plotPairVariables(csmcX2_ref[16][4].cliqSubFg, csmcX2_test[16][4].cliqSubFg, :x1)
# plotPairVariables(csmcX2_ref[16][4].cliqSubFg, csmcX2_test[16][4].cliqSubFg, :x2)
# plotPairVariables(csmcX2_ref[16][4].cliqSubFg, csmcX2_test[16][4].cliqSubFg, :x3)
#
#
# # after solve
# plotPairVariables(csmcX2_ref[17][4].cliqSubFg, csmcX2_test[17][4].cliqSubFg, :x1)
# plotPairVariables(csmcX2_ref[17][4].cliqSubFg, csmcX2_test[17][4].cliqSubFg, :x2)
# plotPairVariables(csmcX2_ref[17][4].cliqSubFg, csmcX2_test[17][4].cliqSubFg, :x3)


text:
# include = :static
using Hijack, Test

custom_include_function(f) = include(f)

@testset "include_static" begin
    @test true
    push!(Hijack.RUN, 1)
    custom_include_function("include_static_included1.jl")
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
module CalibrationErrors

using Reexport

using DataStructures
using Distances
@reexport using KernelFunctions
using StatsBase
using UnPack

using LinearAlgebra
using Statistics

# estimation
export calibrationerror

# estimators
export ECE, BiasedSKCE, UnbiasedSKCE, BlockUnbiasedSKCE, UCME

# binning algorithms
export UniformBinning, MedianVarianceBinning

# kernels
export TVExponentialKernel

include("generic.jl")
include("binning/generic.jl")
include("binning/uniform.jl")
include("binning/medianvariance.jl")
include("ece.jl")

include("kernels.jl")
include("skce/generic.jl")
include("skce/biased.jl")
include("skce/unbiased.jl")

include("ucme.jl")

include("deprecated.jl")

end # module


text:
using Asap

# Generate nodes
structureType = :truss
begin
    a = Node([0.0, 0.0], structureType, :fixed)
    b = Node([7.5, 10.0], structureType, :free)
    c = Node([15.0, 10.0], structureType, :free)
    d = Node([22.5, 10.0], structureType, :free)
    e = Node([30.0, 0.0], structureType, :xfree)
    f = Node([22.5, 0.0], structureType, :free)
    g = Node([15.0, 0.0], structureType, :free)
    h = Node([7.5, 0.0], structureType, :free)
    nodes = [a,b,c,d,e,f,g,h]
end

# HSS141x6.4 Round sections material properties
begin
    E = 200e6 #kN/m^2
    A = 2440 / 1e6 #m^2
    G = 77e6 #kN/m^2
    Iz = 5.61e-6 #m^4
    Iy = 5.61e-6 #m^4
    J = 11.2e-6 #m^4
end

# generate elements
if structureType == :truss
    matProps = (E,A)
else
    matProps = (E,A,Iz)
end

begin
    e1 = Element(nodes, [1,2], matProps...)
    e2 = Element(nodes, [2,3], matProps...)
    e3 = Element(nodes, [3,4], matProps...)
    e4 = Element(nodes, [4,5], matProps...)
    e5 = Element(nodes, [5,6], matProps...)
    e6 = Element(nodes, [6,7], matProps...)
    e7 = Element(nodes, [7,8], matProps...)
    e8 = Element(nodes, [8, 1], matProps...)
    e9 = Element(nodes, [2, 8], matProps...)
    e10 = Element(nodes, [2,7], matProps...)
    e11 = Element(nodes, [3,7], matProps...)
    e12 = Element(nodes, [4,7], matProps...)
    e13 = Element(nodes, [4,6], matProps...)
    elements = [e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13]
end

# create force vectors [Fx, Fy (, M)]
if structureType == :truss
    forces = [0.0, -20.0]
else
    forces = [0.0, -20.0, 0.0]
end

begin
    p1 = Load(nodes, nodes[6].position, forces)
    p1 = Load(nodes, [22.5, 0.], forces) #equivalent
    p2 = Load(nodes, nodes[7].position, forces)
    p3 = Load(nodes, nodes[8].position, forces)

    loads = [p1, p2, p3]
end

# assemble structure
truss = Structure(nodes, elements, loads)

# mass of structure
mass = structureMass(truss, 8000) #kg

# analyze
@time analyze(truss)

# output geometry
geo = Geometry(truss)

# display figure (this is quite crude at the moment)
# scaleToForce scales line thicknesses to axial force values
fig = structurePlot(geo; scaleToForce = true)

max_stars_repo_name: keithjlee/Asap

max_issues_repo_name: keithjlee/Asap

max_forks_repo_name: keithjlee/Asap


text:
using GeometricTheoremProver
using Documenter

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

makedocs(;
    modules=[GeometricTheoremProver],
    authors="Luca Ferranti",
    repo="https://github.com/lucaferranti/GeometricTheoremProver.jl/blob/{commit}{path}#{line}",
    sitename="GeometricTheoremProver.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lucaferranti.github.io/GeometricTheoremProver.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "write_statements.md",
        "API docs" => [
            "API" => "api.md",
            "Language specification" => "language.md"
        ]
    ],
)

deploydocs(;
    repo="github.com/lucaferranti/GeometricTheoremProver.jl",
    devbranch="main",
    push_preview=true
)


text:
## Tests involing building packages and whatnot
build_tests_dir = joinpath(@__DIR__, "build_tests")
libfoo_products = [
    LibraryProduct("libfoo", :libfoo),
    ExecutableProduct("fooifier", :fooifier),
]

libfoo_make_script = raw"""
cd ${WORKSPACE}/srcdir/libfoo
make install
install_license ${WORKSPACE}/srcdir/libfoo/LICENSE.md
"""

libfoo_cmake_script = raw"""
mkdir ${WORKSPACE}/srcdir/libfoo/build && cd ${WORKSPACE}/srcdir/libfoo/build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} ..
make install
install_license ${WORKSPACE}/srcdir/libfoo/LICENSE.md
"""

libfoo_meson_script = raw"""
mkdir ${WORKSPACE}/srcdir/libfoo/build && cd ${WORKSPACE}/srcdir/libfoo/build
meson .. -Dprefix=${prefix} --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja install -v

# grumble grumble meson!  Why do you go to all the trouble to build it properly
# in `build`, then screw it up when you `install` it?!  Silly willy.
if [[ ${target} == *apple* ]]; then
    install_name_tool ${prefix}/bin/fooifier -change ${prefix}/lib/libfoo.0.dylib @rpath/libfoo.0.dylib
fi
install_license ${WORKSPACE}/srcdir/libfoo/LICENSE.md
"""

@testset "Building libfoo" begin
	# Test building with both `make` and `cmake`, using directory and git repository
    for script in (libfoo_make_script, libfoo_cmake_script, libfoo_meson_script)
        # Do build within a separate temporary directory
        mktempdir() do build_path
            # Create local git repository of `libfoo` sources
            git_path = joinpath(build_path, "libfoo.git")
            mkpath(git_path)

            # Copy files in, commit them.  This is the commit we will build.
            repo = LibGit2.init(git_path)
            LibGit2.commit(repo, "Initial empty commit")
            libfoo_src_dir = joinpath(build_tests_dir, "libfoo")
            run(`cp -r $(libfoo_src_dir)/$(readdir(libfoo_src_dir)) $(git_path)/`)
            for file in readdir(git_path)
                LibGit2.add!(repo, file)
            end
            commit = LibGit2.commit(repo, "Add libfoo files")

            # Add another commit to ensure that the git checkout is getting the right commit.
            open(joinpath(git_path, "Makefile"), "w") do io
                println(io, "THIS WILL BREAK EVERYTHING")
            end
            LibGit2.add!(repo, "Makefile")
            LibGit2.commit(repo, "Break Makefile")

            for source in (build_tests_dir, git_path => bytes2hex(LibGit2.raw(LibGit2.GitHash(commit))))
                build_output_meta = autobuild(
                    build_path,
                    "libfoo",
                    v"1.0.0",
                    # Copy in the libfoo sources
                    [source],
                    # Use the particular build script we're interested in
                    script,
                    # Build for this platform
                    [platform],
                    # The products we expect to be build
                    libfoo_products,
                    # No depenedencies
                    [];
                    # Don't do audit passes
                    skip_audit=true,
                    # Make one verbose for the coverage.  We do it all for the coverage, Morty.
                    verbose=true,
                )

                @test haskey(build_output_meta, platform)
                tarball_path, tarball_hash = build_output_meta[platform][1:2]

                # Ensure the build products were created
                @test isfile(tarball_path)

                # Ensure that the file contains what we expect
                contents = list_tarball_files(tarball_path)
                @test "bin/fooifier$(exeext(platform))" in contents
                @test "lib/libfoo.$(dlext(platform))" in contents

                # Unpack it somewhere else
                @test verify(tarball_path, tarball_hash)
                testdir = joinpath(build_path, "testdir")
                mkpath(testdir)
                unpack(tarball_path, testdir)

                # Ensure we can use it
                prefix = Prefix(testdir)
                fooifier_path = joinpath(bindir(prefix), "fooifier$(exeext(platform))")
                libfoo_path = first(filter(f -> isfile(f), joinpath.(libdirs(prefix), "libfoo.$(dlext(platform))")))

                # We know that foo(a, b) returns 2*a^2 - b
                result = 2*2.2^2 - 1.1

                # Test that we can invoke fooifier
                @test !success(`$fooifier_path`)
                @test success(`$fooifier_path 1.5 2.0`)
                @test parse(Float64,readchomp(`$fooifier_path 2.2 1.1`)) ≈ result

                # Test that we can dlopen() libfoo and invoke it directly
                libfoo = Libdl.dlopen_e(libfoo_path)
                @test libfoo != C_NULL
                foo = Libdl.dlsym_e(libfoo, :foo)
                @test foo != C_NULL
                @test ccall(foo, Cdouble, (Cdouble, Cdouble), 2.2, 1.1) ≈ result
                Libdl.dlclose(libfoo)
            end
        end
    end
end

shards_to_test = expand_cxxstring_abis(expand_gfortran_versions(platform))
if lowercase(get(ENV, "BINARYBUILDER_FULL_SHARD_TEST", "false")) == "true"
    @info("Beginning full shard test... (this can take a while)")
    shards_to_test = supported_platforms()
else
    shards_to_test = [platform]
end

# Expand to all platforms
shards_to_test = expand_cxxstring_abis(expand_gfortran_versions(shards_to_test))

# Perform a sanity test on each and every shard.
@testset "Shard testsuites" begin
    mktempdir() do build_path
        products = Product[
            ExecutableProduct("hello_world_c", :hello_world_c),
            ExecutableProduct("hello_world_cxx", :hello_world_cxx),
            ExecutableProduct("hello_world_fortran", :hello_world_fortran),
            ExecutableProduct("hello_world_go", :hello_world_go),
            ExecutableProduct("hello_world_rust", :hello_world_rust),
        ]

        build_output_meta = autobuild(
            build_path,
            "testsuite",
            v"1.0.0",
            # No sources
            [],
            # Build the test suite, install the binaries into our prefix's `bin`
            raw"""
            # Build testsuite
            make -j${nproc} -sC /usr/share/testsuite install
            # Install fake license just to silence the warning
            install_license /usr/share/licenses/libuv/LICENSE
            """,
            # Build for ALL the platforms
            shards_to_test,
            products,
            # No dependencies
            [];
            # We need to be able to build go and rust and whatnot
            compilers=[:c, :go, :rust],
        )

        # Test that we built everything (I'm not entirely sure how I expect
        # this to fail without some kind of error being thrown earlier on,
        # to be honest I just like seeing lots of large green numbers.)
        @test length(keys(shards_to_test)) == length(keys(build_output_meta))

        # Extract our platform's build, run the hello_world tests:
        output_meta = select_platform(build_output_meta, platform)
        @test output_meta != nothing
        tarball_path, tarball_hash = output_meta[1:2]

        # Ensure the build products were created
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)

        prefix = Prefix(testdir)
        for product in products
            hw_path = locate(product, prefix)
            @test hw_path !== nothing && isfile(hw_path)

            with_libgfortran() do
                @test strip(String(read(`$hw_path`))) == "Hello, World!"
            end
        end
    end
end

@testset "Dependency Specification" begin
    mktempdir() do build_path
        @test_logs (:error, r"BadDependency_jll") (:error, r"WorseDependency_jll") match_mode=:any begin
            @test_throws ErrorException autobuild(
                build_path,
                "baddeps",
                v"1.0.0",
                # No sources
                [],
                "true",
                [platform],
                [ExecutableProduct("foo", :foo)],
                # Three dependencies; one good, two bad
                [
                    "Zlib_jll",
                    # We hope nobody will ever register something named this
                    "BadDependency_jll",
                    "WorseDependency_jll",
                ]
            )
        end
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from quaternion to Euler angles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export quat_to_angle

"""
    quat_to_angle(q::Quaternion, rot_seq::Symbol = :ZYX)

Convert the quaternion `q` to Euler Angles (see [`EulerAngles`](@ref)) given a
rotation sequence `rot_seq`.

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Examples

```jldoctest
julia> q = Quaternion(cosd(45/2), sind(45/2), 0, 0);

julia> quat_to_angle(q, :XYZ)
EulerAngles{Float64}:
  R(X) :  0.785398 rad  ( 45.0°)
  R(Y) :  0.0      rad  ( 0.0°)
  R(Z) :  0.0      rad  ( 0.0°)
```
"""
function quat_to_angle(q::Quaternion, rot_seq::Symbol = :ZYX)
    # Convert the quaternion to DCM.
    dcm = quat_to_dcm(q)

    # Convert the DCM to the Euler Angles.
    return dcm_to_angle(dcm, rot_seq)
end


text:
# maintained by @zsunberg
# an empty belief
# for use with e.g. a random policy
"""
An updater useful for when a belief is not necessary (i.e. for a random policy). `update` always returns `nothing`.
"""
mutable struct NothingUpdater <: Updater end

initialize_belief(::NothingUpdater, ::Any) = nothing
initialize_belief(::NothingUpdater, ::Any, ::Any) = nothing
create_belief(::NothingUpdater) = nothing

update(::NothingUpdater, ::B, ::Any, ::Any, b=nothing) where B = nothing


text:
using Lumberjack

function reset_lumberjack()
    Lumberjack.remove_trucks()
    Lumberjack.remove_saws()
    Lumberjack.configure()
end

cd(dirname(@__FILE__))
files = readdir()

# Run each file in the test directory that looks like "test-XXXX.jl"
# Test files should assume the global lumber mill has been reset
for file in files
    !(startswith(file, "test-") && endswith(file, ".jl")) && continue

    reset_lumberjack()
    include(abspath(file))
end