text:
using PerformanceEval
using Base.Test

max_stars_repo_name: slundberg/PerformanceEval

max_issues_repo_name: slundberg/PerformanceEval

max_forks_repo_name: slundberg/PerformanceEval


text:
using Flux, Statistics, DelimitedFiles
using Flux: Params, gradient
using Flux.Optimise: update!
using DelimitedFiles, Statistics
using Parameters: @with_kw

# This replicates the housing data example from the Knet.jl readme. Although we
# could have reused more of Flux (see the mnist example), the library's
# abstractions are very lightweight and don't force you into any particular
# strategy.

# Struct to define hyperparameters
@with_kw mutable struct Hyperparams
    lr::Float64 = 0.1		# learning rate
    split_ratio::Float64 = 0.1	# Train Test split ratio, define percentage of data to be used as Test data
end

function get_processed_data(args)

    isfile("housing.data") ||
        download("https://raw.githubusercontent.com/MikeInnes/notebooks/master/housing.data",
            "housing.data")

    rawdata = readdlm("housing.data")'

    # The last feature is our target -- the price of the house.
    split_ratio = args.split_ratio # For the train test split

    x = rawdata[1:13,:]
    y = rawdata[14:14,:]

    # Normalise the data
    x = (x .- mean(x, dims = 2)) ./ std(x, dims = 2)

    # Split into train and test sets
    split_index = floor(Int,size(x,2)*split_ratio)
    x_train = x[:,1:split_index]
    y_train = y[:,1:split_index]
    x_test = x[:,split_index+1:size(x,2)]
    y_test = y[:,split_index+1:size(x,2)]

    train_data = (x_train, y_train)
    test_data = (x_test, y_test)

    return train_data,test_data
end

# Struct to define model
mutable struct model
    W::AbstractArray
    b::AbstractVector
end

# Function to predict output from given parameters
predict(x, m) = m.W*x .+ m.b

# Mean Squared Error
meansquarederror(ŷ, y) = sum((ŷ .- y).^2)/size(y, 2)

function train(; kws...)
    # Initialize the Hyperparamters
    args = Hyperparams(; kws...)

    # Load the data
    (x_train,y_train),(x_test,y_test) = get_processed_data(args)

    # The model
    m = model((randn(1,13)),[0.])

    loss(x, y) = meansquarederror(predict(x, m), y)

    ## Training
    η = args.lr
    θ = params([m.W, m.b])

    for i = 1:500
      g = gradient(() -> loss(x_train, y_train), θ)
      for x in θ
        update!(x, -g[x]*η)
      end
      if i%100==0
          @show loss(x_train, y_train)
        end
    end

    # Predict the RMSE on the test set
    err = meansquarederror(predict(x_test, m),y_test)
    println(err)
end

cd(@__DIR__)
train()

ext: jl

max_stars_repo_name: thomascherickal/model-zoo

max_issues_repo_name: thomascherickal/model-zoo

max_forks_repo_name: thomascherickal/model-zoo


text:
# A super-simple command line program that just prints hello.
# Build this with the `commandline_app=true` flag in `ApplicationBuilder.build_app_bundle`.

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    println("Hi what's your name?")
    name = readline()
    println("Oh hi, $name. It's a pleasure to meet you.")
    println("By the way, here's the current working directory:\n'$(pwd())'")

    println("\nGoodbye! (Press enter to exit)")
    readline()
    return 0
end


text:
friction_coefficients(model::Cartpole{T,Friction}) where T = model.friction

function RoboDojo.indices_z(model::Cartpole{T,Friction}) where T
    nq = model.nq
    nc = model.nc
    q = collect(1:nq)
    γ = collect(nq .+ (1:0))
    sγ = collect(nq .+ (1:0))
    ψ = collect(nq .+ (1:nc))
    b = collect(nq + nc .+ (1:nc))
    sψ = collect(nq + nc + nc .+ (1:nc))
    sb = collect(nq + nc + nc + nc .+ (1:nc))
    IndicesZ(q, γ, sγ, ψ, b, sψ, sb)
end

RoboDojo.num_var(model::Cartpole{T,Friction}) where T = model.nq + 4 * model.nc

function RoboDojo.indices_optimization(model::Cartpole{T,Friction}) where T
    nq = model.nq
    nc = model.nc
    nz = nq + 4 * nc
    IndicesOptimization(
        nz,
        nz,
        [collect(1:0), collect(1:0)],
        [collect(1:0), collect(1:0)],
        [[collect([nq + i, nq + nc + i]), collect([nq + nc + nc + i, nq + nc + nc + nc + i])] for i = 1:nc],
        [[collect([nq + i, nq + nc + i]), collect([nq + nc + nc + i, nq + nc + nc + nc + i])] for i = 1:nc],
        collect(1:(nq + 2 * nc)),
        collect(1:0),
        collect(nq + 2 * nc .+ (1:(2 * nc))),
        [collect(nq + 2 * nc + (i - 1) * 2 .+ (1:2)) for i = 1:nc],
        collect(nq + 2 * nc .+ (1:(2 * nc))))
end

function RoboDojo.initialize_z!(z, model::Cartpole{T,Friction}, idx::IndicesZ, q) where T
    z[idx.q] .= q
    z[idx.ψ] .= 1.0
    z[idx.b] .= 0.1
    z[idx.sψ] .= 1.0
    z[idx.sb] .= 0.1
end

RoboDojo.nominal_configuration(model::Cartpole) = zeros(model.nq)


text:
module CommonMsgs
include("geometry_msgs.jl")
export geometry_msgs
include("actionlib_msgs.jl")
export actionlib_msgs
include("diagnostic_msgs.jl")
export diagnostic_msgs
include("nav_msgs.jl")
export nav_msgs
include("sensor_msgs.jl")
export sensor_msgs
include("shape_msgs.jl")
export shape_msgs
include("stereo_msgs.jl")
export stereo_msgs
include("trajectory_msgs.jl")
export trajectory_msgs
include("visualization_msgs.jl")
export visualization_msgs
end


text:
using Random
using StatsBase
using Distributions

import StatsBase.mean

#increment is the factor by which bin increases (casted to integer)
#normalize at end
#interpret P[i] as number of elements in bin i (can be normalized in advance ...)
#this is approximate because int P(x)dx ~ sum P(x) delta x (what we do here)
"""
calculate a log-binned distribution from an evenly space distribution `P(x)`
where bins increase with `inrement_factor`

# Remark
`P` hast to be a probability density
"""
function logbin(x, P; increment_factor=sqrt(2), bin_width_start=1)
    @assert length(x)==length(P)

    dx = x[2]-x[1]

    counter   = 0
    index     = 1
    x_bin = zeros(1)
    P_bin = zeros(1)
    w_bin = ones(1)*floor(Int, bin_width_start)

    bin_width::Float64 = w_bin[end]
    for (i,p) in enumerate(P)
        #print(i, " ", index, " ", x[i], " ", P[i], " ", x_bin[index], " ", P_bin[index], " ",  w_bin[index], "\n")
        # window full; create new window
        if counter == w_bin[end]
            index += 1
            counter = 0
            bin_width *= increment_factor
            push!(x_bin, 0)
            push!(P_bin, 0)
            push!(w_bin, floor(Int, bin_width))
        end

        counter += 1
        # towards expectation value of x in log bin
        x_bin[index] += x[i]*P[i]
        # probabilities of values over bin (sum of equally spaced probabilities)
        P_bin[index] += P[i]
    end
    # normalize x bins (end of expectation value)
    x_bin ./= P_bin

    # transform P into propert probability DENSITY such that <x>=\sum xP(x)dx/sum P(x)dx
    P_bin ./= w_bin

    return x_bin, P_bin, w_bin.*dx
end
logbin(dist; increment_factor=sqrt(2), bin_width_start=1) = logbin(collect(dist.edges[1])[1:end-1], dist.weights, increment_factor=increment_factor, bin_width_start=bin_width_start)

function calculate_cdf(dist)
    x = collect(dist.edges[1])
    dx = x[2]-x[1]
    cdf = similar(x, Float64)
    cdf[1] = 0.0
    for i in 2:length(cdf)
        cdf[i] = cdf[i-1] + dist.weights[i-1]*dx
    end
    return x, cdf
end

###############################################################################
# convenient Histogram access as hist[value] or even hist[value] += number entries
Base.getindex(h::AbstractHistogram{T,1}, x::Real) where {T} = getindex(h, (x,))
function Base.getindex(h::Histogram{T,N}, xs::NTuple{N,Real}) where {T,N}
    idx = StatsBase.binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        return @inbounds h.weights[idx...]
    else
        return missing
    end
end

Base.setindex!(h::AbstractHistogram{T,1}, value::Real, x::Real) where {T} = setindex!(h, value, (x,))
function Base.setindex!(h::Histogram{T,N}, value::Real, xs::NTuple{N,Real}) where {T,N}
    h.isdensity && error("Density histogram must have float-type weights")
    idx = StatsBase.binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        @inbounds h.weights[idx...] = value
    end
end


