text:
# Needed packages
using DataFrames
using CSV
using Distributions
using Plots
using Random
using Images, CoordinateTransformations, TestImages, OffsetArrays, CoordinateTransformations # For data augmentation
using PaddedViews # for cropping images

#=
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 = hcat(train_label,m_train)
scaled_test_inputs = 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

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

# 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

    """
    Leaky ReLU
    f(x)=1(x<0)(αx)+1(x>=0)(x) where α = some constant
    """
    function ReLuL(x)
        if (x > 0.0)
            return x
        elseif (x <= 0.0)
            return x * 0.01
        end
    end

     function ReLUl_derivative(x, alpha=0.01)
          if x <= 0
          end
          return alpha
      end

    # Create the training function
    # Step 1 - feed the data through the 3 layers whilst adjusting the inputs by link weights + squishing with sigmoid function at each output node.
    # Layer 1 = inputs (inputs * weights)
    # Layer 2 = hidden layer = squish with sigmoid the output of layer 1 - the output of hidden layer then is * by link weights (who) which then becomes the input to the output layer
    # Layer 3 = output = squish output of hidden layer by sigmoid to obtain final output
    # Step 2
    # Calculate the error - the network output (what the network thinks the answer is) - the target
    # Step 3 - Update the link weights in proporation to the error (ie if a link weight is higher than another link weight to same node - adjust the larger weight by the bigger fraciton)
    # Output Errors = targets - final_outputs
    # Hidden layer output errors = who' * output_errors
    # Using chain rule - update the link weights from the hidden layer to output (who) and again for the input to the hidden layer (wih)
    # This updating of weights is how the network learns!

    function train(wih::Array{Float64,2}, who::Array{Float64,2},inputs::Array{Float64,2}, targets::Array{Float64,2}; lr::Float64=0.1, drop_out::Bool=true, drop_out_p::Float64=.8) #wih, who, lr, inputs, targets)
        # Matrix multiplication Rules - https://www.mathsisfun.com/algebra/matrix-multiplying.html
        # [m * n] * [n * p] = [m * p] # the last n must match the first n
        # [row, column]
        # wih dim = [200,784]
        # input dim = [784, 1]
        hidden_inputs = (wih * inputs)  # hidd_input dim = [200, 1]

        # Apply sigmoid function to each element the hidden layer output
        # hidden_inputs dim = [200, 1]
        if drop_out == false
        hidden_outputs = sigmoid.(hidden_inputs)
        elseif drop_out == true
        hidden_outputs = sigmoid.(hidden_inputs)
        hidden_outputs .= (hidden_outputs .* rand(Binomial(1,drop_out_p), hnodes)) ./ drop_out_p # 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

        # Apply the weights to the hidden layer outputs which then are the final layers inputs
        # who dims = [10,200]
        # hidden_outputs dims = [200, 1]
        final_inputs = (who * hidden_outputs)

        # Apply sigmoid function to each element the final layer output
        # final_inputs dim = [10,1]
        if drop_out == false
        final_outputs = sigmoid.(final_inputs)
        elseif drop_out == true
        final_outputs = sigmoid.(final_inputs)
        final_outputs .= (final_outputs .* rand(Binomial(1,drop_out_p), onodes)) ./ drop_out_p # 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

        # output layer error is the (target - actual)
        # output_errors dim = [10,1]
        output_errors = targets - final_outputs
        # hidden_errors dim = [200,1]
        # hidden layer error is the output_errors, split by weights, recombined at hidden nodes
        hidden_errors = who' * output_errors

        # Update weights of the neural network - this is how the network learns by adjusting the link weights in propotation to the link weight and error.
        # update the weights for the links between the hidden and output layers
        who .= who .+ (lr .* (output_errors .* final_outputs .* (1.0 .- final_outputs)) * hidden_outputs')
        # update the weights for the links between the input and hidden layers
        wih .= wih .+ (lr .* (hidden_errors .* hidden_outputs .* (1.0 .- hidden_outputs)) * inputs')

        return(wih, who, final_outputs)
    end

    lr = .01
    #=
    # Example of link weight updating (prior to using calclus chain rule)
    # Refine weight on output node (1x node example with two w11 and w21 connected links)
    # Splits the error in proportion to the error ie if w11 is contributing more to error, adjust that by bigger weight.
    # Fraction of e1 to refine w11
    w11 = 6.0  # weight 1, 1 means node connection 1 to node 1
    w21 = 3.0 # weight 2, 1 means node connection 2 to 1
    refine_w11 = w11 / (w11+w21) # 2.3

    # Fraction of e1 to refine w21
    refine_w11 = w21 / (w21+w11) # 1/3

    # Recombine split errors for the links using the error backpropogation
    # ehidden,1 = sum of split errors on links w11 and w12
    # Link weights beween hidden and output layer example
    w11 = 2.0
    w12 = 1.0
    w21 = 3.0
    w22 = 4.0
    e1 = 0.8 # final output error
    e2 = 0.5 # final output error
    # Find the sum of the errors on e1 and e2
    error_out_e1 = e1 * (w11 / (w11 + w21)) + e2 * (w12 / (w12 + w22)) # Sum of split errors
    # Flip to find the other summed link error
    error_out_e2 = e2 * (w22 / (w22 + w12)) + e1 * (w21 / (w21 + w11)) # Sum of split errors

    # Back propogate error_out_e1 and error_out_e2 to the hidden layer input weights
    w11 = 3.0
    w12 = 1.0
    w21 = 2.0
    w22 = 7.0
    e1 = error_out_e1 # Hidden layer output sum of proportional errors
    e2 = error_out_e2 # Hidden layer output sum of proportional errors
    error_out_e1 = e1 * (w11 / (w11 + w21)) + e2 * (w12 / (w12 + w22)) # Sum of split errors
    error_out_e2 = e2 * (w22 / (w22 + w12)) + e1 * (w21 / (w21 + w11))

    # simply in matrix
    (w11 / (w11 + w21)) (w12 / (w12 + w22))     e1
    (w21 / (w21 + w11)) (w22 / (w22 + w12))     e2
    =#

    # Function for applying the network to unseen or test data not seen (trained) by the network before
    function query(wih::Array{Float64,2}, who::Array{Float64,2},inputs::Array{Float64,2})
       # calculate signals into hidden layer
       hidden_inputs = wih * inputs
       # calculate the signals emerging from hidden layer
       hidden_outputs = sigmoid.(hidden_inputs)

       # calculate signals into final output layer
       final_inputs = who * hidden_outputs
       # calculate the signals emerging from final output layer
       final_outputs = sigmoid.(final_inputs)

       return(final_outputs)
    end

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

    # learning rate
    learning_rate = 0.1 # small learning rate for ReLU 0.01

    # Create instance of neural network
    inodes = input_nodes
    hnodes = hidden_nodes
    onodes = output_nodes
    # As randn generates numbers from a standard normal distribution we have to multiply by the standard deviation
    wih = randn(hnodes,inodes) * inodes^(-0.5)
    who = randn(onodes,hnodes) * hnodes^(-0.5)

    # Train the neural network
    # 4.22.2019 - Add functinality to save the numbers which the network did not get correct
    # Epochs is the number of times the training data set is used for training
    epochs = 5
    # Initialize training progress
    learning_curve = fill(0,size(scaled_train_inputs,1))
    training_accuracy = fill(0.0,size(scaled_train_inputs,1))
    # Check that the neural network is converging
    i = 1
    j = 1
    @inbounds for i in 1:epochs
       # go through all records in the training data set
       @inbounds for j in 1:size(scaled_train_inputs,1)
           # Subset by row [j,1] = label [j,2:length(scaled_train_inputs)] is the scaled image data
           label = Int64.(scaled_train_inputs[j,1])
           inputs = reshape(Float64.(scaled_train_inputs[j,2:size(scaled_train_inputs,2)]),784,1) # rotr90() for changing dim to [784,1] from [1,784] to meet matrix multiplication rules
           # create the target output values (all 0.01, except the desired label which is 0.99)
           targets = zeros(10,1)
           targets .= 0.01 .+ targets
           # Set the target .99 to the correct index position within the target array
           targets[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 wih, who, final_outputs = train(wih, who, inputs, targets; lr=learning_rate,drop_out=false,drop_out_p=.8)
           # Check what the network thought the number was
           network_train_output = argmax(final_outputs)[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,size(scaled_train_inputs,1))
           end
           print("This is epoch ",i,"\nThis is iteration ",j,"\nCheck hidden layer output weight links are updating ", who[1],"\n")
           print("Training Accuracy ",(sum(learning_curve) / j),"\n")
           training_accuracy[j] = (sum(learning_curve) / j)
       end
    end

    # Plot the learning curve
    outie = training_accuracy[1:500]
    plot(outie,title="Training Accuracy")

    # Call the network on an indivdual entry
    example_call = Float64.(scaled_train_inputs[500,2:size(scaled_train_inputs,2)])
    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_train_inputs[500,1]
    # Apply the trained network to the input example
    outputs = query(wih, who, example_call)
    # Check what the label thought the number was
    network_label = argmax(outputs)[1]-1
    print("Network thinks the number is ",network_label)
    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

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

    # Input the test set data over the trained neural network to see how well it does
    i=1
    @inbounds for i in 1:size(scaled_test_inputs,1)
        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,2)]),784,1) # rotr90() for changing dim to [784,1] from [1,784] to meet matrix multiplication rules
        # query the network
        outputs = query(wih, who, 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
        # Kaggle out
        test_accuracy[i] = sum(scorecard) / i
        # 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")
        end
    end

ext: jl


text:
using BenchmarkTools, Jets, LinearAlgebra

const SUITE = BenchmarkGroup()

s = JetSpace(Float64, 100, 200)
x = rand(20_000)
SUITE["JetSpace"] = BenchmarkGroup()
SUITE["JetSpace"]["construct"] = @benchmarkable JetSpace(Float64, 2, 3)
SUITE["JetSpace"]["size"] = @benchmarkable size($s)
SUITE["JetSpace"]["eltype"] = @benchmarkable eltype($s)
SUITE["JetSpace"]["ndims"] = @benchmarkable ndims($s)
SUITE["JetSpace"]["rand"] = @benchmarkable rand($s)
SUITE["JetSpace"]["zeros"] = @benchmarkable zeros($s)
SUITE["JetSpace"]["ones"] = @benchmarkable ones($s)
SUITE["JetSpace"]["reshape"] = @benchmarkable reshape($x, $s)
SUITE["JetSpace"]["construct from array"] = @benchmarkable space(x)

dom = JetSpace(Float64, 20)
rng = JetSpace(Float64, 10, 2)
f!(d,m;a) = d .= a .* m.^2
df!(δd,δm;a,mₒ) = δd .= 2 .* a .* mₒ .* δm
a = rand(20)
✈ = Jet(dom=dom, rng=rng, f! = f!, df! = df!, df′! = df!, s = (a = a,))
SUITE["✈"] = BenchmarkGroup()
SUITE["✈"]["construct"] = @benchmarkable Jet(dom=$dom, rng=$rng, f! = $(f!), df! = $(df!), df′! = $(df!), s = (a = $a,))
SUITE["✈"]["domain"] = @benchmarkable domain($✈)
SUITE["✈"]["range"] = @benchmarkable range($✈)
SUITE["✈"]["point"] = @benchmarkable point($✈)
SUITE["✈"]["state"] = @benchmarkable state($✈)
SUITE["✈"]["state!"] = @benchmarkable state!($✈, (a=$(rand(20)),))
SUITE["✈"]["shape"] = @benchmarkable shape($✈)
SUITE["✈"]["size"] = @benchmarkable size($✈)

JopFoo_df!(d,m;diagonal,kwargs...) = d .= diagonal .* m
function JopFoo(diag)
    spc = JetSpace(Float64, length(diag))
    JopLn(;df! = JopFoo_df!, dom = spc, rng = spc, s = (diagonal=diag,))
end
A = JopFoo(rand(100))
m = rand(100)
d = rand(100)
SUITE["JopLn"] = BenchmarkGroup()
SUITE["JopLn"]["construct"] = @benchmarkable JopFoo($(rand(100)))
SUITE["JopLn"]["mul!"] = @benchmarkable mul!($d, $A, $m)
SUITE["JopLn"]["mul!,adjoint"] = @benchmarkable mul!($m, $(A)', $d)
SUITE["JopLn"]["mul"] = @benchmarkable $(A) * $(m)
SUITE["JopLn"]["mul, adjoint"] = @benchmarkable $(A)' * $(d)
SUITE["JopLn"]["adjoint"] = @benchmarkable $(A)'
SUITE["JopLn"]["size"] = @benchmarkable size($A)
SUITE["JopLn"]["shape"] = @benchmarkable shape($A)
SUITE["JopLn"]["domain"] = @benchmarkable domain($A)
SUITE["JopLn"]["range"] = @benchmarkable range($A)

JopBar_f!(d,m;kwargs...) = d .= m.^2
JopBar_df!(δd,δm;mₒ,kwargs...) = δd .= 2 .* mₒ .* δm
function JopBar(n)
    spc = JetSpace(Float64, n)
    JopNl(f! = JopBar_f!, df! = JopBar_df!, dom = spc, rng = spc)
end
F = JopBar(100)
m = rand(100)
d = rand(100)
SUITE["JopNl"] = BenchmarkGroup()
SUITE["JopNl"]["construct"] = @benchmarkable JopBar(100)
SUITE["JopNl"]["mul!"] = @benchmarkable mul!($d, $F, $m)
SUITE["JopNl"]["mul"] = @benchmarkable $F * $m
SUITE["JopNl"]["jacobian"] = @benchmarkable jacobian!($F, $m)
SUITE["JopNl"]["size"] = @benchmarkable size($F)
SUITE["JopNl"]["shape"] = @benchmarkable shape($F)
SUITE["JopNl"]["domain"] = @benchmarkable domain($F)
SUITE["JopNl"]["range"] = @benchmarkable range($F)

G = F ∘ A ∘ F ∘ A
J = jacobian!(G, m)
SUITE["Composition"] = BenchmarkGroup()
SUITE["Composition"]["construct"] = @benchmarkable $F ∘ $A ∘ $F ∘ $A
SUITE["Composition"]["mul!"] = @benchmarkable mul!($d, $G, $m)
SUITE["Composition"]["mul"] = @benchmarkable $G * $m
SUITE["Composition"]["jacobian"] = @benchmarkable jacobian!($G, $m)
SUITE["Composition"]["mul!,adjoint"] = @benchmarkable mul!($m, ($J)', $d)
SUITE["Composition"]["mul, adjoint"] = @benchmarkable ($J)' * $d
SUITE["Composition"]["adjoint"] = @benchmarkable ($J)'
SUITE["Composition"]["size"] = @benchmarkable size($G)
SUITE["Composition"]["shape"] = @benchmarkable shape($G)
SUITE["Composition"]["domain"] = @benchmarkable domain($G)
SUITE["Composition"]["range"] = @benchmarkable range($G)

_F = [JopBar(100) JopBar(100) JopBar(100) ; JopBar(100) JopBar(100) JopBar(100)]
F = @blockop _F
domainF = domain(F)
rangeF = range(F)
m = rand(domain(F))
d = rand(range(F))
e = rand(range(F))
f = rand(range(F))
d′ = rand(Float64,size(range(F)))
e′ = rand(Float64,size(range(F)))
f′ = rand(Float64,size(range(F)))
J = jacobian!(F, m)
SUITE["Block, homogeneous"] = BenchmarkGroup()
SUITE["Block, homogeneous"]["construct"] = @benchmarkable @blockop $_F
SUITE["Block, homogeneous"]["mul!"] = @benchmarkable mul!($d, $F, $m)
SUITE["Block, homogeneous"]["mul"] = @benchmarkable $F * $m
SUITE["Block, homogeneous"]["jacobian"] = @benchmarkable jacobian!($F, $m)
SUITE["Block, homogeneous"]["mul!, adjoint"] = @benchmarkable mul!($m, ($J)', $d)
SUITE["Block, homogeneous"]["mul, adjoint"] = @benchmarkable ($J)' * $d
SUITE["Block, homogeneous"]["adjoint"] = @benchmarkable ($J)'
SUITE["Block, homogeneous"]["shape"] = @benchmarkable shape($F)
SUITE["Block, homogeneous"]["size"] = @benchmarkable size($F)
SUITE["Block, homogeneous"]["domain"] = @benchmarkable domain($F)
SUITE["Block, homogeneous"]["range"] = @benchmarkable range($F)
SUITE["Block, homogeneous"]["block"] = @benchmarkable getblock($m, 2)
SUITE["Block, homogeneous"]["block!"] = @benchmarkable setblock!($m, 2, $(rand(100)))
SUITE["Block, homogeneous"]["broadcast"] = @benchmarkable $f .= $d .+ $e
SUITE["Block, homogeneous"]["broadcast (base-case)"] = @benchmarkable $f′ .= $d′ .+ $e′
SUITE["Block, homogeneous"]["fill!"] = @benchmarkable $f .= $(rand())
SUITE["Block, homogeneous"]["fill! (base-case)"] = @benchmarkable $f′ .= $(rand())
SUITE["Block, homogeneous"]["dot"] = @benchmarkable dot($d, $e)
SUITE["Block, homogeneous"]["dot (base-case)"] = @benchmarkable dot($d′, $e′)
SUITE["Block, homogeneous"]["norm"] = @benchmarkable norm($d)
SUITE["Block, homogeneous"]["norm (base-case)"] = @benchmarkable norm($d′)
SUITE["Block, homogeneous"]["extrema"] = @benchmarkable extrema($d)
SUITE["Block, homogeneous"]["extrema (base-case)"] = @benchmarkable extrema($d′)
SUITE["Block, homogeneous"]["reshape"] = @benchmarkable reshape($d′, $rangeF)

x = rand(100)
_F = [JopBar(100) JopFoo(x) JopBar(100) ; JopBar(100) JopFoo(x) JopBar(100)]
F = @blockop _F
domainF = domain(F)
rangeF = range(F)
m = rand(domain(F))
d = rand(range(F))
e = rand(range(F))
f = rand(range(F))
d′ = rand(Float64,size(range(F)))
e′ = rand(Float64,size(range(F)))
f′ = rand(Float64,size(range(F)))
J = jacobian!(F, m)
SUITE["Block, heterogeneous"] = BenchmarkGroup()
SUITE["Block, heterogeneous"]["construct"] = @benchmarkable @blockop $_F
SUITE["Block, heterogeneous"]["mul!"] = @benchmarkable mul!($d, $F, $m)
SUITE["Block, heterogeneous"]["mul"] = @benchmarkable $F * $m
SUITE["Block, heterogeneous"]["jacobian"] = @benchmarkable jacobian!($F, $m)
SUITE["Block, heterogeneous"]["mul!, adjoint"] = @benchmarkable mul!($m, ($J)', $d)
SUITE["Block, heterogeneous"]["mul, adjoint"] = @benchmarkable ($J)' * $d
SUITE["Block, heterogeneous"]["adjoint"] = @benchmarkable ($J)'
SUITE["Block, heterogeneous"]["shape"] = @benchmarkable shape($F)
SUITE["Block, heterogeneous"]["size"] = @benchmarkable size($F)
SUITE["Block, heterogeneous"]["domain"] = @benchmarkable domain($F)
SUITE["Block, heterogeneous"]["range"] = @benchmarkable range($F)
SUITE["Block, heterogeneous"]["block"] = @benchmarkable getblock($m, 2)
SUITE["Block, heterogeneous"]["block!"] = @benchmarkable setblock!($d, 2, $(rand(100)))
SUITE["Block, heterogeneous"]["broadcast"] = @benchmarkable $f .= $d .+ $e
SUITE["Block, heterogeneous"]["broadcast (base-case)"] = @benchmarkable $f′ .= $d′ .+ $e′
SUITE["Block, heterogeneous"]["fill!"] = @benchmarkable $f .= $(rand())
SUITE["Block, heterogeneous"]["fill! (base-case)"] = @benchmarkable $f′ .= $(rand())
SUITE["Block, heterogeneous"]["reshape"] = @benchmarkable reshape($d′, $rangeF)

JopBaz_df!(d,m;diagonal,kwargs...) = d .= diagonal .* m
function JopBaz(diag)
    spc = JetSpace(Float64, size(diag))
    JopLn(;df! = JopBaz_df!, dom = spc, rng = spc, s = (diagonal=diag,))
end
F = JopBaz(rand(10,10))
rng = range(F)
Fvec = vec(F)
J = jacobian(F, m)
Jvec = vec(J)
m = rand(domain(F))
d = rand(range(F))
mvec = vec(m)
dvec = vec(d)
SUITE["vec"] = BenchmarkGroup()
SUITE["vec"]["construct"] = @benchmarkable vec($F)
SUITE["vec"]["mul!"] = @benchmarkable mul!($dvec, $Fvec, $mvec)
SUITE["vec"]["mul"] = @benchmarkable $Fvec * $mvec
SUITE["vec"]["jacobian"] = @benchmarkable jacobian!($Fvec, $mvec)
SUITE["vec"]["mul!, adjoint"] = @benchmarkable mul!($mvec, ($Jvec)', $dvec)
SUITE["vec"]["mul, adjoint"] = @benchmarkable ($Jvec)' * $dvec
SUITE["vec"]["adjoint"] = @benchmarkable ($Jvec)'
SUITE["vec"]["shape"] = @benchmarkable shape($Fvec)
SUITE["vec"]["size"] = @benchmarkable size($Fvec)
SUITE["vec"]["domain"] = @benchmarkable domain($Fvec)
SUITE["vec"]["range"] = @benchmarkable range($Fvec)
SUITE["vec"]["vec range"] = @benchmarkable vec($rng)

SUITE


text:
"""
    p(contents)

Wrap some text in a paragraph tag. Helps with layout of raw text.
"""
p(contents...) = tag_wrap("p", contents...)

"""
    h1(contents...)

Show contents as a level 1 heading
"""
h1(contents...) = tag_wrap("h1", contents...)

"""
    h2(contents...)

Show contents as a level 2 heading
"""
h2(contents...) = tag_wrap("h2", contents...)

"""
    h3(contents...)

Show contents as a level 3 heading
"""
h3(contents...) = tag_wrap("h3", contents...)

"""
    h4(contents...)

Show contents as a level 4 heading
"""
h4(contents...) = tag_wrap("h4", contents...)

"""
    br()

Insert a line break (html `<br>`) into your UI
"""
br() = UIElement("
")

"""
    dialog(id, title, content, width = 500)

Display a dialog box (which pops over the rest of the content, forcing user acknowledgement).
Visibility is controlled by `id`. A common use case is to define this output function to have
sole input from a button, and simply return whether the button has been clicked.

`title` and `content` can be dynamically set.

If you need to include multiple elements inside `content`, wrap them as a tuple.

!!! note
    For best results `dialog`s should be placed at the _end_ of your UI, after the layout.
    Placing them inside a layout can cause them to be hidden.
"""
function dialog(id, title, content; width = 500)
    (UIElement("""
    <v-dialog
      v-model="$id"
      width="$width">

      <v-card>
        <v-card-title
          primary-title>"""),
    title,
    UIElement("""
    </v-card-title>
        <v-card-text>"""),
    content,
    UIElement("""
        </v-card-text>

        <v-divider></v-divider>

        <v-card-actions>
          <v-spacer></v-spacer>
          <v-btn
            color="primary"
            text
            @click="$id = false"
          >
            Close
          </v-btn>
        </v-card-actions>
      </v-card>
    </v-dialog>
    """),
    UIModel(id, "false"))
end

"""
    snackbar(id, content; color = "error", width = 500)

Display a snackbar at the bottom of the page to notify users of events or other info.
Visibility of the snackbar is controlled by `id`, which should be a boolean. A common use
case is to define this output function to have sole input from a button, and simply return
whether the button has been clicked.

`content` can be dynamically set. Snackbars are small, so you should only include (raw) text
i.e. (`text_output`) in the content, not other `UIElement`s.

If you need to include multiple elements inside `content`, wrap them as a tuple.

Snackbars automatically close based on the timeout set (in milliseconds). Set to 0 to keep
open indefinitely.

!!! note
    For best results `snackbar`s should be placed at the _end_ of your UI, after the layout.
    Placing them inside a layout can cause them to be hidden.
"""
function snackbar(id, content; color = "error", timeout = 6000)
    (UIElement("""
    <v-snackbar
          v-model="$id"
          color = "$color"
          :timeout="$timeout">"""),
    content,
    UIElement("""<v-btn
            text
            @click="$id = false"
          >
            Close
          </v-btn>
        </v-snackbar>"""),
    UIModel(id, "false"))
end

"""
    tooltip(content, tip)

Add a mouseover-activated tooltip to the UI elements contained in `content`. Works best for
small elements like buttons.

If you need to include multiple elements inside `content`, wrap them as a tuple.
"""
function tooltip(content, tip)
    (UIElement("""<v-tooltip bottom>
      <template v-slot:activator="{ on }">
        """),
    content,
    UIElement("""

      </template>
      $tip
    </v-tooltip>
    """))
end

"""
    show_if(id, content...)

Create a span that only shows `content` if the variable `id` is `true`. As with all Matte
logic `id` should be a function defined in the `Server` module of your app.
"""
function show_if(id, content...)
    UIElement(""""""),
    UIModel(id, "false"),
    content,
    UIElement("")
end

"""
    hide_if(id, content...)

Companion to `show_if`. Creates a span that only shows `content` if the variable `id` is
`false`. As with all Matte logic `id` should be a function defined in the `Server` module
of your app.
"""
function hide_if(id, content...)
    UIElement(""""""),
    UIModel(id, "false"),
    content,
    UIElement("")
end

"""
    circular_loader(; width = 3, color = "primary")

Add a circular spinning loading animation to your UI. Best wrapped in a `show_if` so you
can hide it once the relevant content has finished loading.
"""
function circular_loader(; width::Integer = 3, color::AbstractString = "primary")
    UIElement("""<v-progress-circular
      :width="$width"
      color="$color"
      indeterminate
    ></v-progress-circular>""")
end

"""
    divider(; inset = false, vertical = false)

Create a thin line to separate sections of a layout. `inset` adds 72px of indentation to the
left for horizontal dividers, and reduces max-height for vertical dividers.
"""
function divider(; inset = false, vertical = false)
    UIElement("""
    <v-divider
      $(inset ? "inset" : "")
      $(vertical ? "vertical" : "")
    ></v-divider>
    """)
end

"""
    function icon(name)

Add an icon to part of your UI. Supports material design icons. The `name` is
the name of icon - e.g. "home", "heart", etc

For a list of available icons visit [materialdesignicons.com](https://materialdesignicons.com).
You can find their names by clicking on them.

Size can be one of `x-small`, `small`, `medium` (default), `large`, `x-large`.

Works well embedded in buttons.
"""
function icon(name; size = "medium", color = "primary")
    UIElement("""<v-icon $size color = "$color">mdi-$name</v-icon>""")
end


text:
using InfiniteArrays
using Distributions
using DifferentialEquations
using Random
using Plots

# find index of last nonzero element
function find_end(I)
    findfirst(x -> isequal(x, sum(I)), cumsum(I))
end

# find indices of nonzero elements
function find_nonzero(I)
    last = find_end(I)
    findall(>(0), I[1:last])
end

struct SIR_struct
    I::AbstractArray
    R::AbstractArray
end

SIR_struct(I0) = SIR_struct(I0, zeros(Int64, ∞))

function infection_rate(u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    β*c*I/N*S
end

function infection!(integrator, SIR::SIR_struct)

    I_elements = find_nonzero(SIR.I)
    infector_bin = wsample(I_elements, SIR.I[I_elements], 1)[1]

    # infector increases their count of infections by one
    SIR.I[infector_bin] -= 1
    SIR.I[infector_bin + 1] += 1

    # add a 0-infections infector
    SIR.I[1] += 1

    # update S and I
    integrator.u[1] -= 1
    integrator.u[2] = sum(SIR.I)

end

const infection_jump = ConstantRateJump(infection_rate, (integrator) -> infection!(integrator, SIR))

function recovery_rate(u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    γ*I
end

function recovery!(integrator, SIR::SIR_struct)

    I_elements = find_nonzero(SIR.I)
    recovery_bin = wsample(I_elements, SIR.I[I_elements], 1)[1]

    SIR.I[recovery_bin] -= 1
    SIR.R[recovery_bin] += 1

    integrator.u[2] = sum(SIR.I)
    integrator.u[3] = sum(SIR.R)
end

const recovery_jump = ConstantRateJump(recovery_rate, (integrator) -> recovery!(integrator, SIR))

tmax = 40.0
tspan = (0.0,tmax);

δt = 0.1
t = 0:δt:tmax;

u0 = [990,10,0]; # S,I,R

I0 = zeros(Int64, ∞)
I0[1] = u0[2]

SIR = SIR_struct(I0)

p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

prob_discrete = DiscreteProblem(u0,tspan,p);

prob_jump = JumpProblem(prob_discrete,Direct(),infection_jump,recovery_jump);

sol_jump = solve(prob_jump,SSAStepper());

out_jump = sol_jump(t);

plot(
    out_jump,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)

infectors = find_nonzero(SIR.R)
infectors_counts = zeros(Int64, infectors[end])
infectors_counts[infectors] = SIR.R[infectors]

plot(
    infectors_counts ./ sum(infectors_counts),
    seriestype = :bar,
    xlabel="Number",
    ylabel="Proportion",
    color = 1:length(infectors_counts),
    legend = false,
    xticks = 1:length(infectors_counts),
    xformatter = x -> Int(x - 1)
)


text:
using Cameras
using Test

@testset "Acquisition" begin
    image_size = (2, 2)

    function produce_images!(image_source::Channel)
        while true
            put!(image_source, zeros(image_size))
        end
    end

    @testset "Manually triggered" begin
        image_source = Channel(produce_images!)
        camera = SimulatedCamera(image_source)

        @test !isopen(camera)
        open!(camera)
        @test isopen(camera)

        @test !isrunning(camera)
        start!(camera)
        @test isrunning(camera)

        @testset "Synchronous" begin
            trigger!(camera)
            img = take!(camera)
            @test size(img) == image_size
        end

        @testset "Asynchronous" begin
            let img = zeros(0, 0)
                @sync begin
                    @async img = take!(camera)
                    @async trigger!(camera)
                end
                @test size(img) == image_size
            end
        end

        @assert isrunning(camera)
        stop!(camera)
        @test !isrunning(camera)

        @assert isopen(camera)
        close!(camera)
        @test !isopen(camera)
    end

    @testset "Continuously triggered" begin
        period = 0.04 # seconds

        function produce_triggers!(trigger_source::Channel{UInt64})
            while true
                put!(trigger_source, time_ns())
                sleep(period)
            end
        end
        trigger_source = Channel(produce_triggers!; ctype = UInt64)
        image_source = Channel(produce_images!)

        continuous_camera = SimulatedCamera(trigger_source, image_source)

        start!(continuous_camera)

        let i = 0
            while i < 5
                i += 1
                img = take!(continuous_camera)
                @assert size(img) == image_size
            end
            @test i == 5
        end

        @testset "Stopping camera terminates acquisition task" begin
            acquisition_task_started = Condition()
            acquisition_task_terminated = Condition()
            acquisition_task = @task begin
                notify(acquisition_task_started)
                try
                    while true
                        take!(continuous_camera)
                    end
                finally
                    notify(acquisition_task_terminated)
                end
            end
            schedule(acquisition_task)
            wait(acquisition_task_started)
            @assert istaskstarted(acquisition_task)

            @assert isrunning(continuous_camera)
            stop!(continuous_camera)
            @test !isrunning(continuous_camera)
            @test !isopen(trigger_source)
            @test !isopen(image_source)

            wait(acquisition_task_terminated)
            @test istaskdone(acquisition_task)
            @static if VERSION >= v"1.3"
                @test_throws TaskFailedException wait(acquisition_task)
                (e, _) = first(Base.catch_stack(acquisition_task))
                @test e isa InvalidStateException
            else
                @test_throws InvalidStateException wait(acquisition_task)
            end
        end
    end
end


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Durbin-Watson" begin
x_sim = [
    0.2972879845354616, 0.3823959677906078, -0.5976344767282311, -0.01044524463737564,
    -0.839026854388764, 0.31111133849833383, 2.2950878238373105, -2.2670863488005306,
    0.5299655761667461, 0.43142152642291204, 0.5837082875687786, 0.9632716050381906,
    0.45879095505371686, -0.5223367574215084, 0.40839583832475224, -0.050451229933665284,
    -0.6936536438038856, -1.7738332534080556, 0.12076596493743928, -0.7576330377205277,
    -1.7297561813906863, 0.7959486220046159, 0.6700619812560624, 0.5508523775007609,
    -0.06337459242956062, 1.3369427822509223, -0.07314863333773934, -0.7454643166407556,
    -1.2200551285346526, -0.05317733752769253, -0.1651363578028337, -2.115370248427311,
    -0.0667679865065892, 1.2224618484680618, 0.5676954597384435, 0.7634989283233873,
    0.37859887985335006, -0.645597056636597, -0.6646462443910314, -1.8030263242108941,
    -1.2890476038497443, -0.3352700795181175, 0.07046760136115796, 0.34179376737982303,
    1.7351656660543822, 1.299915964310192, 0.20636408140261675, -1.0088599555376034,
    -0.8500556703153979, 1.1294094358843267, -1.348533456491827, 1.2402909004290597,
    -0.055647093207324264, 0.7593902139454598, -0.030266886646762546,
    -0.36108266775609366, -2.0073032025325626, -0.5688591273466379, -1.1496275243969976,
    1.880622174631185, -1.4706245413873802, -0.5352109701657319, -0.9635439598255011,
    -1.3851080676759138, 0.13430023756127, -0.6161166660088192, -1.719989356319742,
    0.32076927998741056, -1.4473728978815112, -0.4922711880695664, -0.4130873839831828,
    1.001928512455953, 1.174681325975111, -1.326970345062047, -1.593214868044067,
    -0.43870997575879667, 1.225797497997867, 1.4890902780696147, 2.1594331871008534,
    -0.9094790180375028, 1.1478203006125018, -0.5050721965223731, -0.2670670352181297,
    1.499082393930729, 0.7973041601753904, -0.17164044330709693, -0.46907960993009634,
    0.21762402120860017, 0.35914640760234967, 0.3200247157080577, 0.2592163280556861,
    0.459695985041628, 0.20998636004214272, 1.5231753911772645, 1.2597946036121963,
    -0.16528310111426064, 0.15619449488763584, -0.8143152943766987, -1.7098682060045913,
    1.4999378761140325, -0.1289156385300165
]

y_p1_sim = [ # simulated with 0.9 error correlation
    -2.201341073669063, -4.307803337512109, -0.09857411277132777, -1.541105730821739,
    -0.027145670802778188, -2.686792369309708, 1.0325029592455204, -4.546430983473314,
    -0.05574145987936152, -1.0484672373273398, -0.0727407422074402, 2.928588488596894,
    2.431548457465139, -1.6824694270670264, -1.1425390408702634, -0.05703441237172502,
    -3.7519965683524736, -2.093296650507312, -3.017127214946707, -5.023733247576096,
    -6.99222421705045, -1.5297654016216695, -7.079699435755043, -3.653535679506873,
    -7.199787077141922, -8.99013125698203, -7.629907186756253, -9.932939819777665,
    -5.792469964809508, -6.07806145541466, -3.499680333201193, -4.202234791648873,
    -5.778858317935674, -4.24373192427373, -2.085017620018707, -3.014335715988717,
    -0.1471613303199173, 1.8578539252085335, -0.23403766619329314, -1.6531056420525947,
    -3.788213810884947, 2.4223665300397057, -0.7782405228814611, -1.446646277164314,
    -0.439089647647926, -0.6568645300510487, -0.1479791111026496, -0.37816853906616543,
    -1.5546979990960859, 0.7934110284198954, 0.5250598397936805, 0.715179599821307,
    -2.6619881622077957, 0.6738306385872681, -0.5829579084083533, -1.6036851727710095,
    1.7676938677132092, 2.1178719441137925, -0.9822453774326113, 2.644398248982846,
    -0.4222365679492326, -0.08802559167948826, 0.004781669065527239, 0.772037726995288,
    -0.589376818822907, -0.009942853176962085, -1.2940732511754347, 3.3672146186572194,
    0.5720957571272467, 1.9994832545202734, -0.2408374457735749, 1.6235559488713027,
    4.765160983845028, 2.2423062077647176, 2.018531795146693, 3.6912668457665974,
    7.052100921101374, 2.891089489994272, 5.012011484448072, 2.9058206661961257,
    4.74433866373019, 4.522254758963782, 4.483117294212917, 7.161931707757217,
    6.608363617973065, 8.54801536915603, 5.975318310443797, 4.404349461089431,
    6.93366550050009, 4.447502046635322, 5.24747472180249, 6.072362283142831,
    5.9553560646539605, 3.2489334527924516, 3.5194287418405414, -0.24943173222751147,
    2.619648178234713, 4.7339671460783475, -1.5079676768571109, 3.805732205983656,
    2.384418539000519
]

# regression
n = length(y_p1_sim)
X = [ones(n, 1) x_sim]
b = X \ y_p1_sim
resid = y_p1_sim - X * b

t = DurbinWatsonTest(X, resid) # asymptotic p-values

@test t.n == 101
@test t.DW ≈ 0.40499929824035813
@test t.p_compute == :ndep
@test pvalue(t) ≈ 5.158964292454297e-16
@test pvalue(t, tail = :left) ≈ 0.9999999999999998
@test pvalue(t, tail = :right) ≈ 2.5794821462271485e-16
@test default_tail(t) == :both
show(IOBuffer(), t)

t = DurbinWatsonTest(X, resid; p_compute = :exact) # exact p-values

@test t.n == 101
@test t.DW ≈ 0.40499929824035813
@test t.p_compute == :exact
@test pvalue(t) ≈ 3.6994415420651786e-24
@test pvalue(t, tail = :left) ≈ 1.00
@test pvalue(t, tail = :right) ≈ 1.8497207710325893e-24
@test default_tail(t) == :both
show(IOBuffer(), t)

y_m1_sim = [ # simulated with -0.9 error correlation
    -1.3910099248613454, -1.5236862020915085, 4.097936666328973, -1.8511224691327444,
    4.089543164538723, -4.495552189456728, 7.111817193652225, -5.978343737983263,
    6.002151715790337, -2.927307499513737, 5.458083999297893, 1.169591175532601,
    4.594296073450898, -4.243581518221936, 2.8299935873046445, -0.5860831006926581,
    -0.9991240800029624, 0.10615464305915895, -1.7044498818807359, -0.9995564116656306,
    -3.5319792844484295, 4.179137229344947, -5.254163800722545, 5.711864125579977,
    -5.990925805935972, 1.5816393282765917, -1.7032484999745658, -0.617286443723021,
    1.5845107369357039, -0.6096898012505031, 3.3496253127300184, -1.1686422345229122,
    -1.5918280931408098, 1.1415907584907719, 3.006908809246799, 0.16279484914286346,
    4.785411299104295, 2.2981183100018683, 0.212429391905895, -1.163256462733893,
    -3.062755396990661, 5.4435487477430335, -3.425834811822814, 2.6321872628475,
    -0.837244387584767, 4.031062468957392, 0.24230351459424804, 2.563723587223389,
    -2.388107267343436, 3.894089606181092, -0.005377707121046793, 1.4234173532692622,
    -1.8229221872128627, 3.8287942020863945, -0.6858661742876399, 1.1178057096772471,
    2.5196873774438067, 0.38418130890537316, -1.4155851844623182, 2.8866125226188566,
    -0.4941340798489638, -0.0011241065963577546, 0.43409066782110806, 0.6425507130120189,
    -1.3986687332385075, 1.8531056183925938, -2.1057647202687244, 5.003685916104994,
    -2.3348521765098624, 3.0528589106154387, -3.234749402470284, 4.142408525714717,
    2.75967828902951, 1.8334770219703578, -0.9220644601902186, 2.498387343549808,
    3.580976074314626, -0.37369360814264785, 4.713559024928862, 0.05056665062182075,
    3.374560008823081, 1.6454107505584676, 1.5983424828203905, 3.9197926529917604,
    2.7024280115597086, 4.662666562844124, -0.805868878530069, 0.9823593012348448,
    2.8591273650074895, -0.12192899074208174, 2.8388011459132807, 1.716213267783003,
    2.7887106550396696, -1.3611491982524841, 3.6582922884633065, -3.3126745262864654,
    5.030080939945388, 1.0915507567022984, -2.9861092295546072, 3.7140067907674967,
    0.6518327501226435
]

# regression
n = length(y_m1_sim)
X = [ones(n, 1) x_sim]
b = X \ y_m1_sim
resid = y_m1_sim - X * b

t = DurbinWatsonTest(X, resid) # asymptotic p-values

@test t.n == 101
@test t.DW ≈ 3.1196941838878423
@test t.p_compute == :ndep
@test pvalue(t) ≈ 1.3691928363065143e-8
@test pvalue(t, tail = :left) ≈ 6.8459641815325715e-9
@test pvalue(t, tail = :right) ≈ 0.9999999931540359
@test default_tail(t) == :both
show(IOBuffer(), t)

t = DurbinWatsonTest(X, resid; p_compute = :exact) # exact p-values

@test t.n == 101
@test t.DW ≈ 3.1196941838878423
@test t.p_compute == :exact
@test pvalue(t) ≈ 6.608698033261362e-10
@test pvalue(t, tail = :left) ≈ 3.304349016630681e-10
@test pvalue(t, tail = :right) ≈ 0.9999999996695651
@test default_tail(t) == :both
show(IOBuffer(), t)

y_0_sim = [ # simulated with 0 error correlation
    -1.7961754992652044, -2.5510957528383353, 3.2525339877180928, 0.19231575061789374,
    1.89169121462802, -1.7386623034795425, 3.258218157382714, -2.526695955245271,
    2.3288443884260106, 0.7381645606308258, 1.8471935105613477, 4.537960965742147,
    2.7213734745790865, -1.9897890454508897, -0.30877316780251873, 1.466080926146517,
    -2.613632233922138, 0.24522161603320347, -1.3710354663088093, -2.4209400297411765,
    -3.4512221745897302, 2.881796133532548, -3.597925434303817, 1.8506552588011762,
    -2.380926529249864, -3.160258392310042, 0.09071892000097082, -2.6081167226985844,
    2.088064405287688, -0.024234312547235515, 2.385739734138283, 0.3967490275831522,
    -2.3202265548315593, 0.33309301826620985, 2.8843408018580714, 0.8655964597465511,
    3.7488337387014004, 4.297643800846097, 0.1873148360133016, -1.207270876248609,
    -3.205052473244388, 4.259413925143798, -0.7425056693856398, -0.598646937182016,
    1.1973080753889698, 1.507929336481593, 2.1567293512995978, 1.2684047056422159,
    -0.6475511763894612, 1.968716146589186, 1.6551464263288553, 0.8306015804336573,
    -1.9237481856587493, 2.628892109084551, 0.7853215622266103, -0.28924845119256015,
    3.3683615196802235, 1.5894237058883518, -1.9790760667912535, 2.570502472637483,
    -0.34918890076289344, -0.07692872949280205, 0.2585418367307264, 0.9004832694436649,
    -1.0522919323231783, 0.5574000211207956, -0.8615471735157796, 3.820189106289126,
    -0.14496612583980942, 1.218044512431157, -1.2637743788791052, 1.5357218567794906,
    4.895903296016806, 1.1354244022005544, 0.36426053387077517, 1.7715587797565924,
    4.779742721710445, -0.30330824012822477, 3.393632860526853, 1.3438900516253283,
    2.7745850292681986, 2.4674323600529258, 1.7461500847342628, 4.242713515247852,
    3.1964332401219315, 4.847669943114067, 0.8363177531165059, -0.35817985387610163,
    3.3565008608192253, 0.32924436697495, 1.9868939670380534, 2.810384666312773,
    2.4117663029348924, -0.4810983070564472, 1.5143233221817027, -1.7185645332767443,
    2.4464053017635212, 3.997453694160127, -3.8861258284250813, 3.094705799661703,
    1.4768492077143098
]

# regression
n = length(y_0_sim)
X = [ones(n, 1) x_sim]
b = X \ y_0_sim
resid = y_0_sim - X * b

t = DurbinWatsonTest(X, resid)

@test t.n == 101
@test t.DW ≈ 2.0315450779948154
@test t.p_compute == :ndep
@test pvalue(t) ≈ 0.8794168553233327
@test pvalue(t, tail = :left) ≈ 0.43970842766166635
@test pvalue(t, tail = :right) ≈ 0.5602915723383337
@test default_tail(t) == :both
show(IOBuffer(), t)

t = DurbinWatsonTest(X, resid; p_compute = :exact) # exact p-values

@test t.n == 101
@test t.DW ≈ 2.0315450779948154
@test t.p_compute == :exact
@test pvalue(t) ≈ 0.8794168553233327 # exact p-values outsid [0,1], use approx.
@test pvalue(t, tail = :left) ≈ 0.43970842766166635
@test pvalue(t, tail = :right) ≈ 0.5602915723383337
@test default_tail(t) == :both
show(IOBuffer(), t)

# test corner cases
t1 = DurbinWatsonTest(X, size(X,1), 0.00, :exact)

@test pvalue(t1) == 0.00
@test pvalue(t1, tail = :left) == 1.00
@test pvalue(t1, tail = :right) == 0.00
show(IOBuffer(), t1)

t1 = DurbinWatsonTest(X, size(X,1), 4.00, :exact)

@test pvalue(t1) == 0.00
@test pvalue(t1, tail = :left) == 0.00
@test pvalue(t1, tail = :right) == 1.00
show(IOBuffer(), t1)
end


text:
"""
Concatenate Stings.
"""
# Write your package code he
function cat_str end

"""
    cat_str(a, b)
Cat two string.
"""
function cat_str(a, b)
    string(a, b)
end

"""
    cat_str(a, b, c)
Cat three string.
"""
function cat_str(a, b, c)
    string(a, b, c)
end


text:
# Autogenerated wrapper script for QDLDL_jll for x86_64-apple-darwin
export libqdldl

JLLWrappers.@generate_wrapper_header("QDLDL")
JLLWrappers.@declare_library_product(libqdldl, "@rpath/libqdldl.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libqdldl,
        "lib/libqdldl.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function dropload(value)
    if value == "make_gaussian_quantiles"
        return gaussian
    elseif value == "make_blobs"
        return blobs
    elseif value == "make_moons"
        return moons
    elseif value == "make_twospirals"
        return spirals
    elseif value == "make_circles"
        return circles
    end
end

function gen_plot(n_samples, cluster_std)
    blobs = SyntheticDatasets.make_blobs(   n_samples = n_samples,
                                            n_features = 2,
                                            centers = [-1 1; -0.5 0.5],
                                            cluster_std = Float64(cluster_std),
                                            center_box = (-2.0, 2.0),
                                            shuffle = true,
                                            random_state = nothing);
    p = @df blobs scatter(:feature_1, :feature_2, group = :label)
    (data = Plots.plotly_series(p), layout = Plots.plotly_layout(p))
end

function gen_gaussian(n_samples, cov,n_classes)
    gaussian = SyntheticDatasets.make_gaussian_quantiles(   n_samples = n_samples,
                                                            cov = Float64(cov),
                                                            n_classes = n_classes,
                                                            random_state = nothing);
    p = @df gaussian scatter(:feature_1, :feature_2, group = :label)
    (data = Plots.plotly_series(p), layout = Plots.plotly_layout(p))
end

function gen_moons(n_samples, noise)
    moons = SyntheticDatasets.make_moons(   n_samples = n_samples,
                                            noise = Float64(noise),
                                            random_state = nothing);
    p = @df moons scatter(:feature_1, :feature_2, group = :label)
    (data = Plots.plotly_series(p), layout = Plots.plotly_layout(p))
end

function gen_two_spirals(n_samples, noise, start_degrees, total_degrees)
    spirals = SyntheticDatasets.make_twospirals(   n_samples = n_samples,
                                                    noise = Float64(noise),
                                                    start_degrees = start_degrees,
                                                    total_degrees = total_degrees);
    p = @df spirals scatter(:feature_1, :feature_2, group = :label)
    (data = Plots.plotly_series(p), layout = Plots.plotly_layout(p))
end

function gen_circles(n_samples, noise, factor)
    circles = SyntheticDatasets.make_circles(   n_samples = n_samples,
                                                noise = Float64(noise),
                                                factor = Float64(factor));
    p = @df circles scatter(:feature_1, :feature_2, group = :label)
    (data = Plots.plotly_series(p), layout = Plots.plotly_layout(p))
end


text:
using Test

@testset "./src/scenarios.jl" begin

using BellScenario

@testset "BlackBox()" begin
    BB = BlackBox(2,3)

    @test BB.num_in == 3
    @test BB.num_out == 2

    @test BB isa BlackBox

    @test_throws DomainError BlackBox(0,2)
    @test_throws DomainError BlackBox(2,0)
    @test_throws DomainError BlackBox(0,-1)
end

@testset "BipartiteNonSignaling()" begin
    scenario = BipartiteNonSignaling(1,2,3,4)
    @test scenario isa Scenario
    @test scenario isa BipartiteNonSignaling

    @test scenario.A == 1
    @test scenario.B == 2
    @test scenario.X == 3
    @test scenario.Y == 4

    @testset "DomainErrors" begin
        @test_throws DomainError BipartiteNonSignaling(0,2,2,2)
        @test_throws DomainError BipartiteNonSignaling(2,0,2,2)
        @test_throws DomainError BipartiteNonSignaling(2,2,0,2)
        @test_throws DomainError BipartiteNonSignaling(2,2,2,0)
    end
end

@testset "BipartiteSignaling()" begin
    @testset "calling w/ BlackBox parameters, no Comm" begin
        BS = BipartiteSignaling((2,3),(4,5))

        @test BS isa BipartiteSignaling
        @test BS isa Scenario

        @test BS.A == BlackBox(2,3)
        @test BS.B == BlackBox(4,5)
        @test BS.dits == 1
        @test BS.bidirectional == false
    end

    @testset "calling w/ BlackBox parameters, no Comm" begin
        BS = BipartiteSignaling((2,2),(2,2), dits =3, bidirectional=true)

        @test BS isa BipartiteSignaling
        @test BS isa Scenario

        @test BS.A == BlackBox(2,2)
        @test BS.B == BlackBox(2,2)
        @test BS.dits == 3
        @test BS.bidirectional == true
    end

    @test_throws DomainError BipartiteSignaling((2,3),(4,5),dits=0)
end

@testset "LocalSignaling()" begin
    @testset "construction" begin
        scenario = LocalSignaling(3,4,2)

        @test scenario isa Scenario
        @test scenario isa LocalSignaling

        @test scenario.X == 3
        @test scenario.Y == 4
        @test scenario.d ==  2
    end

    @testset "DomainError" begin
        @test_throws DomainError LocalSignaling(0,2,2)
    end
end
end


text:
init_grad(x) = zero(x)
zero_grad!(x) = zero(x)
zero_grad!(x::AbstractArray) = (x .= 0)

scan(c::Call) = foreach(scan, c.args)

function scan(x::Tracked)
  x.isleaf && return
  ref = x.ref += 1
  if ref == 1
    scan(x.f)
    isdefined(x, :grad) && (x.grad = zero_grad!(x.grad))
  end
  return
end

function scan(x)
  istracked(x) && scan(tracker(x))
  return
end

function back_(c::Call, Δ, once)
  Δs = c.func(Δ)
  (Δs isa Tuple && length(Δs) >= length(c.args)) ||
    error("Gradient is not a tuple of length $(length(c.args))")
  foreach((x, d) -> back(x, d, once), c.args, data.(Δs))
end

back_(::Call{Nothing}, Δ, once) = nothing
back_(::Call{Missing}, Δ, once) = error("`back!` was already used")

accum!(x, Δ) = x .+ Δ
accum!(x::AbstractArray, Δ) = (x .+= Δ)

function back(x::Tracked, Δ, once)
  x.isleaf && (x.grad = accum!(x.grad, Δ); return)
  ref = x.ref -= 1
  grad = if isdefined(x, :grad)
    x.grad = accum!(x.grad, Δ)
  elseif ref > 0
    x.grad = Δ
  else
    Δ
  end
  if ref == 0
    back_(x.f, grad, once)
    once && !x.isleaf && (x.f = Call(missing, ()))
  end
  return
end

back(::Nothing, Δ, once) = return

# Interface methods

# TODO: if an error occurs in `back` the refcounts will be broken
# and `back` will silently fail to update.
# (but only if you re-use intermediate values between passes)
# Refcounts are also probably not safe in some situations (e.g. back called
# from within a backpropagator)

function back!(x, Δ; once = true)
  istracked(x) || return
  scan(x)
  back(tracker(x), Δ, once)
  return
end

function gradient_(f, xs...)
  xs = param.(data.(xs))
  l = f(xs...)
  losscheck(l)
  back!(l)
  nobacksies("Use `gradient(...; nest = true)` for nested derivatives",
             grad.(xs))
end

# Out-of-place gradients

struct Params
  order::Vector{Any}
  params::IdSet{Any}
  Params() = new([], IdSet())
end

@forward Params.order Base.iterate, Base.length

function Base.push!(ps::Params, x)
  if !(x in ps.params)
    push!(ps.order, x)
    push!(ps.params, x)
  end
  return ps
end

Base.push!(ps::Params, x...) = (foreach(x -> push!(ps, x), x); ps)

Params(xs) = push!(Params(), xs...)

function Base.show(io::IO, ps::Params)
  print(io, "Params([")
  join(io, ps.order, ", ")
  print(io, "])")
end

struct Grads
  grads::IdDict{Any,Any}
end

Base.show(io::IO, ps::Grads) = println(io, "Grads(...)")

Grads() = Grads(IdDict())

@forward Grads.grads Base.setindex!, Base.haskey, Base.length, Base.iterate

Grads(ps::Params) = Grads(IdDict(tracker(p) => init_grad(data(p)) for p in ps))

Base.getindex(g::Grads, x::Tracked) = g.grads[x]

function Base.getindex(g::Grads, x)
  istracked(x) || error("Object not tracked: $x")
  g[tracker(x)]
end

accum!(g::Grads, x, Δ) = g[x] = haskey(g, x) ? g[x] .+ Δ : Δ

function back_(g::Grads, c::Call, Δ)
  Δs = c.func(Δ)
  (Δs isa Tuple && length(Δs) >= length(c.args)) ||
    error("Gradient is not a tuple of length $(length(c.args))")
  foreach((x, Δ) -> back(g, x, Δ), c.args, Δs)
end

back_(g::Grads, ::Call{Nothing}, Δ) = nothing

function back(g::Grads, x::Tracked, Δ)
  x.isleaf && (accum!(g, x, Δ); return)
  ref = x.ref -= 1
  if ref > 0 || haskey(g, x)
    accum!(g, x, Δ)
    ref == 0 && back_(g, x.f, g[x])
  else
    ref == 0 && back_(g, x.f, Δ)
  end
  return
end

back(::Grads, ::Nothing, _) = return

collectmemaybe(xs) = xs

function forward(f, ps::Params)
  y = collectmemaybe(f())
  y, function (Δ)
    g = Grads(ps)
    if istracked(y)
      scan(y)
      back(g, tracker(y), Δ)
    end
    return g
  end
end

function forward(f, args...)
  args = param.(args)
  y, back = forward(() -> f(args...), Params(args))
  y, Δ -> getindex.(Ref(back(Δ)), args)
end

function losscheck(x)
  x isa Real || error("Function output is not scalar")
  isinf(x) && error("Loss is infinite")
  isnan(x) && error("Loss is NaN")
end

function gradient_nested(f, args...)
  y, back = forward(f, args...)
  losscheck(y)
  return back(1)
end

gradient(f, xs...; nest = false) =
  nest ? gradient_nested(f, xs...) : gradient_(f, xs...)

gradient(f, ps::Params) = gradient_nested(f, ps)

# Jacobians and Hessians

import ..Flux

"""
    J = jacobian(m,x)

Calculate the output jacobian `J = d/dx m(x)` such that each row `i` of `J` corresponds to the gradient `J[i,:] = ∇ₓ(m(x)[i])`
"""
function jacobian(m,x)
    xp = param(x)
    y  = m(xp)
    k  = length(y)
    n  = length(x)
    J  = Matrix{eltype(x)}(undef,k,n)
    for i = 1:k
        Flux.back!(y[i], once = false) # Populate gradient accumulator
        J[i,:] = xp.grad
        xp.grad .= 0 # Reset gradient accumulator
    end
    J
end

hessian(f, x) = jacobian(x -> gradient(f, x, nest=true)[1], x)


text:
function myfunc()
    A = rand(200, 200, 400)
    maximum(A)
end

myfunc() # run once to force compilation

using Profile

@profile myfunc()

Profile.print()

Profile.print(format = :flat, sortedby = :overhead)

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


text:
using ADFPEPS
using ADFPEPS:double_ipeps_energy
using CUDA
using Random

CUDA.allowscalar(false)
Random.seed!(100)
model = Hubbard(1.0,12.0,6.0)
folder = "E:/1 - research/4.9 - AutoDiff/data/ADFPEPS/Hubbard/Z2/"
ipeps, key = init_ipeps(model; Ni = 2, Nj = 2, symmetry = :Z2, atype = CuArray, folder = folder, D=2, χ=10, tol=1e-10, maxiter=10)
# folder, model, Ni, Nj, atype, D, χ, tol, maxiter = key
# key = folder, model, Ni, Nj, atype, D, χ, tol, maxiter
# double_ipeps_energy(atype(ipeps), key)
optimiseipeps(ipeps, key; f_tol = 1e-6, opiter = 0, verbose = true)


text:
struct CompressedWordVectors{Q,U,D,T,S,H}
    vocab::Vector{S}                   # vocabulary
    vectors::QuantizedMatrix{Q,U,D,T}  # quantized vectors
    vocab_hash::Dict{S,H}              # word to vector index (column)
end

# Constructors
function CompressedWordVectors(vocab::AbstractVector{S},
                               vectors::QuantizedMatrix{Q,U,D,T}) where {Q,U,D,T,S}
    length(vocab) == size(vectors, 2) ||
        throw(DimensionMismatch("Dimension of vocab and vectors are inconsistent."))
    vocab_hash = Dict{S, Int}()
    for (i, word) in enumerate(vocab)
        vocab_hash[word] = i
    end
    CompressedWordVectors(vocab, vectors, vocab_hash)
end

function Base.show(io::IO, cwv::CompressedWordVectors{Q,U,D,T,S,H}) where {Q,U,D,T,S,H}
    len_vecs, num_words = size(cwv.vectors)
    print(io, "Compressed WordVectors $(num_words) words, $(len_vecs)-element $(T) vectors")
end

"""
    compress(wv [;kwargs...])

Compresses `wv::WordVectors` by using array quantization.

# Keyword arguments
  * `sampling_ratio::AbstractFloat` specifies the percentage of vectors to use
for quantization codebook creation
  * `k::Int` number of quantization values for a codebook
  * `m::Int` number of codebooks to use
  * `method::Symbol` specifies the array quantization method
  * `distance::PreMetric` is the distance

Other keyword arguments specific to the quantization methods can also be provided.
"""
function compress(wv::WordVectors{S,T,H};
                  sampling_ratio::AbstractFloat=DEFAULT_QUANTIZATION_SAMPLING_RATIO,
                  k::Int=DEFAULT_QUANTIZAION_K,
                  m::Int=DEFAULT_QUANTIZATION_M,
                  method::Symbol=DEFAULT_QUANTIZATION_METHOD,
                  distance::Distances.PreMetric=DEFAULT_QUANTIZATION_DISTANCE,
                  kwargs...) where {S,T,H}
    # Checks, initializations#
    @assert 0.0 < sampling_ratio <= 1.0 "The sampling ratio must be in >0 and <=1"
    _, n = size(wv)
    ns = clamp(round(Int, sampling_ratio * n), 1, n)

    # Sample vectors, build quantizer, quantize everything and return
    @debug "Building quantizer using $ns vectors..."
    svecs = wv.vectors[:, sample(1:n, ns, replace=false)]
    aq = build_quantizer(svecs, k=k, m=m, method=method, distance=distance; kwargs...)
    @debug "Quantizing $n vectors..."
    qvecs = quantize(aq, wv.vectors)
    return CompressedWordVectors(wv.vocab, qvecs)
end

"""
    vocabulary(cwv)

Return the vocabulary as a vector of words of the CompressedWordVectors `cwv`.
"""
vocabulary(cwv::CompressedWordVectors) = cwv.vocab

"""
    in_vocabulary(cwv, word)

Return `true` if `word` is part of the vocabulary of the CompressedWordVector `cwv` and
`false` otherwise.
"""
in_vocabulary(cwv::CompressedWordVectors, word::AbstractString) = word in cwv.vocab

"""
    size(cwv)

Return the word vector length and the number of words as a tuple.
"""
size(cwv::CompressedWordVectors) = size(cwv.vectors)

"""
    index(cwv, word)

Return the index of `word` from the CompressedWordVectors `cwv`.
"""
index(cwv::CompressedWordVectors, word) = cwv.vocab_hash[word]

"""
    get_vector(cwv, word)

Return the vector representation of `word` from the CompressedWordVectors `cwv`.
"""
get_vector(cwv::CompressedWordVectors, word) =
      (idx = cwv.vocab_hash[word]; cwv.vectors[:,idx])

"""
    cosine(cwv, word, n=10)

Return the position of `n` (by default `n = 10`) neighbors of `word` and their
cosine similarities.
"""
function cosine(cwv::CompressedWordVectors, word, n=10)
    metrics = cwv.vectors'*get_vector(cwv, word)
    topn_positions = sortperm(metrics[:], rev = true)[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    similarity(cwv, word1, word2)

Return the cosine similarity value between two words `word1` and `word2`.
"""
function similarity(cwv::CompressedWordVectors, word1, word2)
    return get_vector(cwv, word1)'*get_vector(cwv, word2)
end

"""
    cosine_similar_words(cwv, word, n=10)

Return the top `n` (by default `n = 10`) most similar words to `word`
from the CompressedWordVectors `cwv`.
"""
function cosine_similar_words(cwv::CompressedWordVectors, word, n=10)
    indx, metr = cosine(cwv, word, n)
    return vocabulary(cwv)[indx]
end

"""
    analogy(cwv, pos, neg, n=5)

Compute the analogy similarity between two lists of words. The positions
and the similarity values of the top `n` similar words will be returned.
For example,
`king - man + woman = queen` will be
`pos=[\"king\", \"woman\"], neg=[\"man\"]`.
"""
function analogy(cwv::CompressedWordVectors{Q,U,D,T,S,H},
                 pos::AbstractArray, neg::AbstractArray, n=5
                ) where {Q,U,D,T,S,H}
    m, n_vocab = size(cwv)
    n_pos = length(pos)
    n_neg = length(neg)
    anal_vecs = Matrix{T}(undef, m, n_pos + n_neg)

    for (i, word) in enumerate(pos)
        anal_vecs[:,i] = get_vector(cwv, word)
    end
    for (i, word) in enumerate(neg)
        anal_vecs[:,i+n_pos] = -get_vector(cwv, word)
    end
    mean_vec = mean(anal_vecs, dims=2)
    metrics = cwv.vectors'*mean_vec
    top_positions = sortperm(metrics[:], rev = true)[1:n+n_pos+n_neg]
    for word in [pos;neg]
        idx = index(cwv, word)
        loc = findfirst(x->x==idx, top_positions)
        if loc != nothing
            splice!(top_positions, loc)
        end
    end
    topn_positions = top_positions[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    analogy_words(cwv, pos, neg, n=5)

Return the top `n` words computed by analogy similarity between
positive words `pos` and negaive words `neg`. from the
CompressedWordVectors `cwv`.
"""
function analogy_words(cwv::CompressedWordVectors, pos, neg, n=5)
    indx, metr = analogy(cwv, pos, neg, n)
    return vocabulary(cwv)[indx]
end

"""
    compressedwordvectors(filename [,type=Float64][; kind=:text])

Generate a `CompressedWordVectors` type object from a file.

# Arguments
  * `filename::AbstractString` the embeddings file name
  * `type::Type` type of the embedding vector elements; default `Float64`

# Keyword arguments
  * `kind::Symbol` specifies whether the embeddings file is textual (`:text`)
or binary (`:binary`); default `:text`
"""
function compressedwordvectors(filename::AbstractString,
                               ::Type{T};
                               kind::Symbol=:text) where T <: Real
    if kind == :binary
        return _from_binary(T, filename)
    elseif kind == :text
        return _from_text(T, filename)
    else
        throw(ArgumentError("Unknown embedding file kind $(kind)"))
    end
end

compressedwordvectors(filename::AbstractString; kind::Symbol=:text) =
    compressedwordvectors(filename, Float64, kind=kind)

# Generate a WordVectors object from binary file
function _from_binary(::Type{T}, filename::AbstractString) where T<:Real
    open(filename) do fid
        nrows, vocab_size = map(x -> parse(Int, x), split(readline(fid), ' '))
        d, k, m = map(x -> parse(Int, x), split(readline(fid), ' '))
        _module, _val = split(readline(fid), ".")
        Q = eval(Expr(:., Symbol(_module), QuoteNode(Symbol(_val))))
        _module, _val = split(readline(fid), ".")
        D = eval(Expr(:., Symbol(_module), QuoteNode(Symbol(_val))))
        U = eval(Symbol(readline(fid)))
        T0 = eval(Symbol(readline(fid)))

        # Read vocabulary and compressed data
        vocab = Vector{String}(undef, vocab_size)
        data = zeros(U, m, vocab_size)
        binary_length = sizeof(U) * m
        for i in 1:vocab_size
            vocab[i] = strip(readuntil(fid, ' '))
            data[:,i] = collect(reinterpret(U, read(fid, binary_length)))
        end

        # Read codebooks
        cbooks = Vector{CodeBook{U,T}}(undef, m)
        codes_length = sizeof(U) * k
        codevecs_length = sizeof(T0) * k
        for i = 1:m
            codes = collect(reinterpret(U, read(fid, codes_length)))
            vectors = zeros(T, d, k)
            for j in 1:d
                vectors[j,:] = T.(reinterpret(T0, read(fid, codevecs_length)))
            end
            cbooks[i] = CodeBook(codes, vectors)
        end
        rotmat = Matrix{T}(undef, nrows, nrows)
        binary_length = sizeof(T) * nrows
        for i in 1:nrows
            rotmat[:,i] = collect(reinterpret(T, read(fid, binary_length)))
        end
        quantizer = ArrayQuantizer(Q(), (nrows, vocab_size), cbooks, k, D(), rotmat)
        qa = QuantizedArray(quantizer, data)
        return CompressedWordVectors(vocab, qa)
    end
end

# Generate a WordVectors object from text file
function _from_text(::Type{T}, filename::AbstractString) where T<:Real
    open(filename) do fid
        nrows, vocab_size = map(x -> parse(Int, x), split(readline(fid), ' '))
        d, k, m = map(x -> parse(Int, x), split(readline(fid), ' '))
        _module, _val = split(readline(fid), ".")
        Q = eval(Expr(:., Symbol(_module), QuoteNode(Symbol(_val))))
        _module, _val = split(readline(fid), ".")
        D = eval(Expr(:., Symbol(_module), QuoteNode(Symbol(_val))))
        U = eval(Symbol(readline(fid)))
        T0 = eval(Symbol(readline(fid)))

        # Read vocabulary and compressed data
        vocab = Vector{String}(undef, vocab_size)
        data = zeros(U, m, vocab_size)
        for i in 1:vocab_size
            line = readline(fid)
            parts = split(line, ' ')
            vocab[i] = parts[1]
            data[:,i] = map(x->parse(U, x), parts[2:end])
        end

        # Read codebooks
        cbooks = Vector{CodeBook{U,T}}(undef, m)
        for i = 1:m
            codes = map(x->parse(U, x), split(readline(fid),' '))
            vectors = zeros(T, d, k)
            for j in 1:d
                vectors[j,:] = map(x->parse(T, x), split(readline(fid), ' '))
            end
            cbooks[i] = CodeBook(codes, vectors)
        end
        rotmat = Matrix{T}(undef, nrows, nrows)
        for i in 1:nrows
            rotmat[:,i] = map(x->parse(T, x), split(readline(fid), ' '))
        end
        quantizer = ArrayQuantizer(Q(), (nrows, vocab_size), cbooks, k, D(), rotmat)
        qa = QuantizedArray(quantizer, data)
        return CompressedWordVectors(vocab, qa)
    end
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
while true
    initialize(_PhysicalDeviceFeatures2, _PhysicalDeviceVulkan12Features, _PhysicalDeviceVulkanMemoryModelFeatures)
    f1 = initialize(PhysicalDeviceFeatures2, PhysicalDeviceVulkan12Features, PhysicalDeviceVulkanMemoryModelFeatures)

    chain_args = [PhysicalDeviceFeatures2(PhysicalDeviceFeatures()), PhysicalDeviceVulkan12Features(), PhysicalDeviceVulkanMemoryModelFeatures(false, false, false)]
    chained = chain(chain_args...)
    unchain(chained)

    f2 = _PhysicalDeviceFeatures2(f1)
    from_vk(PhysicalDeviceFeatures2, f2.vks, PhysicalDeviceVulkan12Features, PhysicalDeviceVulkanMemoryModelFeatures)
    f3 = Base.unsafe_convert(VkCore.VkPhysicalDeviceFeatures2, f2)
    f4 = PhysicalDeviceFeatures2(f3, PhysicalDeviceVulkan12Features, PhysicalDeviceVulkanMemoryModelFeatures)
    break
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
function plotregion(P, U, q_U, P_y, y, P_cost, y_cost, display_threshold_factor, display_reduction_factor,
    save_folder, title_string, file_name;
    save_plot_flag = true,
    display_plot_flag = true,
    canvas_size = (1000, 400))

    # # reduce the plotting positions for low signal regions. Otherwise the plot store size will be too large, and the time to load the plot will be long.
    # inds, _ = NMRSpectraSimulator.prunelowsignalentries(q_U, display_threshold_factor, display_reduction_factor)
    # P_display = P[inds]
    # U_display = U[inds]
    # q_U_display = q_U[inds]
    P_display = P
    U_display = U
    q_U_display = q_U

    # plot.
    isdir(save_folder) || mkpath(save_folder)

    plots_save_path = joinpath(save_folder, file_name)
    #title_string = "$(project_name) alignment results, region $(r), real part"

    plot_obj = Plots.plot( P_display,
        real.(q_U_display),
        title = title_string,
        label = "model",
        seriestype = :line,
        ticks = :native,
        xlims = (P_display[1],P_display[end]),
        hover = P_display,
        linewidth = 4,
        xlabel = "ppm",
        ylabel = "real part of spectrum",
        xflip = true,
        size = canvas_size)

    Plots.plot!(plot_obj, P_y, real.(y), label = "full data",
        seriestype = :line,
        linestyle = :dot,
        xflip = true,
        linewidth = 4)

    Plots.plot!(plot_obj, P_cost, real.(y_cost), label = "fit data",
        markershape = :circle,
        seriestype = :scatter,
        xflip = true)

    if save_plot_flag
        Plots.savefig(plot_obj, plots_save_path)
    end

    if display_plot_flag
        display(plot_obj)
    end

    return nothing
end

function plotquantificationresults(As, Es, ws, save_folder,
    P_y, y,
    region_min_dist,
    obj_funcs, minfs, minxs, rets,
    display_reduction_factor, display_threshold_factor,
    cost_inds_set, loop_range;
    canvas_size = (1000, 400),
    display_flag = false,
    save_plot_flag = true,
    N_viz = 50000)

    U = LinRange(u_min, u_max, N_viz)
    P = hz2ppmfunc.(U)
    U_rad = U .* (2*π)

    for r in loop_range

        y_cost = y[cost_inds_set[r]]
        P_cost = P_y[cost_inds_set[r]]

        q2 = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As, Es; w = ws[r])

        cost = obj_funcs[r](minxs[r])
        #ws[r][1] = 0.0
        q_U = q2.(U_rad)

        file_name = "quantify_real_$(r).html"
        title_string = "$(project_name) quantify results, region $(r), real part, cost = $(cost)"

        plotregion(P, U, q_U, P_y, y, P_cost, y_cost, display_threshold_factor, display_reduction_factor,
            save_folder, title_string, file_name;
            save_plot_flag = save_plot_flag,
            display_plot_flag = display_flag,
            canvas_size = canvas_size)
    end

    return nothing
end


text:
###############################################################################
#
# Shortcut to do all steps within one function
#
###############################################################################
"""
    process_RAW!(
                dt::MOD15A2Hv006LAI{FT},
                year::Int;
                days::Int = 8,
                zooms::Int = 12,
                nthread::Int = 4
    ) where {FT<:AbstractFloat}

Grid the data and save it to a csv file, given
- `param` Parameter sets to pass to different threads
- `dt` [`AbstractUngriddedData`](@ref) type ungridded data type
- `params` Array of `param`
- `nthread` Number of threads to run the code in parallel
"""
function process_RAW!(
            dt::AbstractUngriddedData{FT},
            year::Int;
            days::Int = 8,
            zooms::Int = 12,
            nthread::Int = 4
) where {FT<:AbstractFloat}
    fetch_RAW!(dt, year);
    params = query_RAW(dt, year, days);
    if nthread > 1
        grid_RAW!(dt, params, nthread);
    else
        grid_RAW!(dt, params);
    end
    compile_RAW!(dt, params, year, days, zooms);

    return nothing
end


text:
#!/usr/bin/env julia

import Celeste.ParallelRun: one_node_joint_infer, infer_init, BoundingBox, find_neighbors
import Celeste.SDSSIO: RunCamcolField, load_field_images, PlainFITSStrategy

const datadir = joinpath(Pkg.dir("Celeste"), "test", "data")
if ccall(:jl_generating_output, Cint, ()) == 0
    cd(datadir) do
        run(`make RUN=7713 CAMCOL=3 FIELD=152`)
    end
end

"""
This benchmark operates on a box of the sky that contains just
one light source. It visits 1048 pixels per evaluation of the elbo,
and the optimizer runs for 37 iterations, for 38,776 pixel-visits
in total.
"""
function benchmark_one_light_source()
    box = BoundingBox(347.7444, 347.7446, 16.6202, 16.6204)
    rcfs = [RunCamcolField(7713, 3, 152)]

    ctni = infer_init(rcfs, PlainFITSStrategy(datadir); box=box)[1:4]

    # Warm up---this compiles the code
    one_node_joint_infer(ctni...)

    # clear allocations in case julia is running with --track-allocations=user
    Profile.clear_malloc_data()

    if isempty(ARGS)
        @time one_node_joint_infer(ctni...)
    elseif ARGS[1] == "--profile"
        Profile.init(delay=0.01)
        @profile one_node_joint_infer(ctni...)
        Profile.print(format=:flat, sortedby=:count)
    end
end

benchmark_one_light_source()


text:
function _update_dual_constitutive!(H, B, C::Array{T,5}, nx, ny, nz) where T
    @turbo for k = 2:nz-1
        for j = 2:ny-1
            for i = 2:nx-1
                by = 0.25 * (B[i, j, k, 2] + B[i-1, j, k, 2] + B[i-1, j+1, k, 2] + B[i, j+1, k, 2])
                bz = 0.25 * (B[i, j, k, 3] + B[i, j, k+1, 3] + B[i-1, j, k+1, 3] + B[i-1, j, k, 3])
                H[i, j, k, 1] = C[i, j, k, 1, 1] * B[i, j, k, 1] + C[i, j, k, 1, 2] * by + C[i, j, k, 1, 3] * bz

                bx = 0.25 * (B[i, j, k, 1] + B[i, j-1, k, 1] + B[i+1, j-1, k, 1] + B[i+1, j, k, 1])
                bz = 0.25 * (B[i, j, k, 3] + B[i, j-1, k, 3] + B[i, j-1, k+1, 3] + B[i, j, k+1, 3])

                H[i, j, k, 2] = C[i, j, k, 2, 1] * bx + C[i, j, k, 2, 2] * B[i, j, k, 2] + C[i, j, k, 2, 3] * bz

                bx = 0.25 * (B[i, j, k, 1] + B[i+1, j, k, 1] + B[i+1, j, k-1, 1] + B[i, j, k-1, 1])
                by = 0.25 * (B[i, j, k, 2] + B[i, j+1, k, 2] + B[i, j+1, k-1, 2] + B[i, j, k-1, 2])

                H[i, j, k, 3] = C[i, j, k, 3, 1] * bx + C[i, j, k, 3, 2] * by + C[i, j, k, 3, 3] * B[i, j, k, 3]

            end
        end
    end

    @turbo for k = 1:1
        for j = 1:1
            for i = 2:nx-1
                by = 0.25 * (B[i, j, k, 2] + B[i-1, j, k, 2] + B[i-1, j+1, k, 2] + B[i, j+1, k, 2])
                bz = 0.25 * (B[i, j, k, 3] + B[i, j, k+1, 3] + B[i-1, j, k+1, 3] + B[i-1, j, k, 3])
                H[i, j, k, 1] = C[i, j, k, 1, 1] * B[i, j, k, 1] + C[i, j, k, 1, 2] * by + C[i, j, k, 1, 3] * bz

            end
        end
    end
    @turbo for k = 1:1
        for j = 2:ny-1
            for i = 1:1
                bx = 0.25 * (B[i, j, k, 1] + B[i, j-1, k, 1] + B[i+1, j-1, k, 1] + B[i+1, j, k, 1])
                bz = 0.25 * (B[i, j, k, 3] + B[i, j-1, k, 3] + B[i, j-1, k+1, 3] + B[i, j, k+1, 3])
                H[i, j, k, 2] = C[i, j, k, 2, 1] * bx + C[i, j, k, 2, 2] * B[i, j, k, 2] + C[i, j, k, 2, 3] * bz
            end
        end
    end
    @turbo for k = 2:nz-1
        for j = 1:1
            for i = 1:1
                bx = 0.25 * (B[i, j, k, 1] + B[i+1, j, k, 1] + B[i+1, j, k-1, 1] + B[i, j, k-1, 1])
                by = 0.25 * (B[i, j, k, 2] + B[i, j+1, k, 2] + B[i, j+1, k-1, 2] + B[i, j, k-1, 2])
                H[i, j, k, 3] = C[i, j, k, 3, 1] * bx + C[i, j, k, 3, 2] * by + C[i, j, k, 3, 3] * B[i, j, k, 3]
            end
        end
    end

end

function _update_dual_constitutive!(H, B, C::Array{T,4}, nx, ny, nz) where T
    @turbo for p in 1:3
        for k = 2:nz-1
            for j = 2:ny-1
                for i = 2:nx-1
                    H[i, j, k, p] = C[i, j, k, p] * B[i, j, k, p]
                end
            end
        end
    end

    @turbo for k = 1:1
        for j = 1:1
            for i = 2:nx-1
                H[i, j, k, 1] = C[i, j, k, 1] * B[i, j, k, 1]
            end
        end
    end
    @turbo for k = 1:1
        for j = 2:ny-1
            for i = 1:1
                H[i, j, k, 2] = C[i, j, k, 2] * B[i, j, k, 2]
            end
        end
    end
    @turbo for k = 2:nz-1
        for j = 1:1
            for i = 1:1
                H[i, j, k, 3] = C[i, j, k, 3] * B[i, j, k, 3]
            end
        end
    end
end


text:
const EAGER_INIT = Ref{Bool}(false)
const EAGER_THUNK_CHAN = Channel(typemax(Int))
const EAGER_ID_MAP = Dict{UInt64,Int}()
const EAGER_CONTEXT = Ref{Context}()
const EAGER_STATE = Ref{ComputeState}()

function eager_context()
    if !isassigned(EAGER_CONTEXT)
        EAGER_CONTEXT[] = Context([myid(),workers()...])
    end
    return EAGER_CONTEXT[]
end

function init_eager()
    EAGER_INIT[] && return
    EAGER_INIT[] = true
    ctx = eager_context()
    @async try
        sopts = SchedulerOptions(;allow_errors=true)
        topts = ThunkOptions(;single=1)
        Dagger.compute(ctx, Dagger.delayed(eager_thunk; options=topts)(); options=sopts)
    catch err
        iob = IOContext(IOBuffer(), :color=>true)
        println(iob, "Error in eager scheduler:")
        Base.showerror(iob, err)
        Base.show_backtrace(iob, catch_backtrace())
        println(iob)
        seek(iob.io, 0)
        write(stderr, iob)
    finally
        EAGER_INIT[] = false
    end
end

"Adjusts the scheduler's cached pressure indicators for the specified worker by
the specified amount, and signals the scheduler to try scheduling again if
pressure decreased."
function adjust_pressure!(h::SchedulerHandle, proc::Processor, pressure)
    uid = Dagger.get_tls().sch_uid
    lock(TASK_SYNC) do
        PROC_UTILIZATION[uid][proc][] += pressure
        notify(TASK_SYNC)
    end
    exec!(_adjust_pressure!, h, myid(), proc, pressure)
end
function _adjust_pressure!(ctx, state, task, tid, (pid, proc, pressure))
    state.worker_pressure[pid][proc] += pressure
    if pressure < 0
        put!(state.chan, RescheduleSignal())
    end
    nothing
end

"Allows a thunk to safely wait on another thunk, by temporarily reducing its
effective pressure to 0."
function thunk_yield(f)
    if Dagger.in_thunk()
        h = sch_handle()
        tls = Dagger.get_tls()
        proc = tls.processor
        util = tls.utilization
        adjust_pressure!(h, proc, -util)
        try
            f()
        finally
            adjust_pressure!(h, proc, util)
        end
    else
        f()
    end
end

function eager_thunk()
    @assert myid() == 1
    h = sch_handle()
    exec!(h) do ctx, state, task, tid, _
        EAGER_STATE[] = state
        nothing
    end
    tls = Dagger.get_tls()
    # Don't apply pressure from this thunk
    adjust_pressure!(h, tls.processor, -tls.utilization)
    while isopen(EAGER_THUNK_CHAN)
        try
            added_future, future, uid, ref, f, args, opts = take!(EAGER_THUNK_CHAN)
            # preserve inputs until they enter the scheduler
            tid = GC.@preserve args begin
                _args = map(x->x isa Dagger.EagerThunk ? ThunkID(EAGER_ID_MAP[x.uid], x.thunk_ref) : x, args)
                add_thunk!(f, h, _args...; future=future, ref=ref, opts...)
            end
            EAGER_ID_MAP[uid] = tid.id
            put!(added_future, tid.ref)
        catch err
            iob = IOContext(IOBuffer(), :color=>true)
            println(iob, "Error in eager listener:")
            Base.showerror(iob, err)
            Base.show_backtrace(iob, catch_backtrace())
            println(iob)
            seek(iob.io, 0)
            write(stderr, iob)
        end
    end
end

eager_cleanup(t::Dagger.EagerThunkFinalizer) =
    @async eager_cleanup(EAGER_STATE[], t.uid)
function eager_cleanup(state, uid)
    lock(state.lock) do
        if !haskey(EAGER_ID_MAP, uid)
            return
        end
        tid = EAGER_ID_MAP[uid]
        delete!(EAGER_ID_MAP, uid)

        # N.B. cache and errored expire automatically
        delete!(state.thunk_dict, tid)
    end
end

_find_thunk(e::Dagger.EagerThunk) =
    unwrap_weak_checked(EAGER_STATE[].thunk_dict[EAGER_ID_MAP[e.uid]])


text:
@concrete struct PeriodicFourierWorkspace
  workvector::Vector{ComplexF64}
  ik::Vector{ComplexF64}
  fft_plan
end

function PeriodicFourierWorkspace(N::Int, L::Float64)
  @assert iseven(N)
  ik = 2π * im * vcat(0:N÷2,-N÷2+1:-1) ./ L;
  return PeriodicFourierWorkspace(zeros(ComplexF64, N), ik, nothing)
end
Base.length(p::PeriodicFourierWorkspace) = length(p.ik)

struct FourierField{BC<:PeriodicGridBC, T} <: AbstractField{BC}
  chargedensity::EquispacedValueGrid{BC,T}
  electricfield::EquispacedValueGrid{BC,T}
  helper::PeriodicFourierWorkspace
  function FourierField(chargedensity::EquispacedValueGrid{BC,T},
      electricfield::EquispacedValueGrid{BC,T}) where {BC<:PeriodicGridBC,T}
    helper = PeriodicFourierWorkspace(first(size(chargedensity)), chargedensity.L)
    zero!(electricfield)
    fourier = new{BC,T}(chargedensity, electricfield, helper)
    solve!(fourier)
    return fourier
  end
end
function FourierField(chargedensity::EquispacedValueGrid{PeriodicGridBC})
  return FourierField(chargedensity, deepcopy(chargedensity))
end
function FourierField(N::Int, L::Real)
  return FourierField(EquispacedValueGrid(N, L, PeriodicGridBC))
end

zero!(f::FourierField) = (zero!(f.chargedensity), zero!(f.electricfield))

function solve!(f::FourierField)
  f.helper.workvector .= f.chargedensity
  fft!(f.helper.workvector)
  f.helper.workvector ./= f.helper.ik
  f.helper.workvector[1] *= false
  f.electricfield .= real.(ifft!(f.helper.workvector))
  return nothing
end


text:
#####################################################
### Project of IN480 - From Python 2.7 to Julia 0.6
### --- SIMPLEXN - speedup ---
### Authors: Fabio Fatelli, Emanuele Loprevite
#####################################################

include("../src/simplexn-serial.jl")
include("../src/simplexn-parallel.jl")

using Plots
Plots.scalefontsizes(0.7)
gr() # loading backend

VOID = [Int64[]],[[0]] # the empty simplicial model
nt, N = 7, 5

# Compute the execution mean time
function timing(f::Function,x,n::Int64)
	t = Array{Float64}(n)
	f(x...)
	for k in 1:n
		t[k] = @elapsed f(x...)
	end
	return mean(t)
end

# Create and save the plots in the folder "figures"
function plotting(name::String,timeS::Array{Float64,1},timeP::Array{Float64,1},flag::Int64)
	l = max(length(timeS),length(timeP))+1
	pathName = "../doc/tex/figures/"*name
	s, p, xlb, ylb, DPI = "Serial", "Parallel", "Input", "Time (seconds)", 150
	plot(timeS,xlims=(1,l),title=name,label=s,xlabel=xlb,ylabel=ylb,dpi=DPI,legend=:topleft)
	savefig(pathName*"serial")
	plot(timeP,xlims=(1,l),title=name,label=p,xlabel=xlb,ylabel=ylb,dpi=DPI,legend=:topleft)
	savefig(pathName*"parallel")
	plot([timeS,timeP],xlims=(1,l),title=name,label=[s,p],xlabel=xlb,ylabel=ylb,dpi=DPI,legend=:topleft)
	savefig(pathName*"both")
end

# Create the plots
function plotting(name::String,timeS::Array{Float64,1},timeP::Array{Float64,1})
	l = max(length(timeS),length(timeP))+1
	s, p, xlb, ylb, DPI = "Serial", "Parallel", "Input", "Time (seconds)", 150
	p1 = plot(timeS,label=s,title=name)
	p2 = plot(timeP,label=p)
	p3 = plot([timeS,timeP],label=[s,p])
	plot(p1,p2,p3,xlims=(1,l),xlabel=xlb,ylabel=ylb,dpi=1.5*DPI,legend=:topleft)
end

timeSer = [timing(larExtrude1,[VOID,repmat([1,2,3,4],k^4)],nt) for k in 1:2*N]
timePar = [timing(plarExtrude1,[VOID,repmat([1,2,3,4],k^4)],nt) for k in 1:2*N]
plotting("larExtrude1",timeSer,timePar)

simp = [collect(1:500)]
sLen = length(simp[1])
for k in 2:2*N^2
	push!(simp,simp[end]+sLen)
end
timeSer = [timing(larSimplexFacets,[simp[1:k]],nt) for k in 1:2*N^2]
timePar = [timing(plarSimplexFacets,[simp[1:k]],nt) for k in 1:2*N^2]
plotting("larSimplexFacets",timeSer,timePar)

verts, quads = [[0,0,0],[0,1,0],[1,0,0],[1,1,0],[2,2,0],[2,3,0],[3,2,0],[3,3,0]], [[0,1,2,3],[4,5,6,7]]
len = length(quads[1])
for k in 2:2*N^2
	append!(verts,verts[end-len+1:end]+1)
	append!(quads,[quads[end]+len])
end
timeSer = [timing(quads2tria,[(verts[1:len*k],quads[1:k])],nt) for k in 1:2*N^2]
timePar = [timing(pquads2tria,[(verts[1:len*k],quads[1:k])],nt) for k in 1:2*N^2]
plotting("quads2tria",timeSer,timePar)

max_stars_repo_name: EmaLoprevite/Simplexn

max_issues_repo_name: EmaLoprevite/Simplexn

max_forks_repo_name: EmaLoprevite/Simplexn


text:
using BinaryBuilder

name = "PCRE"
version = v"8.42"

# Collection of sources required to build Pcre
sources = [
    ArchiveSource("https://ftp.pcre.org/pub/pcre/pcre-$(version.major).$(version.minor).tar.bz2",
               "2cd04b7c887808be030254e8d77de11d3fe9d4505c39d4b15d2664ffe8bf9301")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pcre-*/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-utf8 --enable-unicode-properties
make -j${nproc} VERBOSE=1
make install VERBOSE=1
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libpcre", :libpcre)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using Documenter, GeoAcceleratedArrays

makedocs(;
    modules=[GeoAcceleratedArrays],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/evetion/GeoAcceleratedArrays.jl/blob/{commit}{path}#L{line}",
    sitename="GeoAcceleratedArrays.jl",
    authors="Maarten Pronk, Deltares",
    assets=String[],
)

deploydocs(;
    repo="github.com/evetion/GeoAcceleratedArrays.jl",
)


text:
function bubblesort(A::Array{Int})
    X=deepcopy(A)
    for i in 1:size(A, 1)
        for j in 1:size(A, 1)-i
            if(X[j]>X[j+1])
                X[j],X[j+1]=X[j+1],X[j]
            end
        end
    end
    return X
end

function bubblesort!(A::Array{Int})
    for i in 1:size(A, 1)
        for j in 1:size(A, 1)-i
            if(A[j]>A[j+1])
                A[j],A[j+1]=A[j+1],A[j]
            end
        end
    end
    return A
end

function bubblesortperm(A::Array{Int})
    indexes=collect(1:size(A, 1))
    for i in 1:size(A, 1)
        for j in 1:size(A, 1)-i
            if(A[j]>A[j+1])
                indexes[j],indexes[j+1]=indexes[j+1],indexes[j]
            end
        end
    end
    return indexes
end

function bubblesortperm!(A::Array{Int}, B::Array{Int})
    if(size(A, 1)==size(B, 1))
        B=collect(1:size(A, 1))
        for i in 1:size(A, 1)
            for j in 1:size(A, 1)-i
                if(A[j]>A[j+1])
                    B[j],B[j+1]=B[j+1],B[j]
                end
            end
        end
        return B
    end
end

print(bubblesortperm([2,6,4,1,7,3]))


text:
"""
```
taylor_rule()
```

creates a Taylor rule with monetary policy shocks and coefficients on the
contemporaneous inflation gap, output gap, and growth in the output gap.
"""
function taylor_rule()
    AltPolicy(:taylor_rule, taylor_rule_eqcond, taylor_rule_solve,
              forecast_init = taylor_rule_forecast_init,
              color = RGB(0., 0., 0.5430)) # dark blue
end

"""
```
taylor_replace_eq_entries(m::AbstractDSGEModel,
                          Γ0::Matrix{Float64}, Γ1::Matrix{Float64},
                          C::Vector{Float64}, Ψ::Matrix{Float64}, Π::Matrix{Float64})taylor_replace_eq_entries()
```

implements a Taylor rule with monetary policy shocks and coefficients on the
contemporaneous inflation gap, output gap, and growth in the output gap.
"""
function taylor_rule_replace_eq_entries(m::AbstractDSGEModel,
                                        Γ0::Matrix{Float64}, Γ1::Matrix{Float64},
                                        C::Vector{Float64}, Ψ::Matrix{Float64}, Π::Matrix{Float64})

    # add law of motion for pgap
    eq        = m.equilibrium_conditions
    endo      = m.endogenous_states

    # Zero out old policy rule
    Γ0[eq[:eq_mp], :]              .= 0.
    Γ1[eq[:eq_mp], :]              .= 0.
    Ψ[eq[:eq_mp], :]               .= 0.

    # Replace monetary policy rule
    Γ0[eq[:eq_mp], endo[:R_t]]      = 1.
    Γ1[eq[:eq_mp], endo[:R_t]]      = m[:ρ]
    C[ eq[:eq_mp]]                  = 0.

    Γ0[eq[:eq_mp], endo[:π_t]]      = -(1. - m[:ρ]) * m[:ψ1]
    Γ0[eq[:eq_mp], endo[:π_star_t]] = (1. - m[:ρ]) * m[:ψ1]
    Γ0[eq[:eq_mp], endo[:y_t]]      = -(1. - m[:ρ]) * m[:ψ2] - m[:ψ3]
    Γ0[eq[:eq_mp], endo[:y_f_t]]    = (1. - m[:ρ]) * m[:ψ2] + m[:ψ3]
    Γ1[eq[:eq_mp], endo[:y_t]]      = -m[:ψ3]
    Γ1[eq[:eq_mp], endo[:y_f_t]]    = m[:ψ3]

    # Add MP shocks
    Γ0[eq[:eq_mp], endo[:rm_t]]     = -1.

    return Γ0, Γ1, C, Ψ, Π
end

"""
```
taylor_rule_eqcond(m::AbstractDSGEModel)
```

Solves for the transition equation of `m` under a Taylor Rule.
"""
function taylor_rule_eqcond(m::AbstractDSGEModel, reg::Int = 1)

    # Get equilibrium condition matrices
    Γ0, Γ1, C, Ψ, Π = eqcond(m, reg)

    # Update parameters to regime `reg` since `eqcond` will reset to the first regime
    # Note that this step is necessary (in general) since
    # taylor_replace_eq_entries will use parameters as coefficients for
    # the monetary policy rule
    for para in m.parameters
        if !isempty(para.regimes)
            if (haskey(get_settings(m), :model2para_regimes) ? haskey(get_setting(m, :model2para_regime), para.key) : false)
                ModelConstructors.toggle_regime!(para, reg, get_setting(m, :model2para_regime)[para.key])
            else
                ModelConstructors.toggle_regime!(para, reg)
            end
        end
    end

    Γ0, Γ1, C, Ψ, Π = taylor_rule_replace_eq_entries(m, Γ0, Γ1, C, Ψ, Π)

    # Reset to the first regime
    ModelConstructors.toggle_regime!(m.parameters, 1)

    return Γ0, Γ1, C, Ψ, Π
end

"""
```
taylor_rule_solve(m::AbstractDSGEModel)
```

Solves for the transition equation of `m` under a Taylor Rule.
"""
function taylor_rule_solve(m::AbstractDSGEModel; regime_switching::Bool = false, regimes::Vector{Int} = Int[1])

    # Get equilibrium condition matrices
    if length(regimes) == 1
        Γ0, Γ1, C, Ψ, Π  = taylor_rule_eqcond(m, regimes[1])
        TTT_gensys, CCC_gensys, RRR_gensys, eu = gensys(Γ0, Γ1, C, Ψ, Π, 1+1e-6, verbose = :low)

        # Check for LAPACK exception, existence and uniqueness
        if eu[1] != 1 || eu[2] != 1
            throw(GensysError())
        end

        TTT_gensys = real(TTT_gensys)
        RRR_gensys = real(RRR_gensys)
        CCC_gensys = real(CCC_gensys)

        # Augment states
        TTT, RRR, CCC = DSGE.augment_states(m, TTT_gensys, RRR_gensys, CCC_gensys; regime_switching = regime_switching,
                                            reg = regimes[1])
        return TTT, RRR, CCC
    else
        Γ0s = Vector{Matrix{Float64}}(undef, length(regimes))
        Γ1s = Vector{Matrix{Float64}}(undef, length(regimes))
        Cs = Vector{Vector{Float64}}(undef, length(regimes))
        Ψs = Vector{Matrix{Float64}}(undef, length(regimes))
        Πs = Vector{Matrix{Float64}}(undef, length(regimes))
        for reg in regimes
            Γ0s[reg], Γ1s[reg], Cs[reg], Ψs[reg], Πs[reg]  = taylor_rule_eqcond(m, reg)
        end

        n_regimes = length(regimes)
        TTTs= Vector{Matrix{Float64}}(undef, n_regimes)
        RRRs = Vector{Matrix{Float64}}(undef, n_regimes)
        CCCs = Vector{Vector{Float64}}(undef, n_regimes)

        # Solve model
        for reg in regimes
            TTT_gensys, CCC_gensys, RRR_gensys, eu = gensys(Γ0s[reg], Γ1s[reg], Cs[reg], Ψs[reg], Πs[reg], 1+1e-6)

            if !((eu[1] == 1) & (eu[2] == 1))
                throw(GensysError("Gensys does not give existence"))
            end
            TTT_gensys = real(TTT_gensys)
            RRR_gensys = real(RRR_gensys)
            CCC_gensys = reshape(CCC_gensys, size(CCC_gensys, 1))

            # Augment states
            TTTs[reg], RRRs[reg], CCCs[reg] = DSGE.augment_states(m, TTT_gensys, RRR_gensys, CCC_gensys;
                                                             regime_switching = regime_switching,
                                                             reg = reg)
        end
        return TTTs, RRRs, CCCs
    end
end

"""
```
init_taylor_rule_forecast(m::AbstractDSGEModel, shocks::Matrix{T}, final_state::Vector{T})
```

Adjust shocks matrix and final state vector for forecasting under a Taylor Rule.
"""
function taylor_rule_forecast_init(m::AbstractDSGEModel, shocks::Matrix{T}, final_state::Vector{T};
                                   cond_type::Symbol = :none) where {T<:AbstractFloat}
    return shocks, final_state
end