text:
module SFEM

    using LinearAlgebra
    using DelimitedFiles
    using Pardiso
    using SparseArrays
    using WriteVTK
    #using Gadfly
    #using Cairo
    #using Fontconfig
    #using DataFrames
    #using CSV
    #using PlotlyJS
    #using ORCA

    include("G://Project-SFEM//3_FS_FEM(3D)//area.jl")
    include("G://Project-SFEM//3_FS_FEM(3D)//vectorout.jl")
    include("G://Project-SFEM//3_FS_FEM(3D)//vectorin.jl")
    include("G://Project-SFEM//3_FS_FEM(3D)//volume.jl")
    include("G://Project-SFEM//3_FS_FEM(3D)//assemblehex.jl")
    include("G://Project-SFEM//3_FS_FEM(3D)//assembletet.jl")

    export IOsfem
    export Ksfem
    export MKL
    #export Draw_mesh
    #export Draw_disp
    export WriteDISP
    #export Compare

    function IOsfem(A, B, C, D, E)

        Cen = readdlm(A)
        Cen = Array{Float64,2}(Cen)

        Ele = readdlm(B)
        Ele = Array{Int64,2}(Ele)

        Face_in = readdlm(C)
        Face_in = Array{Int64,2}(Face_in)

        Face_out = readdlm(D)
        Face_out = Array{Int64,2}(Face_out)

        Node = readdlm(E)
        Node = Array{Float64,2}(Node)

        return Node, Ele, Cen, Face_out, Face_in
    end

    function Ksfem(Node, Ele, Cen, Face_out, Face_in, EE, NU, FF)

        Q = zeros(Float64,size(Node,1)*3,1)
        max = maximum(Node,dims = 1)
        min = minimum(Node,dims = 1)
        force_point = Vector{Int64}()

        for i = 1:size(Node,1)
            if Node[i,3] == max[3] && Node[i,1] != min[1]
                force_point = append!(force_point,i)
            end
        end

        for i = force_point
            Q[i*3] = FF/size(force_point,1)
        end
        IK = Vector{Int64}(); JK = Vector{Int64}(); VK = Vector{Float64}()

        for i = 1:size(Face_out,1)
            v = Volume(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                       Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                       Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                       Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
            A1 = area(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                      Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                      Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3])
            A2 = area(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                      Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                      Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
            A3 = area(Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                      Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                      Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
            A4 = area(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                      Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                      Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
            n = Vectorout(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                          Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                          Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                          Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
            bAx = (1/v)*(n[1,1]*(1/3)*A1+n[2,1]*(5/12)*A2+n[3,1]*(1/12)*A3+n[4,1]*(5/12)*A4)
            bAy = (1/v)*(n[1,2]*(1/3)*A1+n[2,2]*(5/12)*A2+n[3,2]*(1/12)*A3+n[4,2]*(5/12)*A4)
            bAz = (1/v)*(n[1,3]*(1/3)*A1+n[2,3]*(5/12)*A2+n[3,3]*(1/12)*A3+n[4,3]*(5/12)*A4)
            bA = [bAx 0 0; 0 bAy 0; 0 0 bAz; bAy bAx 0; 0 bAz bAy; bAz 0 bAx]
            bBx = (1/v)*(n[1,1]*(1/3)*A1+n[2,1]*(5/12)*A2+n[3,1]*(5/12)*A3+n[4,1]*(1/12)*A4)
            bBy = (1/v)*(n[1,2]*(1/3)*A1+n[2,2]*(5/12)*A2+n[3,2]*(5/12)*A3+n[4,2]*(1/12)*A4)
            bBz = (1/v)*(n[1,3]*(1/3)*A1+n[2,3]*(5/12)*A2+n[3,3]*(5/12)*A3+n[4,3]*(1/12)*A4)
            bB = [bBx 0 0; 0 bBy 0; 0 0 bBz; bBy bBx 0; 0 bBz bBy; bBz 0 bBx]
            bCx = (1/v)*(n[1,1]*(1/3)*A1+n[2,1]*(1/12)*A2+n[3,1]*(5/12)*A3+n[4,1]*(5/12)*A4)
            bCy = (1/v)*(n[1,2]*(1/3)*A1+n[2,2]*(1/12)*A2+n[3,2]*(5/12)*A3+n[4,2]*(5/12)*A4)
            bCz = (1/v)*(n[1,3]*(1/3)*A1+n[2,3]*(1/12)*A2+n[3,3]*(5/12)*A3+n[4,3]*(5/12)*A4)
            bC = [bCx 0 0; 0 bCy 0; 0 0 bCz; bCy bCx 0; 0 bCz bCy; bCz 0 bCx]
            bDx = (1/v)*(n[2,1]*(1/12)*A2+n[3,1]*(1/12)*A3+n[4,1]*(1/12)*A4)
            bDy = (1/v)*(n[2,2]*(1/12)*A2+n[3,2]*(1/12)*A3+n[4,2]*(1/12)*A4)
            bDz = (1/v)*(n[2,3]*(1/12)*A2+n[3,3]*(1/12)*A3+n[4,3]*(1/12)*A4)
            bD = [bDx 0 0; 0 bDy 0; 0 0 bDz; bDy bDx 0; 0 bDz bDy; bDz 0 bDx]
            B = hcat(bA,bB,bC,bD)
            D = (EE/((1+NU)*(1-2*NU)))*[1-NU NU NU 0 0 0; NU 1-NU NU 0 0 0; NU NU 1-NU 0 0 0; 0 0 0 (1-2*NU)/2 0 0; 0 0 0 0 (1-2*NU)/2 0; 0 0 0 0 0 (1-2*NU)/2]
            ke = B'*D*B*v
            IK, JK, VK = Assembletet(IK,JK,VK,ke,Face_out[i,1],Face_out[i,2],Face_out[i,3],Face_out[i,5])
        end
        for i = 1:size(Face_in,1)
            v1 = Volume(Node[Face_in[i,1],1],Node[Face_in[i,1],2],Node[Face_in[i,1],3],
                        Node[Face_in[i,2],1],Node[Face_in[i,2],2],Node[Face_in[i,2],3],
                        Node[Face_in[i,3],1],Node[Face_in[i,3],2],Node[Face_in[i,3],3],
                        Cen[Face_in[i,4],1],Cen[Face_in[i,4],2],Cen[Face_in[i,4],3])
            v2 = Volume(Node[Face_in[i,1],1],Node[Face_in[i,1],2],Node[Face_in[i,1],3],
                        Node[Face_in[i,2],1],Node[Face_in[i,2],2],Node[Face_in[i,2],3],
                        Node[Face_in[i,3],1],Node[Face_in[i,3],2],Node[Face_in[i,3],3],
                        Cen[Face_in[i,5],1],Cen[Face_in[i,5],2],Cen[Face_in[i,5],3])
            v = v1 + v2
            A1 = area(Node[Face_in[i,1],1],Node[Face_in[i,1],2],Node[Face_in[i,1],3],
                      Node[Face_in[i,2],1],Node[Face_in[i,2],2],Node[Face_in[i,2],3],
                      Cen[Face_in[i,4],1],Cen[Face_in[i,4],2],Cen[Face_in[i,4],3])
            A2 = area(Node[Face_in[i,2],1],Node[Face_in[i,2],2],Node[Face_in[i,2],3],
                      Node[Face_in[i,3],1],Node[Face_in[i,3],2],Node[Face_in[i,3],3],
                      Cen[Face_in[i,4],1],Cen[Face_in[i,4],2],Cen[Face_in[i,4],3])
            A3 = area(Node[Face_in[i,1],1],Node[Face_in[i,1],2],Node[Face_in[i,1],3],
                      Node[Face_in[i,3],1],Node[Face_in[i,3],2],Node[Face_in[i,3],3],
                      Cen[Face_in[i,4],1],Cen[Face_in[i,4],2],Cen[Face_in[i,4],3])
            A4 = area(Node[Face_in[i,1],1],Node[Face_in[i,1],2],Node[Face_in[i,1],3],
                      Node[Face_in[i,2],1],Node[Face_in[i,2],2],Node[Face_in[i,2],3],
                      Cen[Face_in[i,5],1],Cen[Face_in[i,5],2],Cen[Face_in[i,5],3])
            A5 = area(Node[Face_in[i,2],1],Node[Face_in[i,2],2],Node[Face_in[i,2],3],
                      Node[Face_in[i,3],1],Node[Face_in[i,3],2],Node[Face_in[i,3],3],
                      Cen[Face_in[i,5],1],Cen[Face_in[i,5],2],Cen[Face_in[i,5],3])
            A6 = area(Node[Face_in[i,1],1],Node[Face_in[i,1],2],Node[Face_in[i,1],3],
                      Node[Face_in[i,3],1],Node[Face_in[i,3],2],Node[Face_in[i,3],3],
                      Cen[Face_in[i,5],1],Cen[Face_in[i,5],2],Cen[Face_in[i,5],3])
            n = Vectorin(Node[Face_in[i,1],1],Node[Face_in[i,1],2],Node[Face_in[i,1],3],
                         Node[Face_in[i,2],1],Node[Face_in[i,2],2],Node[Face_in[i,2],3],
                         Node[Face_in[i,3],1],Node[Face_in[i,3],2],Node[Face_in[i,3],3],
                         Cen[Face_in[i,4],1],Cen[Face_in[i,4],2],Cen[Face_in[i,4],3],
                         Cen[Face_in[i,5],1],Cen[Face_in[i,5],2],Cen[Face_in[i,5],3])
            bCx = (1/v)*(n[1,1]*(5/12)*A1+n[2,1]*(1/12)*A2+n[3,1]*(5/12)*A3+n[4,1]*(5/12)*A4+n[5,1]*(1/12)*A5+n[6,1]*(5/12)*A6)
            bCy = (1/v)*(n[1,2]*(5/12)*A1+n[2,2]*(1/12)*A2+n[3,2]*(5/12)*A3+n[4,2]*(5/12)*A4+n[5,2]*(1/12)*A5+n[6,2]*(5/12)*A6)
            bCz = (1/v)*(n[1,3]*(5/12)*A1+n[2,3]*(1/12)*A2+n[3,3]*(5/12)*A3+n[4,3]*(5/12)*A4+n[5,3]*(1/12)*A5+n[6,3]*(5/12)*A6)
            bC = [bCx 0 0; 0 bCy 0; 0 0 bCz; bCy bCx 0; 0 bCz bCy; bCz 0 bCx]
            bDx = (1/v)*(n[1,1]*(5/12)*A1+n[2,1]*(5/12)*A2+n[3,1]*(1/12)*A3+n[4,1]*(5/12)*A4+n[5,1]*(5/12)*A5+n[6,1]*(1/12)*A6)
            bDy = (1/v)*(n[1,2]*(5/12)*A1+n[2,2]*(5/12)*A2+n[3,2]*(1/12)*A3+n[4,2]*(5/12)*A4+n[5,2]*(5/12)*A5+n[6,2]*(1/12)*A6)
            bDz = (1/v)*(n[1,3]*(5/12)*A1+n[2,3]*(5/12)*A2+n[3,3]*(1/12)*A3+n[4,3]*(5/12)*A4+n[5,3]*(5/12)*A5+n[6,3]*(1/12)*A6)
            bD = [bDx 0 0; 0 bDy 0; 0 0 bDz; bDy bDx 0; 0 bDz bDy; bDz 0 bDx]
            bBx = (1/v)*(n[1,1]*(1/12)*A1+n[2,1]*(5/12)*A2+n[3,1]*(5/12)*A3+n[4,1]*(1/12)*A4+n[5,1]*(5/12)*A5+n[6,1]*(5/12)*A6)
            bBy = (1/v)*(n[1,2]*(1/12)*A1+n[2,2]*(5/12)*A2+n[3,2]*(5/12)*A3+n[4,2]*(1/12)*A4+n[5,2]*(5/12)*A5+n[6,2]*(5/12)*A6)
            bBz = (1/v)*(n[1,3]*(1/12)*A1+n[2,3]*(5/12)*A2+n[3,3]*(5/12)*A3+n[4,3]*(1/12)*A4+n[5,3]*(5/12)*A5+n[6,3]*(5/12)*A6)
            bB = [bBx 0 0; 0 bBy 0; 0 0 bBz; bBy bBx 0; 0 bBz bBy; bBz 0 bBx]
            bAx = (1/v)*(n[1,1]*(1/12)*A1+n[2,1]*(1/12)*A2+n[3,1]*(1/12)*A3)
            bAy = (1/v)*(n[1,2]*(1/12)*A1+n[2,2]*(1/12)*A2+n[3,2]*(1/12)*A3)
            bAz = (1/v)*(n[1,3]*(1/12)*A1+n[2,3]*(1/12)*A2+n[3,3]*(1/12)*A3)
            bA = [bAx 0 0; 0 bAy 0; 0 0 bAz; bAy bAx 0; 0 bAz bAy; bAz 0 bAx]
            bEx = (1/v)*(n[4,1]*(1/12)*A4+n[5,1]*(1/12)*A5+n[6,1]*(1/12)*A6)
            bEy = (1/v)*(n[4,2]*(1/12)*A4+n[5,2]*(1/12)*A5+n[6,2]*(1/12)*A6)
            bEz = (1/v)*(n[4,3]*(1/12)*A4+n[5,3]*(1/12)*A5+n[6,3]*(1/12)*A6)
            bE = [bEx 0 0; 0 bEy 0; 0 0 bEz; bEy bEx 0; 0 bEz bEy; bEz 0 bEx]
            B = hcat(bC,bD,bB,bA,bE)
            D = (EE/((1+NU)*(1-2*NU)))*[1-NU NU NU 0 0 0; NU 1-NU NU 0 0 0; NU NU 1-NU 0 0 0; 0 0 0 (1-2*NU)/2 0 0; 0 0 0 0 (1-2*NU)/2 0; 0 0 0 0 0 (1-2*NU)/2]
            ke = B'*D*B*v
            IK, JK, VK = Assemblehex(IK,JK,VK,ke,Face_in[i,1],Face_in[i,2],Face_in[i,3],Face_in[i,6],Face_in[i,7])
        end
        KE = sparse(IK, JK, VK)
        min = minimum(Node, dims = 1)
        node_dex = Vector{Int64}()
        @inbounds for i = 1:size(Node, 1)
            if Node[i,1] == min[1]
                append!(node_dex, i)
            end
        end
        @inbounds for i = 1:size(node_dex, 1)
            KE[node_dex[i]*3, node_dex[i]*3] = KE[node_dex[i]*3, node_dex[i]*3]*1E20
            KE[node_dex[i]*3-1, node_dex[i]*3-1] = KE[node_dex[i]*3-1, node_dex[i]*3-1]*1E20
            KE[node_dex[i]*3-2, node_dex[i]*3-2] = KE[node_dex[i]*3-2, node_dex[i]*3-2]*1E20
        end

        return KE, Q
    end
    function MKL(KE,Q)
        ps = MKLPardisoSolver()
        set_nprocs!(ps,24)
        result = similar(Q)
        result = solve(ps, KE, Q)
        println(get_nprocs(ps))

        Disp = zeros(Float64,Int(size(Q,1)/3),3)
        for i = 1:Int(size(Q,1)/3)
            Disp[i,1] = result[i*3-2,1]
            Disp[i,2] = result[i*3-1,1]
            Disp[i,3] = result[i*3,1]
        end

        return Disp
    end
    function Draw_mesh(Node,Disp,B)
        show_point = Vector{Int64}()
        for i = 1:size(Node,1)
            if Node[i,2] == 0 && Node[i,3] == 0
                show_point = append!(show_point,i)
            end
        end
        show_point_xvalue = Vector{Float64}()
        for i = show_point
            show_point_xvalue = append!(show_point_xvalue,Node[i,1])
        end
        sort_point = sortperm(show_point_xvalue)
        ax = Vector{Float64}()
        for i = sort_point
            ax = append!(ax,Node[show_point[i],1])
        end
        ay = Vector{Float64}()
        for i = sort_point
            ay = append!(ay,Disp[show_point[i],3])
        end

        face = readdlm(B)
        face_num = Int(face[1,1])
        F1 = Vector{Int64}();F2 = Vector{Int64}();F3 = Vector{Int64}()
        for i = 1:face_num
                F1 = append!(F1,Int(face[i+1,2])+1)
                F2 = append!(F2,Int(face[i+1,3])+1)
                F3 = append!(F3,Int(face[i+1,4])+1)
        end
        Face = hcat(F1,F2,F3)

        data = GenericTrace[]
        for i = 1:size(Face,1)
            X = [Node[Face[i,1],1],Node[Face[i,2],1],Node[Face[i,3],1],Node[Face[i,1],1]]
            Y = [Node[Face[i,1],2],Node[Face[i,2],2],Node[Face[i,3],2],Node[Face[i,1],2]]
            Z = [Node[Face[i,1],3],Node[Face[i,2],3],Node[Face[i,3],3],Node[Face[i,1],3]]
            trace = scatter3d(;x=X,y=Y,z=Z,mode="lines",showlegend=false,line=attr(color="#1f77b4", width=1))
            push!(data,trace)
        end

        layout = Layout(title="Beam",
    					autosize="false",
    					width=1280,
    					height=720,
    					margin=attr(l=0,r=0,b=50,t=50),
                        scene=attr(aspectmode="manual",
                                   aspectratio=attr(x=1,y=0.17,z=0.17),
                                   xaxis=attr(title="X",range=[-1,6],showgrid=false,zeroline=false,showline=false,showticklabels=false),
                                   yaxis=attr(title="Y",range=[-1,1.5],showgrid=false,zeroline=false,showline=false,showticklabels=false),
                                   zaxis=attr(title="Z",range=[-1,1.5],showgrid=false,zeroline=false,showline=false,showticklabels=false))
                        )

        plot(data,layout)
        #savehtml(p,"1.html")
    end

    function Draw_disp(Node,Disp,F)
        face = readdlm(F)
        face_num = Int(face[1,1])
        F1 = Vector{Int64}();F2 = Vector{Int64}();F3 = Vector{Int64}()
        for i = 1:face_num
                F1 = append!(F1,Int(face[i+1,2])+1)
                F2 = append!(F2,Int(face[i+1,3])+1)
                F3 = append!(F3,Int(face[i+1,4])+1)
        end
        Face = hcat(F1,F2,F3)

        for i = 1:size(Node,1)
            Node[i,1] = Node[i,1] + Disp[i,1]
            Node[i,2] = Node[i,2] + Disp[i,2]
            Node[i,3] = Node[i,3] + Disp[i,3]
        end

        data = GenericTrace[]
        for i = 1:size(Face,1)
            X = [Node[Face[i,1],1],Node[Face[i,2],1],Node[Face[i,3],1],Node[Face[i,1],1]]
            Y = [Node[Face[i,1],2],Node[Face[i,2],2],Node[Face[i,3],2],Node[Face[i,1],2]]
            Z = [Node[Face[i,1],3],Node[Face[i,2],3],Node[Face[i,3],3],Node[Face[i,1],3]]
            trace = scatter3d(;x=X,y=Y,z=Z,mode="lines",showlegend=false,line=attr(color="#1f77b4", width=1.5))
            push!(data, trace)
        end

        layout = Layout(title="Beam",
                        autosize="false",
                        width=2000,
                        height=1080,
                        margin=attr(l=0,r=0,b=0,t=0),
                        scene=attr(aspectmode="manual",
                                   aspectratio=attr(x=1,y=1/6,z=1/6),
                                   xaxis=attr(title="X",range=[0,6.2],showgrid=true,zeroline=false,showline=true,showticklabels=true,
                                   yaxis=attr(title="Y",range=[0,1.0],showgrid=true,zeroline=false,showline=true,showticklabels=true),
                                   zaxis=attr(title="Z",range=[-0.03,1.0],showgrid=true,zeroline=false,showline=true,showticklabels=true))
                                   )
                        )
        plot(data,layout)
    end
    function WriteDISP(Node,Ele,result)

        Disp = zeros(Float64,size(Node,1),3)
        for i = 1:size(Node,1)
            Disp[i,1] = result[i*3-2,1]
            Disp[i,2] = result[i*3-1,1]
            Disp[i,3] = result[i*3,1]
        end

        pts = Node
        for i = 1:size(pts, 1)
            pts[i,1] = Node[i,1] + Disp[i,1]
            pts[i,2] = Node[i,2] + Disp[i,2]
            pts[i,3] = Node[i,3] + Disp[i,3]
        end

        pts = transpose(pts)
        pdata = Array{Float64}(undef,size(pts, 2))
        for i = 1:size(pts, 2)
            pdata[i] = Disp[i,3]
        end

        celltype = VTKCellTypes.VTK_TETRA
        cells = MeshCell[]

        for i = 1:size(Ele, 1)
            inds = Array{Int64}(undef, 4)
            inds[1] = Ele[i,1]
            inds[2] = Ele[i,2]
            inds[3] = Ele[i,3]
            inds[4] = Ele[i,4]

            c = MeshCell(celltype, inds)
            push!(cells, c)
        end

        vtk = vtk_grid("3D_disp", pts, cells,compress=3)
        vtk_point_data(vtk, pdata, "Z Disp (m)")
        outfiles =  vtk_save(vtk)
        return outfiles::Vector{String}, Disp
    end
    function Compare(Node,result,EE,FF)
        Disp = zeros(Float64,size(Node,1),3)
        for i = 1:size(Node,1)
            Disp[i,1] = result[i*3-2,1]
            Disp[i,2] = result[i*3-1,1]
            Disp[i,3] = result[i*3,1]
        end

        l1 = Vector{Int64}(); l2 = Vector{Int64}(); l3 = Vector{Int64}(); l4 = Vector{Int64}()
        max = maximum(Node, dims = 1)
        min = minimum(Node, dims = 1)
        println(max)
        for i = 1:size(Node, 1)
            if Node[i,2] == min[2] && Node[i,3] == min[3]
                append!(l1, i)
            end
            if Node[i,2] == max[2] && Node[i,3] == min[3]
                append!(l2,i)
            end
            if Node[i,2] == max[2] && Node[i,3] == max[3]
                append!(l3,i)
            end
            if Node[i,2] == min[2] && Node[i,3] == max[3]
                append!(l4,i)
            end
        end
        println(l1)

        XX = Array{Float64}(undef,size(l1, 1))
        for i = 1:size(l1, 1)
            XX[i] = Node[i,1]
        end
        XX = sort(XX)

        v1 = Vector{Float64}(); v2 = Vector{Float64}(); v3 = Vector{Float64}(); v4 = Vector{Float64}()
        for i = size(l1,1)
            append!(v1,Disp[l1[i],3])
        end
        v1 = sort(v1,rev=true)

        for i = size(l2,1)
            append!(v2,Disp[l2[i],3])
        end
        v2 = sort(v2,rev=true)

        for i = size(l3,1)
            append!(v3,Disp[l3[i],3])
        end
        v3 = sort(v3,rev=true)

        for i = size(l4,1)
            append!(v4,Disp[l4[i],3])
        end
        v4 = sort(v4,rev=true)

        Y1 = (v1+v2+v3+v4)/4
        Y2 = map(i -> (FF*i^2)/(24*EE*((max[2]-min[2])*(max[3]-min[3])^3/12))*(i^2-4*(max[1]-min[1])*i+6*(max[1]-min[1])^2),XX)
        line1 = layer(x=XX,y=Y1,Geom.point,Geom.smooth,Theme(default_color="green"))
        line2 = layer(x=XX,y=Y2,Geom.point,Geom.smooth,Theme(default_color="blue"))
        plot(line1,line2)

    end

end

max_stars_repo_path: src/single-core juSFEM/SFEM.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/single-core juSFEM/SFEM.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/single-core juSFEM/SFEM.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
function ComputeFCFV_o2(mesh, sex, sey, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, tau, o2)
    nen   = 3
    ae    = zeros(mesh.nel)
    be    = zeros(mesh.nel,2)
    be_o2 = zeros(mesh.nel,2,nen)
    ze    = zeros(mesh.nel,2,2)
    pe    = zeros(mesh.nel,mesh.nf_el,nen)
    me    = zeros(mesh.nel,nen,nen)
    mei   = zeros(mesh.nel,nen,nen)
    N     = zeros(mesh.nel,mesh.nf_el,nen)
    djx   = zeros(mesh.nel,mesh.nf_el,nen)
    rjx   = zeros(mesh.nel,mesh.nf_el,nen)
    djy   = zeros(mesh.nel,mesh.nf_el,nen)
    rjy   = zeros(mesh.nel,mesh.nf_el,nen)

    # Assemble FCFV elements
    @tturbo for iel=1:mesh.nel

        be[iel,1] += mesh.vole[iel]*sex[iel]
        be[iel,2] += mesh.vole[iel]*sey[iel]

        for ifac=1:mesh.nf_el

            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel],mesh.type,mesh.ke[iel])                              # Stabilisation parameter for the face

            N[iel,ifac,1] = 1.0
            N[iel,ifac,2] = mesh.xf[nodei] - mesh.xc[iel]
            N[iel,ifac,3] = mesh.yf[nodei] - mesh.yc[iel]

            pe[iel,ifac,1] = 1.0
            pe[iel,ifac,2] = mesh.xf[nodei] - mesh.xc[iel]
            pe[iel,ifac,3] = mesh.yf[nodei] - mesh.yc[iel]

            # Assemble
            ze[iel,1,1] += (bc==1) * dAi*ni_x*VxDir[nodei] # Dirichlet
            ze[iel,1,2] += (bc==1) * dAi*ni_x*VyDir[nodei] # Dirichlet
            ze[iel,2,1] += (bc==1) * dAi*ni_y*VxDir[nodei] # Dirichlet
            ze[iel,2,2] += (bc==1) * dAi*ni_y*VyDir[nodei] # Dirichlet
            be[iel,1]   += (bc==1) * dAi*taui*VxDir[nodei] # Dirichlet
            be[iel,2]   += (bc==1) * dAi*taui*VyDir[nodei] # Dirichlet
            ae[iel]     +=           dAi*taui

        end
    end

    if o2==1 # The remaining of this function is only done if second order FCFV is used
        # me matrix
        @tturbo for iel=1:mesh.nel
            for ifac=1:mesh.nf_el

                nodei = mesh.e2f[iel,ifac]
                bc    = mesh.bc[nodei]
                dAi   = mesh.dA[iel,ifac]
                taui  = StabParam(tau,dAi,mesh.vole[iel],mesh.type,mesh.ke[iel])

                for in=1:nen
                    djx[iel,ifac,in] =  dAi * N[iel,ifac,in] * VxDir[nodei]
                    djy[iel,ifac,in] =  dAi * N[iel,ifac,in] * VyDir[nodei]
                    rjx[iel,ifac,in] =  dAi * N[iel,ifac,in]
                    rjy[iel,ifac,in] =  dAi * N[iel,ifac,in]
                    for jn=1:nen
                        me[iel,in,jn] += pe[iel,ifac,jn] * taui * dAi * N[iel,ifac,in]
                    end
                end
            end
        end

        # Inverse of me matrix
        cf11 =  1.0; cf12 =-1.0; cf13 =  1.0;
        cf21 = -1.0; cf22 = 1.0; cf23 = -1.0;
        cf31 =  1.0; cf32 =-1.0; cf33 =  1.0;
        @tturbo for iel=1:mesh.nel
            for ifac=1:mesh.nf_el
                Mm11 = me[iel,2,2] * me[iel,3,3] - me[iel,2,3] * me[iel,3,2]
                Mm12 = me[iel,2,1] * me[iel,3,3] - me[iel,2,3] * me[iel,3,1]
                Mm13 = me[iel,2,1] * me[iel,3,2] - me[iel,2,2] * me[iel,3,1]
                Mm21 = me[iel,1,2] * me[iel,3,3] - me[iel,1,3] * me[iel,3,2]
                Mm22 = me[iel,1,1] * me[iel,3,3] - me[iel,1,3] * me[iel,3,1]
                Mm23 = me[iel,1,1] * me[iel,3,2] - me[iel,1,2] * me[iel,3,1]
                Mm31 = me[iel,1,2] * me[iel,2,3] - me[iel,1,3] * me[iel,2,2]
                Mm32 = me[iel,1,1] * me[iel,2,3] - me[iel,1,3] * me[iel,2,1]
                Mm33 = me[iel,1,1] * me[iel,2,2] - me[iel,1,2] * me[iel,2,1]
                Ad11 = cf11*Mm11; Ad12 = cf21*Mm21; Ad13 = cf31*Mm31;
                Ad21 = cf12*Mm12; Ad22 = cf22*Mm22; Ad23 = cf32*Mm32;
                Ad31 = cf13*Mm13; Ad32 = cf23*Mm23; Ad33 = cf33*Mm33;
                DetM = me[iel,1,1] * Mm11 - me[iel,1,2]*Mm12 + me[iel,1,3]*Mm13
                mei[iel,1,1] = 1.0/DetM*Ad11; mei[iel,1,2] = 1.0/DetM*Ad12; mei[iel,1,3] = 1.0/DetM*Ad13;
                mei[iel,2,1] = 1.0/DetM*Ad21; mei[iel,2,2] = 1.0/DetM*Ad22; mei[iel,2,3] = 1.0/DetM*Ad23;
                mei[iel,3,1] = 1.0/DetM*Ad31; mei[iel,3,2] = 1.0/DetM*Ad32; mei[iel,3,3] = 1.0/DetM*Ad33;
            end
        end

        # Assembly of b vector
        @tturbo for iel=1:mesh.nel

            be_o2[iel,1,1] = mesh.vole[iel]*sex[iel]
            be_o2[iel,2,1] = mesh.vole[iel]*sey[iel]

            for ifac=1:mesh.nf_el

                nodei = mesh.e2f[iel,ifac]
                bc    = mesh.bc[nodei]
                dAi   = mesh.dA[iel,ifac]
                taui  = StabParam(tau,dAi,mesh.vole[iel],mesh.type,mesh.ke[iel])                              # Stabilisation parameter for the face

                for in=1:nen
                    be_o2[iel,1,in] += (bc==1) * taui * djx[iel,ifac,in]
                    be_o2[iel,2,in] += (bc==1) * taui * djy[iel,ifac,in]
                end

            end
        end
    end
    return ae, be, be_o2, ze, pe, mei, pe, rjx, rjy

end

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

function ComputeElementValues_o2(mesh, Vxh, Vyh, Pe, ae, be, be_o2, ze, rjx, rjy, mei, VxDir, VyDir, tau, o2)
    nen         = 3
    Vxe         = zeros(mesh.nel);
    Vye         = zeros(mesh.nel);
    Vxe1        = zeros(mesh.nel,nen);
    Vye1        = zeros(mesh.nel,nen);
    Vxe2        = zeros(mesh.nel,nen);
    Vye2        = zeros(mesh.nel,nen);
    Txxe        = zeros(mesh.nel);
    Tyye        = zeros(mesh.nel);
    Txye        = zeros(mesh.nel);

    @tturbo for iel=1:mesh.nel

        Vxe[iel]    =  be[iel,1]/ae[iel]
        Vye[iel]    =  be[iel,2]/ae[iel]
        for in=1:nen
            Vxe1[iel,in]  =  be_o2[iel,1,in]
            Vye1[iel,in]  =  be_o2[iel,2,in]
        end
        Txxe[iel]   =  mesh.ke[iel]/mesh.vole[iel]*ze[iel,1,1]
        Tyye[iel]   =  mesh.ke[iel]/mesh.vole[iel]*ze[iel,2,2]
        Txye[iel]   =  mesh.ke[iel]/mesh.vole[iel]*0.5*(ze[iel,1,2]+ze[iel,2,1])

        for ifac=1:mesh.nf_el

            # Face
            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel],mesh.type,mesh.ke[iel])      # Stabilisation parameter for the face

            # First order
            Vxe[iel]  += (bc!=1) *  dAi*taui*Vxh[nodei]/ae[iel]
            Vye[iel]  += (bc!=1) *  dAi*taui*Vyh[nodei]/ae[iel]

            # Second order
            for in=1:nen
                Vxe1[iel,in] += (bc!=1) * taui*rjx[iel,ifac,in]*Vxh[nodei]
                Vye1[iel,in] += (bc!=1) * taui*rjy[iel,ifac,in]*Vyh[nodei]
            end

            # Assemble
            Txxe[iel] += (bc!=1) *  mesh.ke[iel]/mesh.vole[iel]*dAi*ni_x*Vxh[nodei]
            Tyye[iel] += (bc!=1) *  mesh.ke[iel]/mesh.vole[iel]*dAi*ni_y*Vyh[nodei]
            Txye[iel] += (bc!=1) *  mesh.ke[iel]*0.5*( 1.0/mesh.vole[iel]*dAi*( ni_x*Vyh[nodei] + ni_y*Vxh[nodei] ) )
        end
        Txxe[iel] *= 2.0
        Tyye[iel] *= 2.0
        Txye[iel] *= 2.0
        # inv(me) * ue
        Vxe2[iel,1] = (o2==1) * (mei[iel,1,1] * Vxe1[iel,1] + mei[iel,1,2] * Vxe1[iel,2] + mei[iel,1,3] * Vxe1[iel,3])
        Vxe2[iel,2] = (o2==1) * (mei[iel,2,1] * Vxe1[iel,1] + mei[iel,2,2] * Vxe1[iel,2] + mei[iel,2,3] * Vxe1[iel,3])
        Vxe2[iel,3] = (o2==1) * (mei[iel,3,1] * Vxe1[iel,1] + mei[iel,3,2] * Vxe1[iel,2] + mei[iel,3,3] * Vxe1[iel,3])
        # inv(me) * ue
        Vye2[iel,1] = (o2==1) * (mei[iel,1,1] * Vye1[iel,1] + mei[iel,1,2] * Vye1[iel,2] + mei[iel,1,3] * Vye1[iel,3])
        Vye2[iel,2] = (o2==1) * (mei[iel,2,1] * Vye1[iel,1] + mei[iel,2,2] * Vye1[iel,2] + mei[iel,2,3] * Vye1[iel,3])
        Vye2[iel,3] = (o2==1) * (mei[iel,3,1] * Vye1[iel,1] + mei[iel,3,2] * Vye1[iel,2] + mei[iel,3,3] * Vye1[iel,3])
    end

    if o2==1 Vxe .= Vxe2[:,1] end # extract first values of the ue solution vector (linear function evaluated at cell centroid)
    if o2==1 Vye .= Vye2[:,1] end # extract first values of the ue solution vector (linear function evaluated at cell centroid)

    return Vxe, Vye, Txxe, Tyye, Txye
end

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

function ElementAssemblyLoop_o2(mesh, ae, be, be_o2, ze, mei, pe, rjx, rjy, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, gbar, tau, o2)
    # Assemble element matrices and rhs
    f   = zeros(mesh.nf)
    Kuu = zeros(mesh.nel, 2*mesh.nf_el, 2*mesh.nf_el)
    fu  = zeros(mesh.nel, 2*mesh.nf_el)
    Kup = zeros(mesh.nel, 2*mesh.nf_el);
    fp  = zeros(mesh.nel)

    @tturbo for iel=1:mesh.nel

        for ifac=1:mesh.nf_el

            nodei = mesh.e2f[iel,ifac]
            bci   = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel], mesh.type,mesh.ke[iel])

            for jfac=1:mesh.nf_el

                nodej = mesh.e2f[iel,jfac]
                bcj   = mesh.bc[nodej]
                dAj   = mesh.dA[iel,jfac]
                nj_x  = mesh.n_x[iel,jfac]
                nj_y  = mesh.n_y[iel,jfac]
                tauj  = StabParam(tau,dAj,mesh.vole[iel], mesh.type,mesh.ke[iel])

                # Delta
                del = 0.0 + (ifac==jfac)*1.0

                # Element matrix
                nitnj   = ni_x*nj_x + ni_y*nj_y;

                res     = 1.0/ae[iel] * dAj
                # inv(me) * be_x
                meidx1  = mei[iel,1,1] * rjx[iel,jfac,1] + mei[iel,1,2] * rjx[iel,jfac,2] + mei[iel,1,3] * rjx[iel,jfac,3]
                meidx2  = mei[iel,2,1] * rjx[iel,jfac,1] + mei[iel,2,2] * rjx[iel,jfac,2] + mei[iel,2,3] * rjx[iel,jfac,3]
                meidx3  = mei[iel,3,1] * rjx[iel,jfac,1] + mei[iel,3,2] * rjx[iel,jfac,2] + mei[iel,3,3] * rjx[iel,jfac,3]
                rx      = (pe[iel,ifac,1] * meidx1 + pe[iel,ifac,2] * meidx2 + pe[iel,ifac,3] * meidx3)
                resx    = (o2==0) * res + (o2==1) * rx
                # inv(me) * be_y
                meidy1  = mei[iel,1,1] * rjy[iel,jfac,1] + mei[iel,1,2] * rjy[iel,jfac,2] + mei[iel,1,3] * rjy[iel,jfac,3]
                meidy2  = mei[iel,2,1] * rjy[iel,jfac,1] + mei[iel,2,2] * rjy[iel,jfac,2] + mei[iel,2,3] * rjy[iel,jfac,3]
                meidy3  = mei[iel,3,1] * rjy[iel,jfac,1] + mei[iel,3,2] * rjy[iel,jfac,2] + mei[iel,3,3] * rjy[iel,jfac,3]
                ry      = (pe[iel,ifac,1] * meidy1 + pe[iel,ifac,2] * meidy2 + pe[iel,ifac,3] * meidy3)
                resy    = (o2==0) * res + (o2==1) * ry

                Kexij =-dAi * (resx * taui*tauj - mesh.ke[iel]/mesh.vole[iel]*dAj*nitnj - taui*del)
                Keyij =-dAi * (resy * taui*tauj - mesh.ke[iel]/mesh.vole[iel]*dAj*nitnj - taui*del)
                yes   = (bci!=1) & (bcj!=1)
                Kuu[iel,ifac,           jfac           ] = yes * Kexij
                Kuu[iel,ifac+mesh.nf_el,jfac+mesh.nf_el] = yes * Keyij
                # Kuu[iel,ifac,           jfac           ] = (bci!=1) * (bcj!=1) * Ke_ij
                # Kuu[iel,ifac+mesh.nf_el,jfac+mesh.nf_el] = (bci!=1) * (bcj!=1) * Ke_ij
            end
            # RHS
            Xi      = 0.0 + (bci==2)*1.0
            Ji      = 0.0 + (bci==3)*1.0
            tix     = ni_x*SxxNeu[nodei] + ni_y*SxyNeu[nodei]
            tiy     = ni_x*SyxNeu[nodei] + ni_y*SyyNeu[nodei]
            nitze_x = ni_x*ze[iel,1,1] + ni_y*ze[iel,2,1]
            nitze_y = ni_x*ze[iel,1,2] + ni_y*ze[iel,2,2]
            #-----------------------------
            ax1    = 1.0/ae[iel]*be[iel,1]
            # inv(me) * be_x
            meidx1  = mei[iel,1,1] * be_o2[iel,1,1] + mei[iel,1,2] * be_o2[iel,1,2] + mei[iel,1,3] * be_o2[iel,1,3]
            meidx2  = mei[iel,2,1] * be_o2[iel,1,1] + mei[iel,2,2] * be_o2[iel,1,2] + mei[iel,2,3] * be_o2[iel,1,3]
            meidx3  = mei[iel,3,1] * be_o2[iel,1,1] + mei[iel,3,2] * be_o2[iel,1,2] + mei[iel,3,3] * be_o2[iel,1,3]
            ax2    = (pe[iel,ifac,1] * meidx1 + pe[iel,ifac,2] * meidx2 + pe[iel,ifac,3] * meidx3)
            ax     = (o2==0) * ax1 + (o2==1) * ax2
            #-----------------------------
            ay1    = 1.0/ae[iel]*be[iel,2]
            # inv(me) * be_y
            meidy1  = mei[iel,1,1] * be_o2[iel,2,1] + mei[iel,1,2] * be_o2[iel,2,2] + mei[iel,1,3] * be_o2[iel,2,3]
            meidy2  = mei[iel,2,1] * be_o2[iel,2,1] + mei[iel,2,2] * be_o2[iel,2,2] + mei[iel,2,3] * be_o2[iel,2,3]
            meidy3  = mei[iel,3,1] * be_o2[iel,2,1] + mei[iel,3,2] * be_o2[iel,2,2] + mei[iel,3,3] * be_o2[iel,2,3]
            ay2    = (pe[iel,ifac,1] * meidy1 + pe[iel,ifac,2] * meidy2 + pe[iel,ifac,3] * meidy3)
            ay     = (o2==0) * ay1 + (o2==1) * ay2
            #-----------------------------
            feix    = (bci!=1) * -dAi * (mesh.ke[iel]/mesh.vole[iel]*nitze_x - tix*Xi - gbar[iel,ifac,1]*Ji - ax*taui)
            feiy    = (bci!=1) * -dAi * (mesh.ke[iel]/mesh.vole[iel]*nitze_y - tiy*Xi - gbar[iel,ifac,2]*Ji - ay*taui)
            # up block
            Kup[iel,ifac]                            -= (bci!=1) * dAi*ni_x;
            Kup[iel,ifac+mesh.nf_el]                 -= (bci!=1) * dAi*ni_y;
            # Dirichlet nodes - uu block
            Kuu[iel,ifac,ifac]                       += (bci==1) * 1e0
            Kuu[iel,ifac+mesh.nf_el,ifac+mesh.nf_el] += (bci==1) * 1e0
            fu[iel,ifac]                             += (bci!=1) * feix + (bci==1) * VxDir[nodei] * 1e0
            fu[iel,ifac+mesh.nf_el]                  += (bci!=1) * feiy + (bci==1) * VyDir[nodei] * 1e0
            # Dirichlet nodes - pressure RHS
            fp[iel]                                  += (bci==1) * -dAi*(VxDir[nodei]*ni_x + VyDir[nodei]*ni_y)
        end
    end
    return Kuu, fu, Kup, fp
end

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

# function CreateTripletsSparse(mesh, Kuu_v, fu_v, Kup_v)
#     # Create triplets and assemble sparse matrix fo Kuu
#     e2fu = mesh.e2f
#     e2fv = mesh.e2f .+ mesh.nf
#     e2f  = hcat(e2fu, e2fv)
#     idof = 1:mesh.nf_el*2
#     ii   = repeat(idof, 1, length(idof))'
#     ij   = repeat(idof, 1, length(idof))
#     Ki   = e2f[:,ii]
#     Kif  = e2f[:,ii[1,:]]
#     Kj   = e2f[:,ij]
#     Kuu  = sparse(Ki[:], Kj[:], Kuu_v[:], mesh.nf*2, mesh.nf*2)
#     # file = matopen(string(@__DIR__,"/results/matrix_uu.mat"), "w" )
#     # write(file, "Ki",       Ki[:] )
#     # write(file, "Kj",    Kj[:] )
#     # write(file, "Kuu",  Kuu_v[:] )
#     # write(file, "nrow",  mesh.nf*2 )
#     # write(file, "ncol",  mesh.nf*2 )
#     # close(file)
#     fu   = sparse(Kif[:], ones(size(Kif[:])), fu_v[:], mesh.nf*2, 1)
#     fu   = Array(fu)
#     droptol!(Kuu, 1e-6)
#     # Create triplets and assemble sparse matrix fo Kup
#     idof = 1:mesh.nf_el*2
#     ii   = repeat(idof, 1, mesh.nel)'
#     ij   = repeat(1:mesh.nel, 1, length(idof))
#     Ki   = e2f
#     Kj   = ij
#     Kup  = sparse(Ki[:], Kj[:], Kup_v[:], mesh.nf*2, mesh.nel  )
#     # file = matopen(string(@__DIR__,"/results/matrix_up.mat"), "w" )
#     # write(file, "Ki",       Ki[:] )
#     # write(file, "Kj",    Kj[:] )
#     # write(file, "Kup",  Kup_v[:] )
#     # write(file, "nrow",  mesh.nf*2 )
#     # write(file, "ncol",  mesh.nel )
#     # close(file)
#     return Kuu, fu, Kup
# end

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


text:
using CUDA
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity,
                MixingLength, SurfaceTKEFlux, VerticallyImplicitTimeDiscretization

parameter_guide = Dict(

     :Cᴰ   => (
          name = "Dissipation parameter (TKE equation)",
          latex = L"C^D",
          default = 2.9079,
          bounds = (0.0, 10.0)),

     :Cᴸᵇ  => (
          name = "Mixing length parameter",
          latex = L"C^{\ell}_b",
          default = 1.1612,
          bounds = (0.0, 7.0)),

     :Cᵟu => (
          name = "Ratio of mixing length to grid spacing",
          latex = L"C^{\delta}u",
          default = 0.5,
          bounds = (0.0, 2.0)),

     :Cᵟc => (
          name = "Ratio of mixing length to grid spacing",
          latex = L"C^{\delta}c",
          default = 0.5,
          bounds = (0.0, 2.0)),

     :Cᵟe => (
          name = "Ratio of mixing length to grid spacing",
          latex = L"C^{\delta}e",
          default = 0.5,
          bounds = (0.0, 2.0)),

     :Cᵂu★  => (
          name = "Mixing length parameter",
          latex = L"C^W_{u\star}",
          default = 3.6188,
          bounds = (0.0, 20.0)),

     :CʷwΔ  => (
          name = "Mixing length parameter",
          latex = L"C^ww\Delta",
          default = 1.3052,
          bounds = (0.0, 5.0)),

     :CᴷRiʷ  => (
          name = "Stability function parameter",
          latex = L"C^KRi^w",
          default = 0.7213,
          bounds = (0.0, 6.0)),

     :CᴷRiᶜ  => (
          name = "Stability function parameter",
          latex = L"C^KRi^c",
          default = 0.7588,
          bounds = (-1.5, 4.0)),

     :Cᴷu⁻  => (
          name = "Velocity diffusivity LB",
          latex = L"C^Ku^-",
          default = 0.1513,
          bounds = (0.0, 1.0)),

     :Cᴷuʳ  => (
          name = "Velocity diffusivity (UB-LB)/LB",
          latex = L"C^Ku^r",
          default = 3.8493,
          bounds = (0.0, 100.0)),

     :Cᴷc⁻  => (
          name = "Velocity diffusivity LB",
          latex = L"C^Kc^-",
          default = 0.3977,
          bounds = (0.0, 3.0)),

     :Cᴷcʳ  => (
          name = "Velocity diffusivity (UB-LB)/LB",
          latex = L"C^Kc^r",
          default = 3.4601,
          bounds = (0.0, 100.0)),

     :Cᴷe⁻  => (
          name = "Velocity diffusivity LB",
          latex = L"C^Ke^-",
          default = 0.1334,
          bounds = (0.0, 3.0)),

     :Cᴷeʳ  => (
          name = "Velocity diffusivity (UB-LB)/LB",
          latex = L"C^Ke^r",
          default = 8.1806,
          bounds = (0.0, 100.0)),

     :Cᴬu  => (
          name = "Convective adjustment velocity parameter",
          latex = L"C^A_U",
          default = 0.0057,
          bounds = (0.0, 0.3)),

     :Cᴬc  => (
          name = "Convective adjustment tracer parameter",
          latex = L"C^A_C",
          default = 0.6706,
          bounds = (0.0, 2.0)),

     :Cᴬe => (
          name = "Convective adjustment TKE parameter",
          latex = L"C^A_E",
          default = 0.2717,
          bounds = (0.0, 2.0)),
)

@free_parameters(CATKEParametersRiDependent,
                 Cᵟu, Cᵟc, Cᵟe,
                 CᴷRiʷ, CᴷRiᶜ,
                 Cᴷu⁻, Cᴷuʳ, Cᴷc⁻, Cᴷcʳ, Cᴷe⁻, Cᴷeʳ,
                 Cᴰ, Cᴸᵇ, Cᵂu★, CᵂwΔ)

@free_parameters(CATKEParametersRiIndependent,
                 Cᵟu, Cᵟc, Cᵟe,
                 Cᴷu⁻, Cᴷc⁻, Cᴷe⁻,
                 Cᴰ, Cᴸᵇ, Cᵂu★, CᵂwΔ)

@free_parameters(CATKEParametersRiDependentConvectiveAdjustment,
                 Cᵟu, Cᵟc, Cᵟe,
                 CᴷRiʷ, CᴷRiᶜ,
                 Cᴬu, Cᴬc, Cᴬe,
                 Cᴷu⁻, Cᴷuʳ, Cᴷc⁻, Cᴷcʳ, Cᴷe⁻, Cᴷeʳ,
                 Cᴰ, Cᴸᵇ, Cᵂu★, CᵂwΔ)

@free_parameters(CATKEParametersRiIndependentConvectiveAdjustment,
                 Cᵟu, Cᵟc, Cᵟe,
                 Cᴷu⁻, Cᴷc⁻, Cᴷe⁻,
                 Cᴬu, Cᴬc, Cᴬe,
                 Cᴰ, Cᴸᵇ, Cᵂu★, CᵂwΔ)

free_parameter_options = Dict(
    "CATKEParametersRiDependent" => CATKEParametersRiDependent,
#     "TKEFreeConvection" => TKEFreeConvection,
#     "TKEBCParameters" => TKEBCParameters,
)

parameter_specific_kwargs = Dict(
   CATKEParametersRiDependent => (mixing_length = MixingLength(Cᴬu=0.0, Cᴬc=0.0, Cᴬe=0.0),
                               ),

   CATKEParametersRiIndependent => (mixing_length = MixingLength(Cᴷuʳ=0.0, Cᴷcʳ=0.0, Cᴷeʳ=0.0,
                                                                Cᴬu=0.0, Cᴬc=0.0, Cᴬe=0.0),
                               ),

   CATKEParametersRiDependentConvectiveAdjustment => (mixing_length = MixingLength(),
                               ),

   CATKEParametersRiIndependentConvectiveAdjustment => (mixing_length = MixingLength(Cᴷuʳ=0.0, Cᴷcʳ=0.0, Cᴷeʳ=0.0),
                               ),
)

set_if_present!(obj, name, field) = name ∈ propertynames(obj) && setproperty!(obj, name, field)

function set_to_my_defaults!(ip::InverseProblem)

     defaults = ip.default_parameters

    for (pname, info) in parameter_guide
        set_if_present!(defaults, pname, info.default)
    end
end


text:
println("\n\n\nRunning tests ...")

using QuickShiftClustering, FactCheck

quickshift(rand(2,1000))

println("  done running tests!")


text:
@testset "Test Constructor I" begin
    mdl = Linear(7, 2)

    yterms = [[1, 4, 5], [1]]
    uterms = [[1, 2], [3]]

    dynmdl = DifferenceEquation(mdl,
                                yterms=yterms,
                                uterms=uterms)
    @test dynmdl.ny == 2
    @test dynmdl.nu == 2
    @test dynmdl.n == 5
    @test dynmdl.m == 3
    @test dynmdl.iodelay == 1
    @test_throws ArgumentError DifferenceEquation(Linear(8, 3),
                                                  yterms=yterms,
                                                  uterms=uterms)
    @test_throws ArgumentError DifferenceEquation(Linear(7, 4),
                                                  yterms=yterms,
                                                  uterms=uterms)
    @test_throws ArgumentError DifferenceEquation(mdl,
                                                  yterms=[[0, 4, 5], [1]],
                                                  uterms=uterms)
end

@testset "Test Constructor II" begin
    mdl = Linear(7, 1)

    yterms = [[1, 2, 3]]
    uterms = [[0, 4, 5, 6]]

    dynmdl = DifferenceEquation(mdl, yterms=yterms,
                                uterms=uterms)
    @test dynmdl.ny == 1
    @test dynmdl.nu == 1
    @test dynmdl.n == 3
    @test dynmdl.m == 7
    @test dynmdl.iodelay == 0
    @test_throws ArgumentError DifferenceEquation(Linear(11, 1),
                                                  yterms=yterms,
                                                  uterms=uterms)
    @test_throws ArgumentError DifferenceEquation(Linear(10, 2),
                                                  yterms=yterms,
                                                  uterms=uterms)
    @test_throws ArgumentError DifferenceEquation(mdl,
                                                  yterms=[[0, 1, 2]],
                                                  uterms=uterms)
end

@testset "Test Type Construction III" begin
    mdl = Linear(4, 1)

    uterms = [[0, 4, 5, 6]]

    dynmdl = DifferenceEquation(mdl,
                                uterms=uterms)
    @test dynmdl.ny == 1
    @test dynmdl.nu == 1
    @test dynmdl.n == 0
    @test dynmdl.m == 7
    @test dynmdl.iodelay == 0
    @test_throws ArgumentError DifferenceEquation(Linear(6, 1),
                                                  uterms=uterms)
end

@testset "Test Type Construction IV" begin
    mdl = Linear(4, 1)

    yterms = [[1, 4, 5, 6]]

    dynmdl = DifferenceEquation(mdl,
                                yterms=yterms)
    @test dynmdl.ny == 1
    @test dynmdl.nu == 0
    @test dynmdl.n == 6
    @test dynmdl.m == 0
    @test dynmdl.iodelay == 0
    @test_throws ArgumentError DifferenceEquation(Linear(5, 1),
                                                  yterms=yterms)
    @test_throws ArgumentError DifferenceEquation(Linear(4, 2),
                                                  yterms=yterms)
    @test_throws ArgumentError DifferenceEquation(mdl,
                                                  yterms=[[0, 1, 2, 6]])
end

@testset "Test Single Step Prediction" begin
    mdl = LinearSparse(13, 2, [1, 1, 1, 2, 2], [1, 6, 11, 1, 2])

    yterms = [[1, 2, 3], [1, 2]]
    uterms = [[0, 1, 2], [1], [2], [0, 1], [1]]

    dynmdl = DifferenceEquation(mdl,
                                yterms=yterms,
                                uterms=uterms)

    y0 = Float64[1 2 3;
                 4 5 6]
    u = [7  8  9;
         10 11 12;
         13 14 15;
         0 16 17;
         0 18 19]
    y = Array{Float64}(2, 1)

    Θ = Float64[1, 1, 1, 1, 2]

    @testset "Simulation" begin
        simulate!(dynmdl, u, Θ, y0, y)

        @test y == reshape([29; 7], 2, 1)

        y2 = simulate(dynmdl, u, Θ, y0)

        @test y2 == reshape([29; 7], 2, 1)
    end

    @testset "Derivatives" begin
        function cal_Θ(theta)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, theta, y0, yy)
            return vec(yy)
        end
        function cal_y0(x0)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, Θ, reshape(x0, size(y0)), yy)
            return vec(yy)
        end

        dy0_numeric = Calculus.finite_difference_jacobian(cal_y0, vec(y0))
        dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

        dΘ = Array{Float64}(2, 1, 5)
        dy0 = Array{Float64}(2, 1, 2, 3)
        simulate!(dynmdl, u, Θ, y0, y, dΘ=dΘ, dy0=dy0)

        @test vec(dΘ) ≈ vec(dΘ_numeric)
        @test vec(dy0) ≈ vec(dy0_numeric)
    end
end

@testset "Test Single Step Prediction empty u" begin
    mdl = LinearSparse(5, 2, [1, 2, 2], [1, 1, 2])

    yterms = [[1, 2, 3], [1, 2]]

    dynmdl = DifferenceEquation(mdl,
                                yterms=yterms)

    y0 = Float64[1 2 3;
                 4 5 6]
    u = Matrix{Float64}(0, 0)
    y = Array{Float64, 2}(2, 1)

    Θ = Float64[1, 1, 1]

    @testset "Simulation" begin
        simulate!(dynmdl, u, Θ, y0, y)

        @test y == reshape([3; 5], 2, 1)

        y2 = simulate(dynmdl, u, Θ, y0)

        @test y2 == reshape([3; 5], 2, 1)
    end

    @testset "Derivatives" begin
        function cal_Θ(theta)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, theta, y0, yy)
            return vec(yy)
        end
        function cal_y0(x0)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, Θ, reshape(x0, size(y0)), yy)
            return vec(yy)
        end

        dy0_numeric = Calculus.finite_difference_jacobian(cal_y0, vec(y0))
        dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

        dΘ = Array{Float64}(2, 1, 3)
        dy0 = Array{Float64}(2, 1, 2, 3)
        simulate!(dynmdl, u, Θ, y0, y, dΘ=dΘ, dy0=dy0)

        @test vec(dΘ) ≈ vec(dΘ_numeric)
        @test vec(dy0) ≈ vec(dy0_numeric)
    end
end

@testset "Test Single Step Prediction Empty y" begin
    mdl = LinearSparse(5, 2, [1, 2, 2], [1, 1, 2])

    uterms = [[1, 2, 3], [1, 2]]

    dynmdl = DifferenceEquation(mdl, uterms=uterms)

    u = [1 2 3;
         4 5 6]
    y0 = Matrix{Float64}(0, 0)
    y = Matrix{Float64}(2, 1)

    Θ = Float64[1, 1, 1]

    @testset "Simulation" begin
        simulate!(dynmdl, u, Θ, y0, y)

        @test y == reshape([3; 5], 2, 1)

        y2 = simulate(dynmdl, u, Θ, y0)

        @test y2 == reshape([3; 5], 2, 1)
    end

    @testset "Derivatives" begin
        function cal_Θ(theta)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, theta, y0, yy)
            return vec(yy)
        end
        function cal_y0(x0)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, Θ, reshape(x0, size(y0)), yy)
            return vec(yy)
        end

        dy0_numeric = Calculus.finite_difference_jacobian(cal_y0, vec(y0))
        dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

        dΘ = Array{Float64}(2, 1, 3)
        dy0 = Array{Float64}(2, 1, 0, 0)
        simulate!(dynmdl, u, Θ, y0, y, dΘ=dΘ, dy0=dy0)

        @test vec(dΘ) ≈ vec(dΘ_numeric)
        @test vec(dy0) ≈ vec(dy0_numeric)
    end
end

@testset "Test SISO Multi Step Simulation" begin
    mdl = Linear(6, 1)

    yterms = [[1, 2, 3]]
    uterms = [[1, 2, 3]]

    dynmdl = DifferenceEquation(mdl, yterms=yterms,  uterms=uterms)
    y0 = reshape([0.0; 4.0; 13.0], 1, 3)

    yref = reshape([24.8; 30.2; 27.96; 11.76;
                    -16.568; -50.384; -64.7776;
                    -57.784; -26.03872; 22.23296;
                    76.507136; 98.513024; 86.7539072;
                    96.2449152; 154.54439424; 244.89924352;
                    290.512711168], 1, 17)
    u = reshape([1; 1; 1;  1;  1;  -2;  -2;
                 -2;  -2;  -2;  3;  3;  3;
                 3;  3;  10;  10;  10;  10], 1, 19)
    y = Array{Float64, 2}(1, 17)

    Θ = [1, -0.8, 0.2, 4, 5, 6]

    @testset "Simulation" begin
        simulate!(dynmdl, u, Θ, y0, y)

        @test y ≈ yref

        y2 = simulate(dynmdl, u, Θ, y0)

        @test y2 ≈ yref
    end

    @testset "Derivatives" begin
        function cal_Θ(theta)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, theta, y0, yy)
            return vec(yy)
        end
        function cal_y0(x0)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, Θ, reshape(x0, size(y0)), yy)
            return vec(yy)
        end

        dy0_numeric = Calculus.finite_difference_jacobian(cal_y0, vec(y0))
        dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

        dΘ = Array{Float64}(1, 17, 6)
        dy0 = Array{Float64}(1, 17, 1, 3)
        simulate!(dynmdl, u, Θ, y0, y, dΘ=dΘ, dy0=dy0)

        @test vec(dΘ) ≈ vec(dΘ_numeric)
        @test vec(dy0) ≈ vec(dy0_numeric)
    end
end

@testset "Test MIMO multistep prediction" begin
    mdl = Linear(12, 2)

    yterms = [[1, 2, 3], [1, 2, 3]]
    uterms = [[1, 2, 3], [1, 2, 3]]

    dynmdl = DifferenceEquation(mdl,
                                yterms=yterms,
                                uterms=uterms)

    y0 = Float64[0 0;
                 4 8;
                 13 26]'
    yref = [24.8 49.6;
            30.2 60.4;
            27.96 63.92;
            11.76 61.52;
            -16.568 55.464;
            -50.384 34.032;
            -64.7776 8.9648;
            -57.784 -10.1680;
            -26.03872 -5.53344;
            22.23296 34.39392;
            76.507136 112.787072;
            98.513024 180.165248;
            86.7539072 216.8143744;
            96.2449152 243.2395904;
            154.54439424 260.82114048;
            244.89924352 271.59234304;
            290.512711168 231.583348736]'
    u = [1 1;
         1 1;
         1 1;
         1 1;
         1 3;
         -2 3;
         -2 3;
         -2 3;
         -2 1;
         -2 1;
         3 1;
         3 1;
         3 5;
         3 5;
         3 5;
         10 5;
         10 -2;
         10 -2;
         10 -2]'
    y = Array{Float64, 2}(2, 17)

    Θ = Float64[1, 0,
                -0.8, 0,
                0.2, 0,
                0, 1,
                0,-0.8,
                0, 0.2,
                4, 4,
                5, 5,
                6, 6,
                0, 4,
                0, 5,
                0, 6]

    @testset "Simulation" begin
        simulate!(dynmdl, u, Θ, y0, y)

        @test y ≈ yref

        y2 = simulate(dynmdl, u, Θ, y0)

        @test y2 ≈ yref
    end

    @testset "Derivatives" begin
        function cal_Θ(theta)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, theta, y0, yy)
            return vec(yy)
        end
        function cal_y0(x0)
            yy = Matrix{Float64}(size(y))
            simulate!(dynmdl, u, Θ, reshape(x0, size(y0)), yy)
            return vec(yy)
        end

        dy0_numeric = Calculus.finite_difference_jacobian(cal_y0, vec(y0))
        dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

        dΘ = Array{Float64}(2, 17, 24)
        dy0 = Array{Float64}(2, 17, 2, 3)
        simulate!(dynmdl, u, Θ, y0, y, dΘ=dΘ, dy0=dy0)

        @test vec(dΘ) ≈ vec(dΘ_numeric)
        @test vec(dy0) ≈ vec(dy0_numeric)
    end
end

@testset "Test Predict Argument Check" begin
    mdl = Linear(12, 2)

    yterms = [[1, 2, 3], [1, 2, 3]]
    uterms = [[1, 2, 3], [1, 2, 3]]

    dynmdl = DifferenceEquation(mdl,
                                yterms=yterms,
                                uterms=uterms)

    y0 = Float64[0 0;
                 4 8;
                 13 26]'
    yref = [24.8 49.6;
            30.2 60.4;
            27.96 63.92;
            11.76 61.52;
            -16.568 55.464;
            -50.384 34.032;
            -64.7776 8.9648;
            -57.784 -10.1680;
            -26.03872 -5.53344;
            22.23296 34.39392;
            76.507136 112.787072;
            98.513024 180.165248;
            86.7539072 216.8143744;
            96.2449152 243.2395904;
            154.54439424 260.82114048;
            244.89924352 271.59234304;
            290.512711168 231.583348736]'
    u = [1 1;
         1 1;
         1 1;
         1 1;
         1 3;
         -2 3;
         -2 3;
         -2 3;
         -2 1;
         -2 1;
         3 1;
         3 1;
         3 5;
         3 5;
         3 5;
         10 5;
         10 -2;
         10 -2;
         10 -2]'
    y = Array{Float64, 2}(2, 17)

    Θ = Float64[1, 0,
                -0.8, 0,
                0.2, 0,
                0, 1,
                0,-0.8,
                0, 0.2,
                4, 4,
                5, 5,
                6, 6,
                0, 4,
                0, 5,
                0, 6]

    dΘ = Array{Float64}(2, 17, 24)
    dy0 = Array{Float64}(2, 17, 2, 3)

    @test_throws ArgumentError simulate!(dynmdl, u, Θ, y0,
                                         y[1:end-1, :], dΘ=dΘ,
                                         dy0=dy0)
    @test_throws ArgumentError simulate!(dynmdl, u[1:end-1, :], Θ, y0,
                                         y, dΘ=dΘ, dy0=dy0)
    @test_throws ArgumentError simulate!(dynmdl, u[:, 1:end-1], Θ, y0,
                                         y, dΘ=dΘ, dy0=dy0)
    @test_throws ArgumentError simulate!(dynmdl, u, Θ[1:end-1], y0,
                                         y, dΘ=dΘ, dy0=dy0)
    @test_throws ArgumentError simulate!(dynmdl, u, Θ, y0[:, 1:end-1],
                                         y, dΘ=dΘ, dy0=dy0)
    @test_throws ArgumentError simulate!(dynmdl, u, Θ, y0[1:end-1, :],
                                         y, dΘ=dΘ, dy0=dy0)
    @test_throws ArgumentError simulate!(dynmdl, u, Θ, y0,
                                         y, dΘ=dΘ[:, :, 1:end-1], dy0=dy0)
    @test_throws ArgumentError simulate!(dynmdl, u, Θ, y0,
                                         y, dΘ=dΘ, dy0=dy0[:, :, :, 1:end-1])
end


text:
"""
    module MaxwellOperators
"""
module MaxwellOperators

export maxwell
export maxwell_lep
export maxwell_nep

using ..Simulations
using ..Dispersions
using Flatten
using LinearAlgebra
using NonlinearEigenproblems
using SparseArrays

const I3 = sparse(I,3,3)

struct Maxwell_LEP{N,TSIM,TΣ}
    sim::TSIM
    Σs::TΣ
    curlcurl::SparseMatrixCSC{ComplexF64,Int}
    αε::SparseMatrixCSC{Complex{Float64},Int}
    αχ::SparseMatrixCSC{Complex{Float64},Int}
    kx::Float64
    ky::Float64
    kz::Float64
    ka::Float64
    kb::Float64
    kc::Float64

    function Maxwell_LEP(
                sim::Simulation{N},
                Σs::Tuple,
                curlcurl::SparseMatrixCSC{ComplexF64,Int},
                αε::SparseMatrixCSC{ComplexF64,Int},
                αχ::SparseMatrixCSC{Complex{Float64},Int},
                kx::Number,
                ky::Number,
                kz::Number,
                ka::Number,
                kb::Number,
                kc::Number) where N

        return new{N,typeof(sim),typeof(Σs)}(sim,Σs,curlcurl,αε,αχ,kx,ky,kz,ka,kb,kc)
    end
end
function Base.show(io::IO,::Maxwell_LEP{N}) where N
    print(io,"Maxwell Linear Eigenproblem (call w/ args ")
    printstyled(io,"ω",color=:cyan)
    print(io,", [")
    printstyled(io,"ωs",color=:cyan)
    print(io,",")
    printstyled(io,"ψs",color=:cyan)
    print(io,"] -> ")
    printstyled(io,"A",color=:cyan)
    print(io,", ")
    printstyled(io,"B",color=:cyan)
    print(io,")")
end

struct Maxwell{N,TM}
    lep::TM
    Maxwell(lep::Maxwell_LEP{N}) where N = new{N,typeof(lep)}(lep)
end
function Base.show(io::IO,::Maxwell{N}) where N
    print(io,"Maxwell Operator (call w/ args ")
    printstyled(io,"ω",color=:cyan)
    print(io,", [")
    printstyled(io,"ωs",color=:cyan)
    print(io,",")
    printstyled(io,"ψs",color=:cyan)
    print(io,"])")
end

include("1D/MaxwellOperators.jl")
# include("2D/MaxwellOperators.jl")
# include("3D/MaxwellOperators.jl")

"""
    maxwell_lep(simulation; ky=0, kz=0)
"""
maxwell_lep

"""
    maxwell(simulation; ky=0, kz=0)
"""
maxwell

"""
    maxwell_nep(simulation; ky=0, kz=0, kwargs...)
"""
maxwell_nep

end #module


text:
version = v"6.1.2"

include("../common.jl")

# Build the tarballs
build_tarballs(ARGS, configure(version)...;
               preferred_gcc_version=v"6", julia_compat="~1.0, ~1.1, ~1.2, ~1.3, ~1.4, ~1.5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# function computeInertiaTensorForTwoBodies!
# computes a common mass, inertia tensor and center of mass of all object3Ds of a super object
# I_steiner = E*I*E' + m * skew(rCM)' * skew(rCM)
function computeInertiaTensorForTwoBodies!(rootSuperObj::Object3D, actualMassSuperObject::Array{Object3D,1})
  if length(actualMassSuperObject) > 0
    for i=1:length(actualMassSuperObject)
      obj = actualMassSuperObject[i]
      if rootSuperObj != obj
        m_child = obj.data.massProperties.m
        m_root  = rootSuperObj.massProperties.m
        rCM_child_old = obj.data.massProperties.rCM
        R_child = obj.R_rel
        r_child = obj.r_rel
        rCM_root = rootSuperObj.massProperties.rCM
        I_child = obj.data.massProperties.I
        I_root  = rootSuperObj.massProperties.I

        # new mass
        m_new = m_root + m_child

        # new center of mass: root to child
        rCM_child_new = r_child + R_child' * rCM_child_old

        # new center of mass of both mass
        @assert(m_new > 0.0)
        rCM_new = (rCM_child_new * m_child + rCM_root * m_root)/m_new

        # new I_child_steiner and I_root_steiner & I_new
        I_child_steiner = ModiaMath.NullRotation*I_child*ModiaMath.NullRotation' + m_child * ModiaMath.skew(rCM_child_new)' * ModiaMath.skew(rCM_child_new)
        I_root_steiner  =          I_root          + m_root  * ModiaMath.skew(rCM_root)' * ModiaMath.skew(rCM_root)
        I_sum = I_child_steiner + I_root_steiner
        I_new = I_sum - m_new * ModiaMath.skew(rCM_new)' * ModiaMath.skew(rCM_new)

        # assign mass properties to rootSuperObj
        rootSuperObj.massProperties.m = m_new
        rootSuperObj.massProperties.rCM = rCM_new
        rootSuperObj.massProperties.I = I_new

        # other way to compute inertia tensor
        #=
        a1 = rCM_new - rCM_child_new
        a2 = rCM_root - rCM_new
        IchildNew =  I_child + m_child * ModiaMath.skew(a1)' * ModiaMath.skew(a1)
        IrootNew  =  I_root  + m_root  * ModiaMath.skew(a2)' * ModiaMath.skew(a2)
        InewNew = IchildNew + IrootNew
        =#
      end
    end
  else
    rootSuperObj.massProperties = nothing
  end
  return nothing
end

# assigns mass properties to the root of a super object
# differe betweeen root objects has already mass properties or not.
function initializeMassComputation!(scene::Scene)
  if scene.initMassComp != true
    superObjs = scene.superObjs
    buffer    = scene.buffer

    for i = 1:length(superObjs)
      rootSuperObj = buffer[i]
      if dataHasMass(rootSuperObj)
        rootSuperObj.massProperties     = Solids.InternalMassProperties()
        rootSuperObj.massProperties.m   = rootSuperObj.data.massProperties.m
        rootSuperObj.massProperties.rCM = rootSuperObj.data.massProperties.rCM
        rootSuperObj.massProperties.I   = rootSuperObj.data.massProperties.I
        computeInertiaTensorForTwoBodies!(rootSuperObj, superObjs[i].superObjMass.superObj)
      else
        if length(superObjs[i].superObjMass.superObj) > 0
          rootSuperObj.massProperties = Solids.InternalMassProperties()
          computeInertiaTensorForTwoBodies!(rootSuperObj, superObjs[i].superObjMass.superObj)
    end; end; end
    scene.initMassComp = true
  end
  return nothing
end


text:
@testset "Cat Covs" begin

  using CategoricalArrays

  #Test handling categorical data
  df = DataFrame(a=[1,2], b=categorical(["low","high"]))
  levels!(df.b, ["low", "high"])

  @test isequal(Unmarked.n_levels(df.a), 0)
  @test isequal(Unmarked.n_levels(df.b), 2)
  @test isequal(Unmarked.rep_levels(df.a, 5), repeat([0], 5))
  @test isequal(Unmarked.rep_levels(df.b, 3), ["low","high","low"])

  df[!, :catvar] = categorical(["cat1", "cat1"])
  levels!(df.catvar, ["cat1","cat2","cat3"])
  dfn = Unmarked.add_levels(df)

  @test isequal(nrow(dfn), 3)
  @test isequal(dfn.catvar, levels(dfn.catvar))
  @test isequal(names(dfn), names(df))

  #Fit model with categorical
  ψ_formula = @formula(ψ~elev+forest);
  p_formula = @formula(p~precip+wind);
  β_truth = [0, -0.5, 1.2, -0.2, 0, 0.7];

  umd = simulate(Occu, ψ_formula, p_formula, (100, 5), β_truth);

  ψ_formula = @formula(ψ~elev+forest+catvar);
  umd.site_covs[!,:catvar] = categorical(rand(["high", "med", "low"], 100))
  levels!(umd.site_covs.catvar, ["low","med","high"])

  fit = occu(ψ_formula, p_formula, umd);

  @test isequal(coefnames(fit)[4:5], ["catvar: med", "catvar: high"])

  #Predict
  catvar=categorical(["low","low"])
  levels!(catvar, ["low","med","high"])
  pr_df = DataFrame(elev=[0.5, -0.3], forest=[1,-1], catvar=catvar);
  predict(occupancy(fit), pr_df, interval=true)

end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
module CommunalHelperDreamSwitchGate

using ..Ahorn, Maple
using Ahorn.CommunalHelper
using Ahorn.CommunalHelperEntityPresets: CustomDreamBlockData

const entityData = appendkwargs(CustomDreamBlockData, :(
    permanent::Bool=false,
))
@mapdefdata Entity "CommunalHelper/DreamSwitchGate" DreamSwitchGate entityData

const placements = Ahorn.PlacementDict(
    "Dream Switch Gate (Communal Helper)" => Ahorn.EntityPlacement(
        DreamSwitchGate,
        "rectangle",
        Dict{String,Any}(),
        function (entity)
            entity.data["nodes"] = [(
                Int(entity.data["x"]) + Int(entity.data["width"]),
                Int(entity.data["y"]),
            )]
        end,
    ),
)

const iconSprite = Ahorn.getSprite("objects/switchgate/icon00", "Gameplay")

Ahorn.nodeLimits(entity::DreamSwitchGate) = 1, 1

Ahorn.minimumSize(entity::DreamSwitchGate) = 16, 16
Ahorn.resizable(entity::DreamSwitchGate) = true, true

function Ahorn.selection(entity::DreamSwitchGate)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

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

    return [
        Ahorn.Rectangle(x, y, width, height),
        Ahorn.Rectangle(stopX, stopY, width, height),
    ]
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::DreamSwitchGate)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

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

    renderDreamBlock(ctx, stopX, stopY, width, height, entity.data)

    Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)

    Ahorn.drawImage(ctx, iconSprite, stopX + div(width - iconSprite.width, 2), stopY + div(height - iconSprite.height, 2))
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::DreamSwitchGate)
    x, y = Ahorn.position(entity)

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

    renderDreamBlock(ctx, x, y, width, height, entity.data)

    Ahorn.drawImage(ctx, iconSprite, x + div(width - iconSprite.width, 2), y + div(height - iconSprite.height, 2))
end

end


text:
const multiary_operators = Set()

struct MultiaryOperation{LX, LY, LZ, N, O, A, I, G, T} <: AbstractOperation{LX, LY, LZ, G, T}
    op :: O
    args :: A
    ▶ :: I
    grid :: G

    function MultiaryOperation{LX, LY, LZ}(op::O, args::A, ▶::I, grid::G) where {LX, LY, LZ, O, A, I, G}
        T = eltype(grid)
        N = length(args)
        return new{LX, LY, LZ, N, O, A, I, G, T}(op, args, ▶, grid)
    end
end

@inline Base.getindex(Π::MultiaryOperation{LX, LY, LZ, N}, i, j, k)  where {LX, LY, LZ, N} =
    Π.op(ntuple(γ -> Π.▶[γ](i, j, k, Π.grid, Π.args[γ]), Val(N))...)

#####
##### MultiaryOperation construction
#####

function _multiary_operation(L, op, args, Largs, grid) where {LX, LY, LZ}
    ▶ = Tuple(interpolation_operator(La, L) for La in Largs)
    return MultiaryOperation{L[1], L[2], L[3]}(op, Tuple(a for a in args), ▶, grid)
end

# Recompute location of multiary operation
@inline at(loc, Π::MultiaryOperation) = Π.op(loc, Tuple(at(loc, a) for a in Π.args)...)

"""Return an expression that defines an abstract `MultiaryOperator` named `op` for `AbstractField`."""
function define_multiary_operator(op)
    return quote
        function $op(Lop::Tuple,
                     a::Union{Function, Number, Oceananigans.Fields.AbstractField},
                     b::Union{Function, Number, Oceananigans.Fields.AbstractField},
                     c::Union{Function, Number, Oceananigans.Fields.AbstractField},
                     d::Union{Function, Number, Oceananigans.Fields.AbstractField}...)

            args = tuple(a, b, c, d...)
            grid = Oceananigans.AbstractOperations.validate_grid(args...)

            # Convert any functions to FunctionFields
            args = Tuple(Oceananigans.Fields.fieldify_function(Lop, a, grid) for a in args)
            Largs = Tuple(Oceananigans.Fields.location(a) for a in args)

            return Oceananigans.AbstractOperations._multiary_operation(Lop, $op, args, Largs, grid)
        end

        $op(a::Oceananigans.Fields.AbstractField,
            b::Union{Function, Oceananigans.Fields.AbstractField},
            c::Union{Function, Oceananigans.Fields.AbstractField},
            d::Union{Function, Oceananigans.Fields.AbstractField}...) = $op(Oceananigans.Fields.location(a), a, b, c, d...)
    end
end

"""
    @multiary op1 op2 op3...

Turn each multiary operator in the list `(op1, op2, op3...)`
into a multiary operator on `Oceananigans.Fields` for use in `AbstractOperations`.

Note that a multiary operator:
  * is a function with two or more arguments: for example, `+(x, y, z)` is a multiary function;
  * must be imported to be extended if part of `Base`: use `import Base: op; @multiary op`;
  * can only be called on `Oceananigans.Field`s if the "location" is noted explicitly; see example.

Example
=======

```jldoctest
julia> using Oceananigans, Oceananigans.AbstractOperations

julia> harmonic_plus(a, b, c) = 1/3 * (1/a + 1/b + 1/c)
harmonic_plus (generic function with 1 method)

julia> c, d, e = Tuple(CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) for i = 1:3);

julia> harmonic_plus(c, d, e) # before magic @multiary transformation
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
└── tree:
    * at (Center, Center, Center)
    ├── 0.3333333333333333
    └── + at (Center, Center, Center)
        ├── / at (Center, Center, Center)
        │   ├── 1
        │   └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
        ├── / at (Center, Center, Center)
        │   ├── 1
        │   └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
        └── / at (Center, Center, Center)
            ├── 1
            └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU

julia> @multiary harmonic_plus
Set{Any} with 3 elements:
  :+
  :harmonic_plus
  :*

julia> harmonic_plus(c, d, e)
MultiaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
└── tree:
    harmonic_plus at (Center, Center, Center)
    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
```
"""
macro multiary(ops...)
    expr = Expr(:block)

    for op in ops
        defexpr = define_multiary_operator(op)
        push!(expr.args, :($(esc(defexpr))))

        add_to_operator_lists = quote
            push!(Oceananigans.AbstractOperations.operators, Symbol($op))
            push!(Oceananigans.AbstractOperations.multiary_operators, Symbol($op))
        end

        push!(expr.args, :($(esc(add_to_operator_lists))))
    end

    return expr
end

#####
##### Nested computations
#####

function compute_at!(Π::MultiaryOperation, time)
    for a in Π.args
        compute_at!(a, time)
    end
    return Π
end

#####
##### GPU capabilities
#####

"Adapt `MultiaryOperation` to work on the GPU via CUDAnative and CUDAdrv."
Adapt.adapt_structure(to, multiary::MultiaryOperation{LX, LY, LZ}) where {LX, LY, LZ} =
    MultiaryOperation{LX, LY, LZ}(Adapt.adapt(to, multiary.op),
                                  Adapt.adapt(to, multiary.args),
                                  Adapt.adapt(to, multiary.▶),
                                  Adapt.adapt(to, multiary.grid))


text:
using SuperStiff.PeriodizeSC
using NPZ
using Glob

# You have to run convert_converged_green() before running this program properly (produce both selfR and greenR)

filename_to_write = "PER_periodized_zz.dat"

t = 1.0; tpp = 0.0
pwd_ = pwd()
path_ = pwd_*"/"*filename_to_write
Grid_ = 100 # Resolution of k-space that is relevant ONLY for inplane calculations
OPT_ = "PER"  # Options are CUM, PER and TR
AXIS_ = "zz" # Options are xx (a axis), yy (b axis) and zz (c axis)

macro assertion(ex, text)
    :($ex ? nothing : error("Assertion failed: ", $(text)))
end

existing_mu = Array{Float64}(); existing_beta = Array{Float64}(); existing_U = Array{Float64}(); existing_tp = Array{Float64}()
if isfile(path_)
    existing_data = readdlm(path_, header=true)
    existing_tp = existing_data[1][:,2]
    existing_U = existing_data[1][:,3]
    existing_mu = existing_data[1][:,4]
    existing_beta = existing_data[1][:,5]
    #println("length(existing_mu): ", existing_mu, "\n", "length(existing_U): ", existing_U, "\n", "length(existing_beta): ", existing_beta)
    @assertion(length(existing_tp)==length(existing_mu)==length(existing_U)==length(existing_beta), "Length of existing_* differs from each other.")
else
    println("The file containing the values of the superfluid stiffness is not existing. Will be created!")
end
existing_data = collect(zip(existing_tp,existing_U,existing_mu,existing_beta))
println("existing_data array: ", existing_data)

function grasp_mu(mu_file_path::String)
    data_mu = readdlm(mu_file_path)[:,2][end]
    return data_mu
end

function load_iwn_list(w_file_path)
    data_iwn = readdlm(w_file_path)[:,1]
    return convert(Array{Complex{Float64},1},data_iwn)
end

pattern_SE = r"selfR_aver"
pattern_GE = r"greenR_aver"
list_of_npy_SE_files = Array{String,1}()
list_of_npy_GE_files = Array{String,1}()
for (root, dirs, files) in walkdir(".")
    #println("Files in $root")
    for file in files
        if ismatch(pattern_SE,file)
            push!(list_of_npy_SE_files,joinpath(root, file)) # path to files
        elseif ismatch(pattern_GE,file)
            push!(list_of_npy_GE_files,joinpath(root, file))
        end
    end
end
#println(list_of_npy_GE_files,"\n",list_of_npy_SE_files)

@assertion(length(list_of_npy_SE_files)==length(list_of_npy_GE_files), "Length of list_of_npy_SE_files differs from that of list_of_npy_GE_files.")

regex_exp_beta_n = r"(?<=beta)(\d+\.?\d*[n]\Z)"
regex_exp_beta = r"(?<=beta)(\d+\.?\d*[^n]\Z)"
regex_exp_tp = r"(?<=tp)([-+]\d*\.\d+|\d+)"
regex_exp_U = r"(?<=U)([-+]\d*\.\d+|\d+)"
regex_exp_mu = r"(?<=m)(\d*\.\d+|\d+)"

list_of_npy_GE_files_todo  = Array{String,1}()
list_of_npy_SE_files_todo = Array{String,1}()
list_list_iwn = Array{Array{Complex{Float64},1},1}()
list_N = Array{Float64,1}()
list_params = Array{Tuple{Float64,Float64,Float64,Float64},1}()
for l in 1:length(list_of_npy_SE_files)
    beta = 0.0
    dir_name = dirname(list_of_npy_SE_files[l])
    try
        beta = parse(Float64, match(regex_exp_beta,dir_name).match)
    catch e
        println("Caught the following exception: ", e)
    end
    U = parse(Float64, match(regex_exp_U,dir_name).match)
    mu = parse(Float64, match(regex_exp_mu,dir_name).match)
    tp = parse(Float64, match(regex_exp_tp,dir_name).match)
    tuple_data = ((tp,U,mu,beta),)
    if isempty(find(existing_data.==tuple_data)) && beta != 0.0
        println("tuple_data: ",tuple_data)
        N_files = glob(dir_name*"/*N.dat") # Contains the file containing the chemical potentials
        iwn_files = glob(dir_name*"/*w.dat") # Contains the file containing the Matsubara frequencies
        NN = grasp_mu(N_files...)
        iwn_list = load_iwn_list(iwn_files...)
        push!(list_N,NN)
        push!(list_list_iwn,1.0im*iwn_list) ### Added 1.0im
        push!(list_params,(tp,U,mu,beta))
        push!(list_of_npy_GE_files_todo,list_of_npy_GE_files[l])
        push!(list_of_npy_SE_files_todo,list_of_npy_SE_files[l])
    elseif beta == 0.0
        println("beta is set to: ",0.0)
    else
        println("This folder's data has already been treated!!")
        println("Set already done at index: ",find(existing_data.==tuple_data)...); println()
    end
end

println("list_params: ",list_params)

list_Kintegral_stiff = Array{Tuple{Tuple{Float64,Float64,Float64,Float64},Array{Complex{Float64},2}},1}()
list_Kintegral_DOS = Array{Tuple{Tuple{Float64,Float64,Float64,Float64},Array{Complex{Float64},2}},1}()
list_mu = Array{Float64,1}(); list_U = Array{Float64,1}(); list_beta = Array{Float64,1}(); list_gap = Array{Float64,1}()
list_tp = Array{Float64,1}(); list_N = Array{Float64,1}()
@time for l in 1:1:length(list_of_npy_GE_files_todo)
    dir_name_l = dirname(list_of_npy_GE_files_todo[l])
    SEvec = npzread(list_of_npy_SE_files_todo[l])
    println("SEvec: ",size(SEvec),"\n")
    GEvec = npzread(list_of_npy_GE_files_todo[l])
    println("GEvec: ",size(GEvec),"\n")
    len_iwn_G = size(GEvec)[1]; len_iwn_SE = size(SEvec)[1]
    len_iwn_G != len_iwn_SE && throw(ErrorException("Must have same number of Matsubara frequencies for self-energy and Green's function!"))
    N_data = readdlm(dir_name_l*"/"*"N.dat")
    list_green_npy = glob(dir_name_l*"/greenR*.npy")
    list_green_npy = Array{String,1}([dir_name for dir_name in list_green_npy if !contains(dir_name,"greenR_aver")])
    println("list_green_npy: ",list_green_npy)
    len_N_aver = length(list_green_npy) #Important to know on how many N values one averages (comes down to number of Green's function used to average)
    N_aver = 0.0
    try
        N_aver = sum(N_data[:,2][end-len_N_aver:end])/(len_N_aver+1)
        println("N_aver: ",N_aver,"\t data_N: ",N_data[:,2][end-len_N_aver:end],"\t",len_N_aver,"\n")
    catch e
        if isa(e,BoundsError)
            println("Caught BoundsError! Not enough iterations to average over.")
        else
            println("Unidentified error occured at runtime.")
        end
    end
    params_tuple = list_params[l]
    tp = params_tuple[1]; beta = params_tuple[4]; mu = params_tuple[3]
    anomalous_GE = abs((1/beta)*sum(GEvec[:,1,6]))
    println("anomalous GE: ", anomalous_GE)
    push!(list_gap, anomalous_GE)
    k_array_prebuild = Array{Array{Float64,2},1}()
    tktilde_prebuild = Array{Array{Complex{Float64},2},2}(0,0)
    modelvec = PeriodizeSC.ModelVector(t,tp,tpp,mu,list_list_iwn[l],SEvec)
    cal = Matrix{Complex{Float64}}(0,0)
    if AXIS_ == "zz"
        println("zz")
        if OPT_ == "PER"
            println("PER")
            cal = PeriodizeSC.calcintegral_BZ(modelvec, PeriodizeSC.make_stiffness_kintegrand_SC, len_sEvec_c=len_iwn_G)  ## cal variable exits an array whose size is (length(iw_n),2). The second column holds the superfluid stiffness data.
            println("cal: ", cal)
        elseif OPT_ == "CUM"
            println("CUM")
            cal = PeriodizeSC.calcintegral_BZ(modelvec, PeriodizeSC.make_stiffness_kintegrand_cum_SC, len_sEvec_c=len_iwn_G) ## Idem cal
            println("cal: ", cal)
        elseif OPT_ == "TR"
            println("TR")
            cal = PeriodizeSC.calcintegral_BZ(modelvec, PeriodizeSC.make_stiffness_trace_G_kintegrand, len_sEvec_c=len_iwn_G) ## Idem cal
        end
    elseif AXIS_ == "yy"
        println("yy")
        fct_array = [PeriodizeSC.DyDyEpsilonbark]
        for f in fct_array
            push!(k_array_prebuild, PeriodizeSC.k_grid(PeriodizeSC.Model(modelvec,1),Grid_,f))
        end
        tktilde_prebuild = PeriodizeSC.k_grid(PeriodizeSC.Model(modelvec,1),Grid_,PeriodizeSC.tktilde)
        if OPT_ == "PER"
            println("PER")
            cal = PeriodizeSC.stiffness_NOCOEX_Per_Cum_ab_k_grid(modelvec,k_array_prebuild,tktilde_prebuild,Grid_,0,0,len_sEvec_c=len_iwn_G)
        elseif OPT_ == "CUM"
            println("CUM")
            cal = PeriodizeSC.stiffness_NOCOEX_Per_Cum_ab_k_grid(modelvec,k_array_prebuild,tktilde_prebuild,Grid_,1,0,len_sEvec_c=len_iwn_G)
        end
    elseif AXIS_ == "xx"
        println("xx")
        fct_array = [PeriodizeSC.DxDxEpsilonbark]
        for f in fct_array
            push!(k_array_prebuild, PeriodizeSC.k_grid(PeriodizeSC.Model(modelvec,1),Grid_,f))
        end
        tktilde_prebuild = PeriodizeSC.k_grid(PeriodizeSC.Model(modelvec,1),Grid_,PeriodizeSC.tktilde)
        if OPT_ == "PER"
            println("PER")
            cal = PeriodizeSC.stiffness_NOCOEX_Per_Cum_ab_k_grid(modelvec,k_array_prebuild,tktilde_prebuild,Grid_,0,0,len_sEvec_c=len_iwn_G)
        elseif OPT_ == "CUM"
            println("CUM")
            cal = PeriodizeSC.stiffness_NOCOEX_Per_Cum_ab_k_grid(modelvec,k_array_prebuild,tktilde_prebuild,Grid_,1,0,len_sEvec_c=len_iwn_G)
        end
    end
    push!(list_Kintegral_stiff,(params_tuple,cal))
    push!(list_mu,mu)
    push!(list_U,params_tuple[2])
    push!(list_beta,beta)
    push!(list_tp,tp)
    push!(list_N,N_aver)
end
println("list_mu: ", list_mu)
stiffness = Array{Float64,1}()
for j in 1:size(list_Kintegral_stiff)[1]
    stiffness_tmp = 0.0
    println(list_Kintegral_stiff[j][1])
    for ii in 1:size(list_Kintegral_stiff[j][2])[1]
        stiffness_tmp+=list_Kintegral_stiff[j][2][ii, 2]
    end
    beta = list_Kintegral_stiff[j][1][4]
    println("beta: ", beta)
    push!(stiffness,(2.0/beta)*real(stiffness_tmp)) #Factor two to sum over all the imaginary axis
    println((2.0/beta)*imag(stiffness_tmp))
end

@assertion(length(list_mu)==length(stiffness), "Length of list_mu differs from that of stiffness.")

params_stiff = hcat(list_N,list_tp,list_U,list_mu,list_beta,stiffness,list_gap)
map!(x->round.(x,8),params_stiff,params_stiff)
println("Stiffness: ", stiffness)
if isfile(path_)
    open(path_,"a") do io
        writedlm(io, params_stiff, "\t\t")
    end
else
    f = open(filename_to_write, "w")
    write(f, "n"*"\t\t"*"tp"*"\t\t"*"U"*"\t\t"*"mu"*"\t\t"*"beta"*"\t\t"*"Stiff"*"\t\t"*"SC_Gap"*"\n")
    writedlm(f, params_stiff, "\t\t")
    close(f)
end

max_stars_repo_name: oliviersimard/SuperStiff

max_issues_repo_name: oliviersimard/SuperStiff

max_forks_repo_name: oliviersimard/SuperStiff


text:
@doc raw"""
    TranslationAction(
        M::Manifold,
        Rn::TranslationGroup,
        AD::ActionDirection = LeftAction(),
    )

Space of actions of the [`TranslationGroup`](@ref) $\mathrm{T}(n)$ on a Euclidean-like
manifold `M`.

The left and right actions are equivalent.
"""
struct TranslationAction{TM<:Manifold,TRn<:TranslationGroup,TAD<:ActionDirection} <:
       AbstractGroupAction{TAD}
    manifold::TM
    Rn::TRn
end

function TranslationAction(
    M::Manifold,
    Rn::TranslationGroup,
    ::TAD = LeftAction(),
) where {TAD<:ActionDirection}
    return TranslationAction{typeof(M),typeof(Rn),TAD}(M, Rn)
end

function show(io::IO, A::TranslationAction)
    print(io, "TranslationAction($(A.manifold), $(A.Rn), $(direction(A)))")
end

base_group(A::TranslationAction) = A.Rn

g_manifold(A::TranslationAction) = A.manifold

function switch_direction(A::TranslationAction{TM,TRN,TAD}) where {TM,TRN,TAD}
    return TranslationAction(A.manifold, A.Rn, switch_direction(TAD()))
end

apply(A::TranslationAction, a, p) = p + a

apply!(A::TranslationAction{M,G}, q, a, p) where {M,G} = (q .= p .+ a)
apply!(A::TranslationAction{M,G}, q, e::Identity{G}, p) where {M,G} = copyto!(q, p)

inverse_apply(A::TranslationAction, a, p) = p - a

inverse_apply!(A::TranslationAction{M,G}, q, a, p) where {M,G} = (q .= p .- a)
inverse_apply!(A::TranslationAction{M,G}, q, e::Identity{G}, p) where {M,G} = copyto!(q, p)

apply_diff(A::TranslationAction, a, p, X) = X

apply_diff!(A::TranslationAction, Y, a, p, X) = copyto!(Y, X)

inverse_apply_diff(A::TranslationAction, a, p, X) = X

inverse_apply_diff!(A::TranslationAction, Y, a, p, X) = copyto!(Y, X)


text:
function  push!(Bop :: BFGSOperator{T, N, F1, F2, F3},
                sk:: Vector{T},
                yk:: Vector{T}) where {T, N, F1, F2, F3}

    scaling = Bop.data.scaling
    nsk2 =sk'*sk
    denom = yk'*sk
    B = Bop.data.M
    if (denom > 1.0e-20) # && (mod(iter,n)!=0)
        #self-scaled version aka Luenberger
        By = B*yk
        ytBy = (yk'*By)
        γ = 1.0
        if scaling
            γ = denom / (ytBy)
        end

        # math formula, involves matrix multiplications, thus O(n³) complexity
        M = I - sk*yk'/denom
        B1 = γ*M*B*M' + sk*sk'/denom
        B = 0.5*(B1+B1')  # make sure B is symmetric
    else
        @warn "No update, B=I, denom = ", denom
        B = I
    end

    Bop.data.M .= B

    return Bop
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
using Augmentor, ISICArchive, Images, Colors
import ImageTransformations
using Reel, PaddedViews, OffsetArrays

Reel.set_output_type("gif")

_mkpath(fname) = joinpath(dirname(@__FILE__()), fname)
srand(123)

det_ops = (
    :(FlipX()),
    :(FlipY()),
    :(Rotate90()),
    :(Rotate270()),
    :(Rotate180()),
    :(Rotate(15)),
    :(ShearX(10)),
    :(ShearY(10)),
    :(Scale(0.9,0.5)),
    :(Scale(1.2)),
    :(Crop(70:140,25:155)),
    :(CropSize(45, 225)),
    :(Zoom(1.2)),
    :(Resize(100,150)),
    "Crop" => :((Rotate(45),Crop(1:210,1:280))),
    "CropNative" => :((Rotate(45),CropNative(1:210,1:280))),
    :(CropRatio(1)),
    :(ConvertEltype(GrayA)),
)
prob_ops = (
    10 => :(Rotate(-10:10)),
    10 => :(ShearX(-10:10)),
    10 => :(ShearY(-10:10)),
    10 => :(Zoom(0.9:0.05:1.3)),
    10 => :(Scale(0.9:0.05:1.3)),
    10 => :(ElasticDistortion(15,15,0.1)),
    10 => :(ElasticDistortion(10,10,0.2,4,3,true)),
    10 => :(RCropRatio(1)),
)

pattern = load(_mkpath("testpattern_small.png"))
pattern_noalpha = load(_mkpath("testpattern_small_noalpha.png"))

tables = Dict{String,Vector{String}}()

log = Dict{String,Int}()
function makename(symb)
    name = string(symb)
    id = if !haskey(log, name)
        get!(log, name, 1)
    else
        log[name] = log[name]+1
    end
    string(name, id > 1 ? id : "")
end

imgpath = " .. image:: https://raw.githubusercontent.com/JuliaML/FileStorage/master/Augmentor/"
inputpath = string(imgpath, "testpattern_small.png ")

extractname(pair::Pair) = string(pair.first)
extractname(expr::Expr) = string(expr.args[1])
extractexpr(pair::Pair) = pair.second
extractexpr(expr::Expr) = expr

# --------------------------------------------------------------------
# Create pngs

srand(1335)

for op_full in det_ops
    op_expr = extractexpr(op_full)
    descr = string(op_expr)
    op = eval(op_expr)
    name_raw = extractname(op_full)
    name = makename(name_raw)
    if !haskey(tables, name_raw)
        tables[name_raw] = Vector{String}()
    end
    fname = string("operations/", name, ".png")

    outpath = string(imgpath, fname, " ")
    maxlen = max(length(inputpath), length(outpath))
    outpath = rpad(outpath, maxlen, " ")
    inpath = rpad(inputpath, maxlen, " ")
    border = string("+", rpad("", maxlen, "-"), "+", rpad("", maxlen, "-"), "+")
    header = string("|", rpad(" Input", maxlen, " "), "|", rpad(" Output for ``$descr``", maxlen, " "), "|")
    fatborder = string("+", rpad("", maxlen, "="), "+", rpad("", maxlen, "="), "+")
    line = string("|", inpath, "|", outpath, "|")
    table = join((border,header,fatborder,line,border), "\n")
    push!(tables[name_raw], table)

    out = augment(pattern, op)
    save(_mkpath(fname), out)
end

# --------------------------------------------------------------------
# Create gifs

log = Dict{String,Int}()

function drawborder!(img, col)
    img[1:end,   1] .= fill(col, size(img,1))
    img[1:end, end] .= fill(col, size(img,1))
    img[1,   1:end] .= fill(col, size(img,2))
    img[end, 1:end] .= fill(col, size(img,2))
    img
end

centered(img) = OffsetArray(img, convert(Tuple, 1 .- round.(Int, ImageTransformations.center(img))))

srand(1335)

for (nframe, op_full) in prob_ops
    op_expr = extractexpr(op_full)
    descr = string(op_expr)
    op = eval(op_expr)
    name_raw = extractname(op_full)
    name = makename(name_raw)
    if !haskey(tables, name_raw)
        tables[name_raw] = Vector{String}()
    end
    fname = string("operations/", name, ".gif")

    outpath = string(imgpath, fname, " ")
    maxlen = max(length(inputpath), length(outpath))
    outpath = rpad(outpath, maxlen, " ")
    inpath = rpad(inputpath, maxlen, " ")
    border = string("+", rpad("", maxlen, "-"), "+", rpad("", maxlen, "-"), "+")
    header = string("|", rpad(" Input", maxlen, " "), "|", rpad(" Sampled outputs for ``$descr``", maxlen, " "), "|")
    fatborder = string("+", rpad("", maxlen, "="), "+", rpad("", maxlen, "="), "+")
    line = string("|", inpath, "|", outpath, "|")
    table = join((border,header,fatborder,line,border), "\n")
    push!(tables[name_raw], table)

    raw_imgs = [centered(drawborder!(augment(pattern_noalpha, op), colorant"pink")) for i in 1:nframe]
    imgs = map(parent, map(copy, [paddedviews(colorant"#F3F6F6", raw_imgs...)...]))
    insert!(imgs, 1, first(imgs)) # otherwise loop isn't smooth

    film = roll(imgs, fps = 2)
    write(_mkpath(fname), film)
end

# --------------------------------------------------------------------
# Autogenerate .rst tables in correct order.

order = (
    "FlipX", "FlipY",
    "Rotate90", "Rotate180", "Rotate270", "Rotate",
    "ShearX",
    "ShearY",
    "Scale", "Zoom",
    "ElasticDistortion",
    "Crop", "CropNative",
    "CropSize", "CropRatio",
    "RCropRatio",
    "Resize",
    "ConvertEltype",
)
@assert length(unique(order)) == length(keys(tables))

open(_mkpath("tables.txt"), "w") do io
    for name in order
        lines = tables[name]
        write(io, name, "\n", repeat("*", length(name)+5), "\n\n")
        for line in lines
            write(io, line, "\n\n")
        end
        write(io, "\n\n")
    end
end

max_stars_repo_name: JuliaTagBot/FileStorage

max_issues_repo_name: JuliaTagBot/FileStorage

max_forks_repo_name: JuliaTagBot/FileStorage


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: qldb_session

using AWS.Compat
using AWS.UUIDs
"""
    SendCommand()

Sends a command to an Amazon QLDB ledger.  Instead of interacting directly with this API, we recommend that you use the Amazon QLDB Driver or the QLDB Shell to execute data transactions on a ledger.   If you are working with an AWS SDK, use the QLDB Driver. The driver provides a high-level abstraction layer above this qldbsession data plane and manages SendCommand API calls for you. For information and a list of supported programming languages, see Getting started with the driver in the Amazon QLDB Developer Guide.   If you are working with the AWS Command Line Interface (AWS CLI), use the QLDB Shell. The shell is a command line interface that uses the QLDB Driver to interact with a ledger. For information, see Accessing Amazon QLDB using the QLDB Shell.

# Optional Parameters
- `AbortTransaction`: Command to abort the current transaction.
- `CommitTransaction`: Command to commit the specified transaction.
- `EndSession`: Command to end the current session.
- `ExecuteStatement`: Command to execute a statement in the specified transaction.
- `FetchPage`: Command to fetch a page.
- `SessionToken`: Specifies the session token for the current command. A session token is constant throughout the life of the session. To obtain a session token, run the StartSession command. This SessionToken is required for every subsequent command that is issued during the current session.
- `StartSession`: Command to start a new session. A session token is obtained as part of the response.
- `StartTransaction`: Command to start a new transaction.
"""

send_command(; aws_config::AWSConfig=global_aws_config()) = qldb_session("SendCommand"; aws_config=aws_config)
send_command(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb_session("SendCommand", args; aws_config=aws_config)


text:
# Autogenerated wrapper script for OpenSpecFun_jll for x86_64-apple-darwin-libgfortran3
export libopenspecfun

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("OpenSpecFun")
JLLWrappers.@declare_library_product(libopenspecfun, "@rpath/libopenspecfun.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libopenspecfun,
        "lib/libopenspecfun.1.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()