text:
#### SubIndex of Index. Used by SubDataFrame, DataFrameRow, and DataFrameRows

struct SubIndex{I<:AbstractIndex,S<:AbstractVector{Int},T<:AbstractVector{Int}} <: AbstractIndex{}
    parent::I
    cols::S # columns from idx selected in SubIndex
    remap::T # reverse mapping from cols to their position in the SubIndex
end

SubIndex(parent::AbstractIndex, ::Colon) = parent

Base.copy(x::SubIndex) = Index(_names(x))

@inline parentcols(ind::SubIndex) = ind.cols

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Union{Integer,AbstractVector{<:Integer}}) =
    ind.cols[idx]
Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Bool) =
    throw(ArgumentError("column indexing with Bool is not allowed"))

Base.@propagate_inbounds function parentcols(ind::SubIndex, idx::Symbol)
    parentcol = ind.parent[idx]
    @boundscheck begin
        remap = ind.remap
        remap[parentcol] == 0 && throw(ArgumentError("$idx not found"))
    end
    return parentcol
end

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::AbstractVector{Symbol}) =
    [parentcols(ind, i) for i in idx]

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Regex) =
    [parentcols(ind, i) for i in _names(ind) if occursin(idx, String(i))]

Base.@propagate_inbounds parentcols(ind::SubIndex, ::Colon) = ind.cols

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Not) = parentcols(ind, ind[idx])

function SubIndex(parent::AbstractIndex, cols::AbstractUnitRange{Int})
    l = last(cols)
    f = first(cols)
    if !checkindex(Bool, Base.OneTo(length(parent)), cols)
        throw(BoundsError("invalid columns $cols selected"))
    end
    remap = (1:l) .- f .+ 1
    SubIndex(parent, cols, remap)
end

function SubIndex(parent::AbstractIndex, cols::AbstractVector{Int})
    ncols = length(parent)
    remap = zeros(Int, ncols)
    for (i, col) in enumerate(cols)
        if !(1 <= col <= ncols)
            throw(BoundsError("column index must be greater than zero " *
                              "and not larger than number columns in the parent"))
        end
        if remap[col] != 0
            throw(ArgumentError("duplicate selected column detected"))
        end
        remap[col] = i
    end
    SubIndex(parent, cols, remap)
end

@inline SubIndex(parent::AbstractIndex, cols::ColumnIndex) =
    throw(ArgumentError("cols argument must be a vector (got $cols)"))

Base.@propagate_inbounds SubIndex(parent::AbstractIndex, cols) =
    SubIndex(parent, parent[cols])

Base.length(x::SubIndex) = length(x.cols)
Base.names(x::SubIndex) = copy(_names(x))
_names(x::SubIndex) = view(_names(x.parent), x.cols)

function Base.haskey(x::SubIndex, key::Symbol)
    haskey(x.parent, key) || return false
    pos = x.parent[key]
    remap = x.remap
    checkbounds(Bool, remap, pos) || return false
    remap[pos] > 0
end

Base.haskey(x::SubIndex, key::Integer) = 1 <= key <= length(x)
Base.haskey(x::SubIndex, key::Bool) =
    throw(ArgumentError("invalid key: $key of type Bool"))
Base.keys(x::SubIndex) = names(x)

function Base.getindex(x::SubIndex, idx::Symbol)
    remap = x.remap
    remap[x.parent[idx]]
end

function Base.getindex(x::SubIndex, idx::AbstractVector{Symbol})
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    [x[i] for i in idx]
end## SubIndex of Index. Used by SubDataFrame, DataFrameRow, and DataFrameRows

struct SubIndex{I<:AbstractIndex,S<:AbstractVector{Int},T<:AbstractVector{Int}} <: AbstractIndex
    parent::I
    cols::S # columns from idx selected in SubIndex
    remap::T # reverse mapping from cols to their position in the SubIndex
end

SubIndex(parent::AbstractIndex, ::Colon) = parent

Base.copy(x::SubIndex) = Index(_names(x))

@inline parentcols(ind::SubIndex) = ind.cols

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Union{Integer,AbstractVector{<:Integer}}) =
    ind.cols[idx]
Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Bool) =
    throw(ArgumentError("column indexing with Bool is not allowed"))

Base.@propagate_inbounds function parentcols(ind::SubIndex, idx::Symbol)
    parentcol = ind.parent[idx]
    @boundscheck begin
        remap = ind.remap
        remap[parentcol] == 0 && throw(ArgumentError("$idx not found"))
    end
    return parentcol
end

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::AbstractVector{Symbol}) =
    [parentcols(ind, i) for i in idx]

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Regex) =
    [parentcols(ind, i) for i in _names(ind) if occursin(idx, String(i))]

Base.@propagate_inbounds parentcols(ind::SubIndex, ::Colon) = ind.cols

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Not) = parentcols(ind, ind[idx])

function SubIndex(parent::AbstractIndex, cols::AbstractUnitRange{Int})
    l = last(cols)
    f = first(cols)
    if !checkindex(Bool, Base.OneTo(length(parent)), cols)
        throw(BoundsError("invalid columns $cols selected"))
    end
    remap = (1:l) .- f .+ 1
    SubIndex(parent, cols, remap)
end

function SubIndex(parent::AbstractIndex, cols::AbstractVector{Int})
    ncols = length(parent)
    remap = zeros(Int, ncols)
    for (i, col) in enumerate(cols)
        if !(1 <= col <= ncols)
            throw(BoundsError("column index must be greater than zero " *
                              "and not larger than number columns in the parent"))
        end
        if remap[col] != 0
            throw(ArgumentError("duplicate selected column detected"))
        end
        remap[col] = i
    end
    SubIndex(parent, cols, remap)
end

@inline SubIndex(parent::AbstractIndex, cols::ColumnIndex) =
    throw(ArgumentError("cols argument must be a vector (got $cols)"))

Base.@propagate_inbounds SubIndex(parent::AbstractIndex, cols) =
    SubIndex(parent, parent[cols])

Base.length(x::SubIndex) = length(x.cols)
Base.names(x::SubIndex) = copy(_names(x))
_names(x::SubIndex) = view(_names(x.parent), x.cols)

function Base.haskey(x::SubIndex, key::Symbol)
    haskey(x.parent, key) || return false
    pos = x.parent[key]
    remap = x.remap
    checkbounds(Bool, remap, pos) || return false
    remap[pos] > 0
end

Base.haskey(x::SubIndex, key::Integer) = 1 <= key <= length(x)
Base.haskey(x::SubIndex, key::Bool) =
    throw(ArgumentError("invalid key: $key of type Bool"))
Base.keys(x::SubIndex) = names(x)

function Base.getindex(x::SubIndex, idx::Symbol)
    remap = x.remap
    remap[x.parent[idx]]
end

function Base.getindex(x::SubIndex, idx::AbstractVector{Symbol})
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    [x[i] for i in idx]
end


text:
import Observables
import Makie

q1 = Biquaternion(Quaternion(rand(4)), ℝ³(rand(3)))
scene = Makie.Scene()
len = 1 + rand()
width = 5
color = [:red, :green, :blue]
transparency = false
triad = Triad(q1,
              scene,
              length = len,
              width = width,
              color = color)
q2 = q1 * Biquaternion(Quaternion(rand(4)), ℝ³(rand(3)))
value1 = Observables.to_value(triad.observable)
update(triad, q2)
value2 = Observables.to_value(triad.observable)

@test isapprox(triad.q, q2)
@test isapprox(value1, value2) == false

ext: jl

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
#=
 Train a simple network to predict formation energy per atom (downloaded from Materials Project).
=#
using CSV, DataFrames
using Random, Statistics
using Flux
using Flux: @epochs
using ChemistryFeaturization
using AtomicGraphNets

function train_formation_energy(;
    num_pts = 100,
    num_epochs = 5,
    data_dir = joinpath(@__DIR__, "data"),
    verbose = true,
)
    println("Setting things up...")

    # data-related options
    train_frac = 0.8 # what fraction for training?
    num_train = Int32(round(train_frac * num_pts))
    num_test = num_pts - num_train
    prop = "formation_energy_per_atom"
    id = "task_id" # field by which to label each input material

    # set up the featurization
    featurization = GraphNodeFeaturization([
        "Group",
        "Row",
        "Block",
        "Atomic mass",
        "Atomic radius",
        "X",
    ])
    num_features = output_shape(featurization)

    # model hyperparameters – keeping it pretty simple for now
    num_conv = 3 # how many convolutional layers?
    crys_fea_len = 32 # length of crystal feature vector after pooling (keep node dimension constant for now)
    num_hidden_layers = 1 # how many fully-connected layers after convolution and pooling?
    opt = ADAM(0.001) # optimizer

    # dataset...first, read in outputs
    info = CSV.read(string(data_dir, prop, ".csv"), DataFrame)
    y = Array(Float32.(info[!, Symbol(prop)]))

    # shuffle data and pick out subset
    indices = shuffle(1:size(info, 1))[1:num_pts]
    info = info[indices, :]
    output = y[indices]

    # next, make and featurize graphs
    if verbose
        println("Building graphs and feature vectors from structures...")
    end
    inputs = FeaturizedAtoms[]

    for r in eachrow(info)
        cifpath = string(data_dir, prop, "_cifs/", r[Symbol(id)], ".cif")
        gr = AtomGraph(cifpath)
        input = featurize(gr, featurization)
        push!(inputs, input)
    end

    # pick out train/test sets
    if verbose
        println("Dividing into train/test sets...")
    end
    train_output = output[1:num_train]
    test_output = output[num_train+1:end]
    train_input = inputs[1:num_train]
    test_input = inputs[num_train+1:end]
    train_data = zip(train_input, train_output)

    # build the model
    if verbose
        println("Building the network...")
    end
    model = build_CGCNN(
        num_features,
        num_conv = num_conv,
        atom_conv_feature_length = crys_fea_len,
        pooled_feature_length = (Int(crys_fea_len / 2)),
        num_hidden_layers = 1,
    )

    # define loss function and a callback to monitor progress
    loss(x, y) = Flux.Losses.mse(model(x), y)
    evalcb_verbose() = @show(mean(loss.(test_input, test_output)))
    evalcb_quiet() = return nothing
    evalcb = verbose ? evalcb_verbose : evalcb_quiet
    evalcb()

    # train
    if verbose
        println("Training!")
    end
    @epochs num_epochs Flux.train!(
        loss,
        Flux.params(model),
        train_data,
        opt,
        cb = Flux.throttle(evalcb, 5),
    )

    return model
end


text:
using Test

@testset "LookingGlass.jl" begin
    include("LookingGlass.jl")
end


text:
##
using MRIgeneralizedBloch
using DifferentialEquations
using StaticArrays
using Test

## choose random parameters
# this test should only be passed for small α and long TRF
α = rand() * π/2
TRF = 500e-6 + 1000e-6 * rand()
ω1 = α / TRF
B1 = 0.7 + 0.6 * rand()
ω0 = 0 # Graham's spectral model is only implemented for on-resonance pulses
m0s = 0.2 * rand()
m0f = 1 - m0s
R1f = 0.5 + 0.2 * rand()
R1s = 2 + rand()
R2f = 1 / (40e-3 + 100e-3 * rand())
T2s = 5e-6 + 10e-6 * rand()
R2s = 1 / T2s
Rx = 100 * rand()

ms = m0s * rand()
mf = (1 - m0s) * rand()
ϑ = rand() * π / 2
φ = rand() * 2π

## Solve Graham's solution
p = (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s)
alg = Tsit5()
m0 = [mf * sin(ϑ) * cos(φ), mf * sin(ϑ) * sin(φ), mf * cos(ϑ), ms, 1]

sol = solve(
    ODEProblem(MRIgeneralizedBloch.apply_hamiltonian_graham_superlorentzian!, m0, (0.0, TRF), p),
    alg,
)

u_Graham = sol[end]

## Solve generalized Bloch-McConnell with super-Lorentzian lineshape
mfun(p, t; idxs = nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5)

g = interpolate_greens_function(greens_superlorentzian, 0, TRF/T2s)
p = (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g)
alg = MethodOfSteps(DP8())
sol = solve(
    DDEProblem(MRIgeneralizedBloch.apply_hamiltonian_gbloch!, m0, mfun, (0.0, TRF), p),
    alg,
)

u_gBloch = sol[end]

## Test!
max_error = 1e-2
@test u_gBloch[1] ≈ u_Graham[1] atol = max_error
@test u_gBloch[2] ≈ u_Graham[2] atol = max_error
@test u_gBloch[3] ≈ u_Graham[3] atol = max_error
@test u_gBloch[4] ≈ u_Graham[4] atol = max_error
@test u_gBloch[5] ≈ u_Graham[5] atol = max_error


text:
fast_scalar_indexing(::Type{<:CUDA.CuArray}) = false
@inline allowed_getindex(x::CUDA.CuArray,i...) = CUDA.@allowscalar(x[i...])
@inline allowed_setindex!(x::CUDA.CuArray,v,i...) = (CUDA.@allowscalar(x[i...] = v))

function Base.setindex(x::CUDA.CuArray,v,i::Int)
  _x = copy(x)
  allowed_setindex!(_x,v,i)
  _x
end

function restructure(x::CUDA.CuArray,y)
  reshape(Adapt.adapt(parameterless_type(x),y),size(x)...)
end


text:
using HDF5
using Test
using Pkg

println("HDF5 version ", HDF5.h5_get_libversion())

include("plain.jl")
include("compound.jl")
include("custom.jl")
include("reference.jl")
include("hyperslab.jl")
include("readremote.jl")
include("extend_test.jl")
include("gc.jl")
include("external.jl")
include("swmr.jl")
include("mmap.jl")
include("properties.jl")

try
    using MPI
    # basic MPI tests, for actual parallel tests we need to run in MPI mode
    include("mpio.jl")
catch
end


text:
# Aqui tenho o PSO implementado ocnforme equação
# Eq(1): v = ω { v + ϕp . rp . (p -x) + ϕg . rg . (g - x)}
# Eq(2): x = x +v
# podendo utilizar topologia local ou global

function pso(func::Function, particles::Particles, constraints, args, kwargs,
    ψ , ϕp, ϕg, maxiter, minstep, minfunc, verbose, localsearch, n)

    obj = x -> func(x, args...; kwargs...)
    cons = make_constraints(constraints, args, kwargs, verbose)
    is_feasible = x -> all(cons(x) .>= 0)

    S = particles.swarmsize
    D = particles.particlelength
    verbose && println("Initialize current particle function values.")
    particles.fx = [obj(particles.x[i]) for i = 1:S]  # current particle function values
    particles.fs = [is_feasible(particles.x[i]) for i = 1:S]  # feasibility of each particle

    particles.history = []
    particles.fp = copy(particles.fx)  # best particle function values
    particles.order = sortperm(particles.fx, rev=true) #particle neighborhood order

    particles.g = copy(particles.x)  # best neighborhood swarm position
    particles.fg = ones(S) * Inf  # best neighborhood swarm position starting value

    # procurar alguma solução feasible e depois
    # pegar a com o menor valor dentre estas
    if .|(particles.fs...)
        aux1, aux2 = findmin(particles.fx[particles.fs])
        particles.best_position = copy(particles.x[particles.fs][aux2]) # best position
        particles.best_value = aux1 # best position value
    else
        aux1, aux2 = findmin(particles.fx)
        particles.best_position = copy(particles.x[aux2]) # best position
        particles.best_value = aux1 # best position value
    end

    # Store particle's best position (if constraints are satisfied)
    verbose && println("Store particle's best position.")
    update_position!(particles)

    # Update swarm's best position
    verbose && println("Update swarm's best position.")
    if localsearch
        findBestLocal!(particles, n)
    else
        findBestGlobal!(particles)
    end

    it = 1
    it_best = 1
    verbose && println("Start iteration.")
    while it <= maxiter

        # Update the particles' velocities and positions using equations 1 and 2
        # Eq(1): v = ω { v + ϕp . rp . (p -x) + ϕg . rg . (g - x)}
        # Eq(2): x = x +v
        rp = rand(particles.swarmsize) # random value
        rg = rand(particles.swarmsize) # random value
        particles.v = ψ .* ( particles.v .+ (ϕp .* rp .* (particles.p .- particles.x)) .+
        (ϕg .* rg .* (particles.g .- particles.x)) )
        particles.x = particles.x .+ particles.v

        # Correct for bound violations
        maskl = [particles.x[i] .< particles.lb for i = 1:S]
        masku = [particles.x[i] .> particles.ub for i = 1:S]
        mask = [ maskl[i] .| masku[i] for i = 1:S]
        particles.x = map((x,y) -> x.*(.~(y)), particles.x, mask)
        particles.x .+= [particles.lb.*maskl[i] for i = 1:S]
        particles.x .+= [particles.ub.*masku[i] for i=1:S]

        # Update objectives and constraints
        particles.fx = map(x-> obj(x), particles.x)
        particles.fs = map(x-> is_feasible(x), particles.x)

        # Store particle's best position (if constraints are satisfied)
        update_position!(particles)

        # Update swarm's best position
        if localsearch
            findBestLocal!(particles, n)
        else
            findBestGlobal!(particles)
        end

        # Compare swarm's best position with global best position
        ## procura mínimo
        i_min = findmin(particles.fp)
        position_temp = copy(particles.p[i_min[2]])
        ## se existe alguma solução feasible troca o mínimo pelo da feasible
        p_feasible = [is_feasible(particles.p[i]) for i = 1:S]
        if .|(p_feasible...)
            i_min = findmin(particles.fp[p_feasible])
            position_temp = particles.p[p_feasible][i_min[2]]
        end

        if i_min[1] < particles.best_value
            it_best = it
            stepsize = √(sum((particles.best_position .- position_temp).^2))
            if abs.(particles.best_value .- i_min[1]) <= minfunc
                verbose && println("\n Stopping search: Swarm best objective change
                less than $(minfunc)")
                break
            end
            if stepsize <= minstep
                verbose && println("Stopping search: Swarm best position change
                less than $(minstep)")
                break
            end
            # if feasible change value
            if is_feasible(position_temp)
                particles.best_position = copy(position_temp) # best position
                particles.best_value = i_min[1] # best position value
            elseif !is_feasible(position_temp) && !is_feasible(particles.best_position)
                particles.best_position = copy(position_temp) # best position
                particles.best_value = i_min[1] # best position value
            end
        end
        push!(particles.history, particles.best_value)
        verbose && (print("\r Iteration : $(it)/$(maxiter) | best-iteration $(it_best)
        | Best value: $(particles.best_value) "))
        it += 1
    end
    verbose && println("\n Stopping search: maximum iterations reached --> $(maxiter)")
    is_feasible(particles.best_position) || print("However, the optimization couldn't
     find a feasible design. Sorry")
    return nothing
end

# fun => função a ser otimizada
# lb => limite inferior
#ub => limite superior
# swarmsize => total de partículas

function pso(func::Function, lb::Vector, ub::Vector; constraints=nothing, args=(), kwargs=Dict(),
     swarmsize=100, psi=0.5, phip=0.5, phig=0.5, maxiter=100, minstep=1e-8, minfunc=1e-8,
            verbose=false,  localsearch = false, neighborhood = 2)

    @assert length(ub) == length(lb) "ub and lb must have same dimension"
    @assert all(ub .> lb) "Each value of 'ub' must be greater than 'lb'"
    particles = Particles(swarmsize, lb, ub) #create particles

    @assert iseven(neighborhood) "The value of 'neighborhood' must be even"
    @assert (particles.swarmsize >= neighborhood) " 'swarmsize' must be greater
     than 'neighborhood'"

    neighborhood = Integer(neighborhood/2)
    pso(func, particles, constraints, args, kwargs,
    psi, phip, phig, maxiter, minstep, minfunc, verbose,localsearch, neighborhood)
    return particles
end

function pso(particles:: Particles, func::Function; constraints=nothing, args=(), kwargs=Dict(),
    psi=0.5, phip=0.5, phig=0.5, maxiter=100, minstep=1e-8, minfunc=1e-8,
    verbose=false, localsearch = false, neighborhood = 2)

    @assert iseven(neighborhood) "The value of 'neighborhood' must be even"
    @assert (particles.swarmsize >= neighborhood) " 'swarmsize' must be greater
     than 'neighborhood'"

    neighborhood= Integer(neighborhood/2)

    pso(func, particles, constraints, args, kwargs,
        psi, phip, phig, maxiter, minstep, minfunc, verbose, localsearch, neighborhood)

end


text:
using Distributions
using Test, ForwardDiff

# Currently, most of the tests for NegativeBinomail are in the "ref" folder.
# Eventually, we might want to consolidate the tests here

mydiffp(r, p, k) = r/p - k/(1 - p)

@testset "NegativeBinomial r=$r, p=$p, k=$k" for
    p in exp10.(-10:0) .- eps(), # avoid p==1 since it's not differentiable
        r in exp10.(range(-10, stop=2, length=25)),
            k in (0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024)

    @test ForwardDiff.derivative(_p -> logpdf(NegativeBinomial(r, _p), k), p) ≈ mydiffp(r, p, k) rtol=1e-12 atol=1e-12
end

@testset "Check the corner case p==1" begin
    @test logpdf(NegativeBinomial(0.5, 1.0), 0) === 0.0
    @test logpdf(NegativeBinomial(0.5, 1.0), 1) === -Inf
end


text:
# Skipping the argument processing
# Assuming that the input file is quotes.tsv
# and using the test database

Using DBI
Using MySQL

escticks(s) = replace(s,"'","''");
qq = readdlm("quotes.tsv",'\t');
n = size(qq)[1];
conn = connect(MySQL5,"localhost","malcolm","mypasswd","test");

# We can create the categories table and quotes tables:
sql  = "create table categories ";
sql *= "(id integer primary key, ";
sql *= catname varchar(40) not null);
execute(prepare(conn,sql);
if (errcode(conn) > 0)
  error("Can't build categories table")
end;

# Note: Similar code for building the quotes table ...

j = 0;
cats = Dict(String,Int64);

# Main loop to load up the quotes table
for i = 1:n
  cat = qq[i,1];
  if haskey(cats,cat)
    jd = cats[cat];
  else
    j = j + 1; jd = j;
    cats[cat] = jd;
  end
  sql = "insert into quotes values($i,$jd,";
  if (length(qq[i,2]) > 0)
    sql *= string("'", escticks(qq[i,2]), "',");
  else
    sql *= string("null,");
  end
  sql *= string("'", escticks(qq[i,3]), "');");
  stmt = prepare(conn,sql); execute(stmt);
  if errcode(conn) > 0 error(errstring(conn); end;
end

# Now dump the categories
for cat = keys(cats)
  jd = cats[cat];
  sql = "insert into categories($jd,'$cat')";
  stmt = prepare(conn,sql); execute(stmt);
  if errcode(conn) > 0 error(errstring(conn); end;
end;
disconnect(conn);

max_stars_repo_path: Module 3/Chapter08/mysql-dbi.jl

max_issues_repo_path: Module 3/Chapter08/mysql-dbi.jl

max_forks_repo_path: Module 3/Chapter08/mysql-dbi.jl


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_translated)" --startup-file=no %s 2>&1 | FileCheck %s

function calls()
    f = true ? (+) : (-)
    return f(1, 1)
end
emit(calls)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.calls)}, svec(), calls() in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/translate/calls.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.calls)}"(%arg0: !jlir<"typeof(Main.calls)">) -> !jlir.Int64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     "jlir.goto"()[^bb2] : () -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"2">} : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.return"(%0) : (!jlir.Int64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
#=
# A struct containing face info
=#
export FaceData

mutable struct FaceData
    id;                  # face id
    normal;   # normal of this side of this face
    nodes;      	 # nodes on this side of this face
    isbdy;
    FaceData(id, face2v, normals) = (
		(normal, nodes, isbdy) = init(id, face2v, normals);
		new(id,normal,nodes,isbdy);
	)
    FaceData(id,normal,nodes,isbdy) = (
		new(id,normal,nodes,isbdy);
	)

end

function init(id, face2v, nomals)

	nodes = face2v[:,:,id];
	normal = [normals[:,id] -normals[:,id]];
	isbdy = 1;
	return [normal, nodes, isbdy];

end


text:
"""
    This struct keeps track of all uses of some mutable struct allocated
    in the current function. `uses` are all instances of `getfield` on the
    struct. `defs` are all instances of `setfield!` on the struct. The terminology
    refers to the uses/defs of the ``slot bundle'' that the mutable struct represents.

    In addition we keep track of all instances of a foreigncall preserve of this mutable
    struct. Somewhat counterintuitively, we don't actually need to make sure that the
    struct itself is live (or even allocated) at a ccall site. If there are no other places
    where the struct escapes (and thus e.g. where its address is taken), it need not be
    allocated. We do however, need to make sure to preserve any elments of this struct.
"""
struct SSADefUse
    uses::Vector{Int}
    defs::Vector{Int}
    ccall_preserve_uses::Vector{Int}
end
SSADefUse() = SSADefUse(Int[], Int[], Int[])

function try_compute_fieldidx_expr(@nospecialize(typ), @nospecialize(use_expr))
    field = use_expr.args[3]
    isa(field, QuoteNode) && (field = field.value)
    isa(field, Union{Int, Symbol}) || return nothing
    return try_compute_fieldidx(typ, field)
end

function lift_defuse(cfg::CFG, ssa::SSADefUse)
    # We remove from `uses` any block where all uses are dominated
    # by a def. This prevents insertion of dead phi nodes at the top
    # of such a block if that block happens to be in a loop
    ordered = Tuple{Int, Int, Bool}[(x, block_for_inst(cfg, x), true) for x in ssa.uses]
    for x in ssa.defs
        push!(ordered, (x, block_for_inst(cfg, x), false))
    end
    ordered = sort(ordered, by=x->x[1])
    bb_defs = Int[]
    bb_uses = Int[]
    last_bb = last_def_bb = 0
    for (_, bb, is_use) in ordered
        if bb != last_bb && is_use
            push!(bb_uses, bb)
        end
        last_bb = bb
        if last_def_bb != bb && !is_use
            push!(bb_defs, bb)
            last_def_bb = bb
        end
    end
    SSADefUse(bb_uses, bb_defs, Int[])
end

function find_curblock(domtree::DomTree, allblocks, curblock::Int)
    # TODO: This can be much faster by looking at current level and only
    # searching for those blocks in a sorted order
    while !(curblock in allblocks)
        curblock = domtree.idoms[curblock]
    end
    return curblock
end

function val_for_def_expr(ir::IRCode, def::Int, fidx::Int)
    if isexpr(ir[SSAValue(def)], :new)
        return ir[SSAValue(def)].args[1+fidx]
    else
        # The use is whatever the setfield was
        return ir[SSAValue(def)].args[4]
    end
end

function compute_value_for_block(ir::IRCode, domtree::DomTree, allblocks, du, phinodes, fidx, curblock)
    curblock = find_curblock(domtree, allblocks, curblock)
    def = reduce(max, 0, stmt for stmt in du.defs if block_for_inst(ir.cfg, stmt) == curblock)
    def == 0 ? phinodes[curblock] : val_for_def_expr(ir, def, fidx)
end

function compute_value_for_use(ir::IRCode, domtree::DomTree, allblocks, du, phinodes, fidx, use_idx)
    # Find the first dominating def
    curblock = stmtblock = block_for_inst(ir.cfg, use_idx)
    curblock = find_curblock(domtree, allblocks, curblock)
    defblockdefs = [stmt for stmt in du.defs if block_for_inst(ir.cfg, stmt) == curblock]
    def = 0
    if !isempty(defblockdefs)
        if curblock != stmtblock
            # Find the last def in this block
            def = maximum(defblockdefs)
        else
            # Find the last def before our use
            def = mapreduce(x->x >= use_idx ? 0 : x, max, defblockdefs)
        end
    end
    if def == 0
        if !haskey(phinodes, curblock)
            # If this happens, we need to search the predecessors for defs. Which
            # one doesn't matter - if it did, we'd have had a phinode
            return compute_value_for_block(ir, domtree, allblocks, du, phinodes, fidx, first(ir.cfg.blocks[stmtblock].preds))
        end
        # The use is the phinode
        return phinodes[curblock]
    else
        return val_for_def_expr(ir, def, fidx)
    end
end

function walk_to_def(compact::IncrementalCompact, @nospecialize(def), intermediaries=IdSet{Int}(), allow_phinode::Bool=true, phi_locs=Tuple{Int, Int}[])
    if !isa(def, SSAValue)
        return (def, 0)
    end
    orig_defidx = defidx = def.id
    # Step 2: Figure out what the struct is defined as
    def = compact[defidx]
    typeconstraint = types(compact)[defidx]
    ## Track definitions through PiNode/PhiNode
    found_def = false
    ## Track which PhiNodes, SSAValue intermediaries
    ## we forwarded through.
    while true
        if isa(def, PiNode)
            push!(intermediaries, defidx)
            typeconstraint = typeintersect(typeconstraint, def.typ)
            if isa(def.val, SSAValue)
                defidx = def.val.id
                def = compact[defidx]
            else
                def = def.val
            end
            continue
        elseif isa(def, FastForward)
            append!(phi_locs, def.phi_locs)
            def = def.to
        elseif isa(def, PhiNode)
            # For now, we don't track setfields structs through phi nodes
            allow_phinode || break
            push!(intermediaries, defidx)
            possible_predecessors = let def=def, typeconstraint=typeconstraint
                collect(Iterators.filter(1:length(def.edges)) do n
                    isassigned(def.values, n) || return false
                    value = def.values[n]
                    edge_typ = widenconst(compact_exprtype(compact, value))
                    return typeintersect(edge_typ, typeconstraint) !== Union{}
                end)
            end
            # For now, only look at unique predecessors
            if length(possible_predecessors) == 1
                n = possible_predecessors[1]
                pred = def.edges[n]
                val = def.values[n]
                if isa(val, SSAValue)
                    push!(phi_locs, (pred, defidx))
                    defidx = val.id
                    def = compact[defidx]
                elseif def == val
                    # This shouldn't really ever happen, but
                    # patterns like this can occur in dead code,
                    # so bail out.
                    break
                else
                    def = val
                end
                continue
            end
        elseif isa(def, SSAValue)
            push!(intermediaries, defidx)
            defidx = def.id
            def = compact[def.id]
            continue
        end
        found_def = true
        break
    end
    found_def ? (def, defidx) : nothing
end

function process_immutable_preserve(new_preserves::Vector{Any}, compact::IncrementalCompact, def::Expr)
    for arg in (isexpr(def, :new) ? def.args : def.args[2:end])
        if !isbitstype(widenconst(compact_exprtype(compact, arg)))
            push!(new_preserves, arg)
        end
    end
end

struct FastForward
    to::SSAValue
    phi_locs::Vector{Tuple{Int64, Int64}}
end

function getfield_elim_pass!(ir::IRCode, domtree::DomTree)
    compact = IncrementalCompact(ir)
    insertions = Vector{Any}()
    defuses = IdDict{Int, Tuple{IdSet{Int}, SSADefUse}}()
    for (idx, stmt) in compact
        isa(stmt, Expr) || continue
        is_getfield = false
        is_ccall = false
        # Step 1: Check whether the statement we're looking at is a getfield/setfield!
        if is_known_call(stmt, setfield!, compact)
            is_setfield = true
        elseif is_known_call(stmt, getfield, compact)
            is_getfield = true
        elseif isexpr(stmt, :foreigncall)
            nccallargs = stmt.args[5]
            new_preserves = Any[]
            old_preserves = stmt.args[(6+nccallargs):end]
            for (pidx, preserved_arg) in enumerate(old_preserves)
                intermediaries = IdSet()
                isa(preserved_arg, SSAValue) || continue
                def = walk_to_def(compact, preserved_arg, intermediaries, false)
                def !== nothing || continue
                (def, defidx) = def
                if is_tuple_call(compact, def)
                    process_immutable_preserve(new_preserves, compact, def)
                    old_preserves[pidx] = nothing
                    continue
                elseif isexpr(def, :new)
                    typ = def.typ
                    if isa(typ, UnionAll)
                        typ = unwrap_unionall(typ)
                    end
                    if !typ.mutable
                        process_immutable_preserve(new_preserves, compact, def)
                        old_preserves[pidx] = nothing
                        continue
                    end
                else
                    continue
                end
                mid, defuse = get!(defuses, defidx, (IdSet{Int}(), SSADefUse()))
                push!(defuse.ccall_preserve_uses, idx)
                union!(mid, intermediaries)
                continue
            end
            if !isempty(new_preserves)
                old_preserves = filter(ssa->ssa !== nothing, old_preserves)
                new_expr = Expr(:foreigncall, stmt.args[1:(6+nccallargs-1)]...,
                    old_preserves..., new_preserves...)
                new_expr.typ = stmt.typ
                compact[idx] = new_expr
            end
            continue
        else
            continue
        end
        ## Normalize the field argument to getfield/setfield
        field = stmt.args[3]
        isa(field, QuoteNode) && (field = field.value)
        isa(field, Union{Int, Symbol}) || continue

        intermediaries = IdSet()
        phi_locs = Tuple{Int, Int}[]
        def = walk_to_def(compact, stmt.args[2], intermediaries, is_getfield, phi_locs)
        def !== nothing || continue
        (def, defidx) = def

        if !is_getfield
            (defidx == 0) && continue
            mid, defuse = get!(defuses, defidx, (IdSet{Int}(), SSADefUse()))
            push!(defuse.defs, idx)
            union!(mid, intermediaries)
            continue
        end
        # Step 3: Check if the definition we eventually end up at is either
        # a tuple(...) call or Expr(:new) and perform replacement.
        if is_tuple_call(compact, def) && isa(field, Int) && 1 <= field < length(def.args)
            forwarded = def.args[1+field]
        elseif isexpr(def, :new)
            typ = def.typ
            if isa(typ, UnionAll)
                typ = unwrap_unionall(typ)
            end
            isa(typ, DataType) || continue
            if typ.mutable
                @assert defidx != 0
                mid, defuse = get!(defuses, defidx, (IdSet{Int}(), SSADefUse()))
                push!(defuse.uses, idx)
                union!(mid, intermediaries)
                continue
            end
            field = try_compute_fieldidx_expr(typ, stmt)
            field === nothing && continue
            forwarded = def.args[1+field]
        else
            obj = compact_exprtype(compact, def)
            isa(obj, Const) || continue
            obj = obj.val
            isimmutable(obj) || continue
            field = try_compute_fieldidx_expr(typeof(obj), stmt)
            field === nothing && continue
            isdefined(obj, field) || continue
            val = getfield(obj, field)
            is_inlineable_constant(val) || continue
            forwarded = quoted(val)
        end
        # Step 4: Remember any phinodes we need to insert
        if !isempty(phi_locs) && isa(forwarded, SSAValue)
            # TODO: We have have to use BB ids for phi_locs
            # to avoid index invalidation.
            push!(insertions, idx)
            compact[idx] = FastForward(forwarded, phi_locs)
        else
            compact[idx] = forwarded
        end
    end
    ir = finish(compact)
    # Now go through any mutable structs and see which ones we can eliminate
    for (idx, (intermediaries, defuse)) in defuses
        intermediaries = collect(intermediaries)
        # Check if there are any uses we did not account for. If so, the variable
        # escapes and we cannot eliminate the allocation. This works, because we're guaranteed
        # not to include any intermediaries that have dead uses. As a result, missing uses will only ever
        # show up in the nuses_total count.
        nleaves = length(defuse.uses) + length(defuse.defs) + length(defuse.ccall_preserve_uses)
        nuses_total = compact.used_ssas[idx] + mapreduce(idx->compact.used_ssas[idx], +, 0, intermediaries) - length(intermediaries)
        nleaves == nuses_total || continue
        # Find the type for this allocation
        defexpr = ir[SSAValue(idx)]
        isexpr(defexpr, :new) || continue
        typ = defexpr.typ
        if isa(typ, UnionAll)
            typ = unwrap_unionall(typ)
        end
        # Could still end up here if we tried to setfield! and immutable, which would
        # error at runtime, but is not illegal to have in the IR.
        typ.mutable || continue
        # Partition defuses by field
        fielddefuse = SSADefUse[SSADefUse() for _ = 1:fieldcount(typ)]
        ok = true
        for use in defuse.uses
            field = try_compute_fieldidx_expr(typ, ir[SSAValue(use)])
            field === nothing && (ok = false; break)
            push!(fielddefuse[field].uses, use)
        end
        ok || continue
        for use in defuse.defs
            field = try_compute_fieldidx_expr(typ, ir[SSAValue(use)])
            field === nothing && (ok = false; break)
            push!(fielddefuse[field].defs, use)
        end
        ok || continue
        preserve_uses = IdDict{Int, Vector{Any}}((idx=>Any[] for idx in IdSet{Int}(defuse.ccall_preserve_uses)))
        # Everything accounted for. Go field by field and perform idf
        for (fidx, du) in pairs(fielddefuse)
            ftyp = fieldtype(typ, fidx)
            if !isempty(du.uses)
                push!(du.defs, idx)
                ldu = compute_live_ins(ir.cfg, du)
                phiblocks = []
                if !isempty(ldu.live_in_bbs)
                    phiblocks = idf(ir.cfg, ldu, domtree)
                end
                phinodes = IdDict{Int, SSAValue}()
                for b in phiblocks
                    n = PhiNode()
                    phinodes[b] = insert_node!(ir, first(ir.cfg.blocks[b].stmts), ftyp, n)
                end
                # Now go through all uses and rewrite them
                allblocks = sort(vcat(phiblocks, ldu.def_bbs))
                for stmt in du.uses
                    ir[SSAValue(stmt)] = compute_value_for_use(ir, domtree, allblocks, du, phinodes, fidx, stmt)
                end
                if !isbitstype(fieldtype(typ, fidx))
                    for (use, list) in preserve_uses
                        push!(list, compute_value_for_use(ir, domtree, allblocks, du, phinodes, fidx, use))
                    end
                end
                for b in phiblocks
                    for p in ir.cfg.blocks[b].preds
                        n = ir[phinodes[b]]
                        push!(n.edges, p)
                        push!(n.values, compute_value_for_block(ir, domtree,
                            allblocks, du, phinodes, fidx, p))
                    end
                end
            end
            for stmt in du.defs
                stmt == idx && continue
                ir[SSAValue(stmt)] = nothing
            end
            continue
        end
        isempty(defuse.ccall_preserve_uses) && continue
        push!(intermediaries, idx)
        # Insert the new preserves
        for (use, new_preserves) in preserve_uses
            useexpr = ir[SSAValue(use)]
            nccallargs = useexpr.args[5]
            old_preserves = filter(ssa->!isa(ssa, SSAValue) || !(ssa.id in intermediaries), useexpr.args[(6+nccallargs):end])
            new_expr = Expr(:foreigncall, useexpr.args[1:(6+nccallargs-1)]...,
                old_preserves..., new_preserves...)
            new_expr.typ = useexpr.typ
            ir[SSAValue(use)] = new_expr
        end
    end
    for idx in insertions
        # For non-dominating load-store forward, we may have to insert extra phi nodes
        # TODO: Can use the domtree to eliminate unnecessary phis, but ok for now
        ff = ir.stmts[idx]
        ff === nothing && continue # May have been DCE'd if there were no more uses
        ff = ff::FastForward
        forwarded = ff.to
        if isa(forwarded, SSAValue)
            forwarded_typ = ir.types[forwarded.id]
            for (pred, pos) in reverse!(ff.phi_locs)
                node = PhiNode()
                push!(node.edges, pred)
                push!(node.values, forwarded)
                forwarded = insert_node!(ir, pos, forwarded_typ, node)
            end
        end
        ir.stmts[idx] = forwarded
    end
    ir
end

function type_lift_pass!(ir::IRCode)
    type_ctx_uses = Vector{Vector{Int}}[]
    has_non_type_ctx_uses = IdSet{Int}()
    lifted_undef = IdDict{Int, Any}()
    for (idx, stmt) in pairs(ir.stmts)
        if stmt isa Expr && (stmt.head === :isdefined || stmt.head === :undefcheck)
            val = (stmt.head === :isdefined) ? stmt.args[1] : stmt.args[2]
            # undef can only show up by being introduced in a phi
            # node (or an UpsilonNode() argument to a PhiC node),
            # so lift all these nodes that have maybe undef values
            processed = IdDict{Int, Union{SSAValue, Bool}}()
            if !isa(val, SSAValue)
                if stmt.head === :undefcheck
                    ir.stmts[idx] = nothing
                end
                continue
            end
            worklist = Tuple{Int, Int, SSAValue, Int}[(val.id, 0, SSAValue(0), 0)]
            stmt_id = val.id
            while isa(ir.stmts[stmt_id], PiNode)
                stmt_id = ir.stmts[stmt_id].val.id
            end
            def = ir.stmts[stmt_id]
            if !isa(def, PhiNode) && !isa(def, PhiCNode)
                if stmt.head === :isdefined
                    ir.stmts[idx] = true
                else
                    ir.stmts[idx] = nothing
                end
                continue
            end
            if !haskey(lifted_undef, stmt_id)
                first = true
                while !isempty(worklist)
                    item, w_up_id, which, use = pop!(worklist)
                    def = ir.stmts[item]
                    if isa(def, PhiNode)
                        edges = copy(def.edges)
                        values = Vector{Any}(undef, length(edges))
                        new_phi = length(values) == 0 ? false : insert_node!(ir, item, Bool, PhiNode(edges, values))
                    else
                        values = Vector{Any}(undef, length(def.values))
                        new_phi = length(values) == 0 ? false : insert_node!(ir, item, Bool, PhiCNode(values))
                    end
                    processed[item] = new_phi
                    if first
                        lifted_undef[stmt_id] = new_phi
                        first = false
                    end
                    local id::Int = 0
                    for i = 1:length(values)
                        if !isassigned(def.values, i)
                            val = false
                        elseif !isa(def.values[i], SSAValue)
                            val = true
                        else
                            up_id = id = def.values[i].id
                            @label restart
                            if !isa(ir.types[id], MaybeUndef)
                                val = true
                            else
                                if isa(ir.stmts[id], UpsilonNode)
                                    up = ir.stmts[id]
                                    if !isdefined(up, :val)
                                        val = false
                                    elseif !isa(up.val, SSAValue)
                                        val = true
                                    else
                                        id = up.val.id
                                        @goto restart
                                    end
                                else
                                    while isa(ir.stmts[id], PiNode)
                                        id = ir.stmts[id].val.id
                                    end
                                    if isa(ir.stmts[id], Union{PhiNode, PhiCNode})
                                        if haskey(processed, id)
                                            val = processed[id]
                                        else
                                            push!(worklist, (id, up_id, new_phi, i))
                                            continue
                                        end
                                    else
                                        val = true
                                    end
                                end
                            end
                        end
                        if isa(def, PhiNode)
                            values[i] = val
                        else
                            values[i] = insert_node!(ir, up_id, Bool, UpsilonNode(val))
                        end
                    end
                    if which !== SSAValue(0)
                        phi = ir[which]
                        if isa(phi, PhiNode)
                            phi.values[use] = new_phi
                        else
                            phi = phi::PhiCNode
                            ir[which].values[use] = insert_node!(ir, w_up_id, Bool, UpsilonNode(new_phi))
                        end
                    end
                end
            end
            if stmt.head === :isdefined
                ir.stmts[idx] = lifted_undef[stmt_id]
            else
                ir.stmts[idx] = Expr(:throw_undef_if_not, stmt.args[1], lifted_undef[stmt_id])
            end
        end
    end
    ir
end


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

using PyPlot

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

const a=2.9#3.11
const b=3.7#3.62
const d=8.7#4.84
const f0=30e6
const f1=2000e6
rf=.005
nf=int(1/log(1+rf)*log(f1/f0))
f=zeros(nf)
for i=1:nf
  f[i]=f0*(1+rf)^i
end

const V=a*b*d
const Q4=16*pi^2*V*f.^3/(1*c^3)
const Q1=3*V./(2*1*(2./sqrt(2*pi*f.*mur*mu0*1*1e5)))/200
const Q=1./(1./Q1+1./Q4)
modes=zeros(nf)
const x=1
const y=2
const z=1

E0=1
H0=1/(c*mu0*mur)
Ex=zeros(Complex128,nf)
Ey=zeros(Complex128,nf)
Ez=zeros(Complex128,nf)
for u=1:nf
  n_x=int(f[u]*(1+300/Q[u])*2*a/c)
  n_y=int(f[u]*(1+300/Q[u])*2*b/c)
  n_z=int(f[u]*(1+300/Q[u])*2*d/c)
  M=[]
  N=[]
  P=[]
  F=[]
  AMP=[]
  for mm=0:n_x
    for nn=0:n_y
      for pp=0:n_z
        fmnp=c/2*sqrt((mm/a)^2+(nn/b)^2+(pp/d)^2)
        Amp=urc(Q[u],f[u],fmnp)
        if Amp>0.005
          M=vcat(M,mm)
          N=vcat(N,nn)
          P=vcat(P,pp)
          F=vcat(F,fmnp)
          AMP=vcat(AMP,Amp)
        end
      end
    end
  end
  ETMx_x=0
  ETMx_y=0
  ETMx_z=0
  ETMy_x=0
  ETMy_y=0
  ETMy_z=0
  ETMz_x=0
  ETMz_y=0
  ETMz_z=0
  ETEx_y=0
  ETEx_z=0
  ETEy_x=0
  ETEy_z=0
  ETEz_y=0
  ETEz_x=0
  nm=length(M)
  println("$(round(f[u]*100/1e6)/100) MHz, $nm resonant freq.")
  for k=1:nm
    fmnp=F[k]
    Amp=AMP[k]
    kx=M[k]*pi/a
    ky=N[k]*pi/b
    kz=P[k]*pi/d
    k2mnp=kx^2+ky^2+kz^2
    ssc_z=sin(kx*x)*sin(ky*y)*cos(kz*z)
    scs_z=sin(kx*x)*cos(ky*y)*sin(kz*z)
    css_z=cos(kx*x)*sin(ky*y)*sin(kz*z)

    ssc_x=sin(ky*y)*sin(kz*z)*cos(kx*x)
    scs_x=sin(ky*y)*cos(kz*z)*sin(kx*x)
    css_x=cos(ky*y)*sin(kz*z)*sin(kx*x)

    ssc_y=sin(kz*z)*sin(kx*x)*cos(ky*y)
    scs_y=sin(kz*z)*cos(kx*x)*sin(ky*y)
    css_y=cos(kz*z)*sin(kx*x)*sin(ky*y)

    if M[k]>0 && N[k]>0 && P[k]>0
      #TM modes z
      ETMx_z += -E0*Amp*kx*kz/(k2mnp-kz^2)*css_z
      ETMy_z += E0*Amp*ky*kz/(k2mnp-kz^2)*scs_z
      ETMz_z += E0*Amp*ssc_z
      #TE modes z
      ETEx_z += -1im*2*pi*fmnp*mur*mu0*ky*H0*Amp/(k2mnp-kz^2)*css_z
      ETEy_z += 1im*2*pi*fmnp*mur*mu0*kx*H0*Amp/(k2mnp-kz^2)*scs_z
      modes[u] += 2

      #TM modes x
      ETMy_x += -E0*Amp*ky*kx/(k2mnp-kx^2)*css_x
      ETMz_x += E0*Amp*kz*kx/(k2mnp-kx^2)*scs_x
      ETMx_x += E0*Amp*ssc_x
      #TE modes x
      ETEy_x += -1im*2*pi*fmnp*mur*mu0*kz*H0*Amp/(k2mnp-kx^2)*css_x
      ETEz_x += 1im*2*pi*fmnp*mur*mu0*ky*H0*Amp/(k2mnp-kx^2)*scs_x
      modes[u] += 2

      ##TM modes y
      ETMz_y += -E0*Amp*kz*ky/(k2mnp-ky^2)*css_y
      ETMx_y += E0*Amp*kx*ky/(k2mnp-ky^2)*scs_y
      ETMy_y += E0*Amp*ssc_y
      ##TE modes y
      ETEz_y += -1im*2*pi*fmnp*mur*mu0*kx*H0*Amp/(k2mnp-ky^2)*css_y
      ETEx_y += 1im*2*pi*fmnp*mur*mu0*kz*H0*Amp/(k2mnp-ky^2)*scs_y
      modes[u] += 2
    end

    if M[k]>0 && N[k]>0 && P[k]==0
      ETMz_z+=E0*Amp*ssc_z #TM mode z
      ETEz_y+=-1im*2*pi*fmnp*mur*mu0*kx*H0*Amp/(k2mnp-ky^2)*css_y #TE mode y
      ETEz_x+=1im*2*pi*fmnp*mur*mu0*ky*H0*Amp/(k2mnp-kx^2)*scs_x #TE mode x
      modes[u]+=3
    end
    if M[k]==0 && N[k]>0 && P[k]>0
      ETEx_z+=-1im*2*pi*fmnp*mur*mu0*ky*H0*Amp/(k2mnp-kz^2)*css_z #TE mode z
      ETMx_x+=E0*Amp*ssc_x #TM mode x
      ETEx_y+=1im*2*pi*fmnp*mur*mu0*kz*H0*Amp/(k2mnp-ky^2)*scs_y #TE mode y
      modes[u]+=3
    end
    if M[k]>0 && N[k]==0 && P[k]>0
      ETEy_z+=1im*2*pi*fmnp*mur*mu0*kx*H0*Amp/(k2mnp-kz^2)*scs_z #TE mode z
      ETMy_y+=E0*Amp*ssc_y #TM mode y
      ETEy_x+=-1im*2*pi*fmnp*mur*mu0*kz*H0*Amp/(k2mnp-kx^2)*css_x #TE mode x
      modes[u]+=3
    end
  end
  Ex[u]=ETMx_x+ETMx_y+ETMx_z+ETEx_y+ETEx_z
  Ey[u]=ETMy_x+ETMy_y+ETMy_z+ETEy_x+ETEy_z
  Ez[u]=ETMz_x+ETMz_y+ETMz_z+ETEz_y+ETEz_x
end

figure(num=2,figsize=(8,10))
subplot(311)
title("\$ \\vert E_x\\vert\$")
loglog(f,abs(Ex))#,label="\$ \\vert E_x\\vert\$")
ylabel("A.U")
xlim(f0,f1)
grid()
subplot(312)
title("\$ \\vert E_y\\vert\$")
loglog(f,abs(Ey))#,label="\$ \\vert E_x\\vert\$")
ylabel("A.U")
xlim(f0,f1)
grid()
subplot(313)
title("\$ \\vert E_z\\vert\$")
loglog(f,abs(Ez))#,label="\$ \\vert E_x\\vert\$")
ylabel("A.U")
xlim(f0,f1)
grid()
xlabel("\$f\$ in Hz")
#legend()
savefig("fig2.pdf",bbox_inches="tight")


text:
total_mass(n::Integer, mass::Number) = n * mass
total_mass(mass::AbstractVector{<:Number}) = sum(mass)
total_mass(pos::AbstractMatrix{<:Number}, mass::Number) = total_mass(size(pos, 2), mass)
total_mass(pos::AbstractMatrix{<:Number}, mass::AbstractVector{<:Number}) = total_mass(mass)

@doc raw"""
    angular_momentum(pos, vel, mass[, ellipsoid])
    angular_momentum(p::Particles[, ellipsoid])

Returns the total angular momentum ``J = \sum_i m_i \mathbf{r}_i \times \mathbf{v}_i`` within the
ellipsoid if given. Values are assumed to be in kpc, km/s, and solar masses.

```julia
using Unitful, UnitfulAstro

g::Galaxy # assume an already read-in galaxy with :full units
sph = Sphere(8u"kpc")

J = angular_momentum(g.stars, sph)
```
"""
function angular_momentum(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Union{Number,AbstractVector{<:Number}},
)
    return dropdims(sum(mass' .* pos × vel; dims=2); dims=2)
end

function angular_momentum(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::AbstractVector{<:Number},
    ellipsoid::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, ellipsoid)
    return angular_momentum(pos[:, mask], vel[:, mask], mass[mask])
end

function angular_momentum(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Number,
    ellipsoid::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, ellipsoid)
    return angular_momentum(pos[:, mask], vel[:, mask], mass)
end

function angular_momentum(p::Particles, ellipsoid::Union{Ellipsoid,Nothing}=nothing;)
    if isnothing(ellipsoid)
        return angular_momentum(p.pos, p.vel, p.mass)
    else
        return angular_momentum(p.pos, p.vel, p.mass, ellipsoid)
    end
end

@doc raw"""
    specific_angular_momentum(pos, vel, mass; Mtot)
    specific_angular_momentum(pos, vel, mass, ellipsoid)
    specific_angular_momentum(p::Particles[, ellipsoid]; Mtot)

Returns the specific angular momentum ``j = \sum_i m_i \mathbf{r}_i \times \mathbf{v}_i / \sum_i m_i``
within the ellipsoid if given. Values are assumed to be in kpc, km/s, and solar masses.
The total mass can be provided via `Mtot` if already known.

```julia
using Unitful, UnitfulAstro

g::Galaxy # assume an already read-in galaxy with :full units
sph = Sphere(8u"kpc")

j = specific_angular_momentum(g.stars, sph)
```
"""
function specific_angular_momentum(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Union{Number,AbstractVector{<:Number}};
    Mtot::Number=total_mass(pos, mass),
)
    return angular_momentum(pos, vel, mass) ./ Mtot
end

function specific_angular_momentum(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Number,
    ellipsoid::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, ellipsoid)
    return specific_angular_momentum(pos[:, mask], vel[:, mask], mass)
end

function specific_angular_momentum(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::AbstractVector{<:Number},
    ellipsoid::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, ellipsoid)
    return specific_angular_momentum(pos[:, mask], vel[:, mask], mass[mask])
end

function specific_angular_momentum(
    p::Particles,
    ellipsoid::Union{Ellipsoid,Nothing}=nothing;
    Mtot::Union{Number,Nothing}=nothing,
)
    if isnothing(ellipsoid)
        if isnothing(Mtot)
            return specific_angular_momentum(p.pos, p.vel, p.mass)
        else
            return specific_angular_momentum(p.pos, p.vel, p.mass; Mtot)
        end
    else
        return specific_angular_momentum(p.pos, p.vel, p.mass, ellipsoid)
    end
end

@doc raw"""
    b_value(pos, vel, mass[, Mtot])
    b_value(pos, vel, mass, ellipsoid)
    b_value(p::Particles[, ellipsoid])

Returns the b-value ``b = \log(j / (\mathrm{km/s})) - 2/3 \log(M_\mathrm{tot}/\mathrm{M}_\odot)``,
where ``j`` is the [`specific_angular_momentum`](@ref). Values are assumed to be in kpc, km/s, and
solar masses.

```julia
using Unitful, UnitfulAstro

g::Galaxy # assume an already read-in galaxy with :full units
sph = Sphere(8u"kpc")

b = b_value(g.stars, sph)
```
"""
function b_value(
    pos::AbstractMatrix{<:Real},
    vel::AbstractMatrix{<:Real},
    mass::Union{Number,AbstractVector{<:Real}},
    Mtot::Real,
)
    return log10(norm(specific_angular_momentum(pos, vel, mass; Mtot))) - 2 // 3 * log10(Mtot)
end

function b_value(
    pos::AbstractMatrix{<:Quantity},
    vel::AbstractMatrix{<:Quantity},
    mass::Union{Number,AbstractVector{<:Quantity}},
    Mtot::Quantity,
)
    return log10(ustrip(u"kpc * km/s", norm(specific_angular_momentum(pos, vel, mass; Mtot)))) -
           2 // 3 * log10(ustrip(u"Msun", Mtot))
end

function b_value(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Union{Number,AbstractVector{<:Number}},
)
    return b_value(pos, vel, mass, total_mass(pos, mass))
end

function b_value(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Number,
    ellipsoid::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, ellipsoid)
    return b_value(pos[:, mask], vel[:, mask], mass)
end

function b_value(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::AbstractVector{<:Number},
    ellipsoid::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, ellipsoid)
    return b_value(pos[:, mask], vel[:, mask], mass[mask])
end

b_value(p::Particles, ellipsoid::Nothing=nothing) = b_value(p.pos, p.vel, p.mass)
b_value(p::Particles, ellipsoid::Ellipsoid) = b_value(p.pos, p.vel, p.mass, ellipsoid)

@doc raw"""
    circular_velocity(Mtot, r)

Returns the circular velocity ``v_\mathrm{circ} = \sqrt{G M(<r) / r}``. Values are assumed to be in
solar masses and kpc.
"""
function circular_velocity(Mtot::T1, r::T2) where {T1<:AbstractFloat,T2<:AbstractFloat}
    return convert(promote_type(T1, T2), 1e-3sqrt(1.3271244e20Mtot / (3.0856775814913674e19r))) # in km/s
end
circular_velocity(Mtot::Integer, r::T) where {T<:AbstractFloat} = circular_velocity(convert(T, Mtot), r)
circular_velocity(Mtot::T, r::Integer) where {T<:AbstractFloat} = circular_velocity(Mtot, convert(T, r))
circular_velocity(Mtot::Integer, r::Integer) = circular_velocity(float(Mtot), float(r))

function circular_velocity(Mtot::Quantity, r::Quantity)
    vcirc = sqrt(ustrip(u"Msun", Mtot) * u"GMsun" / r)

    # convert to float type of vcirc (conversion to km/s defaults to Float64)
    return convert(typeof(vcirc).parameters[1], ustrip(u"km/s", vcirc)) * u"km/s"
end

@doc raw"""
    spin_parameter(j, Mtot)
    spin_parameter(pos, vel, mass, Mtot, r)
    spin_parameter(pos, vel, mass, Mtot, sphere)
    spin_parameter(p::Particles, Mtot, r; mask=nothing)

Returns the spin parameter ``\lambda(r) = j / (\sqrt{2} r v_\mathrm{circ}(r))``, where ``j`` is
the [`specific_angular_momentum`](@ref). The spin parameter is computed for only one particle component.

A mask can be applied to filter the particles.
"""
function spin_parameter(j::AbstractVector{<:Number}, Mtot::Number, r::Number)
    return norm(j) / (√2 * r * circular_velocity(Mtot, r))
end

function spin_parameter(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Union{Number,AbstractVector{<:Number}},
    Mtot::Number,
    r::Number,
)
    return spin_parameter(specific_angular_momentum(pos, vel, mass), Mtot, r)
end

function spin_parameter(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::Number,
    Mtot::Number,
    sphere::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, sphere)
    return spin_parameter(pos[:, mask], vel[:, mask], mass, Mtot, sphere.radius)
end

function spin_parameter(
    pos::AbstractMatrix{<:Number},
    vel::AbstractMatrix{<:Number},
    mass::AbstractVector{<:Number},
    Mtot::Number,
    sphere::Ellipsoid,
)
    mask = ellipsoidal_mask(pos, sphere)
    return spin_parameter(pos[:, mask], vel[:, mask], mass[mask], Mtot, sphere.radius)
end

function spin_parameter(p::Particles, Mtot::Number, r::Number; mask::Union{AbstractVector,Nothing}=nothing)
    if isnothing(mask)
        return spin_parameter(p.pos, p.vel, p.mass, Mtot, Sphere(r))
    end
    return @views spin_parameter(
        p.pos[:, mask],
        p.vel[:, mask],
        p.mass isa Number ? p.mass : p.mass[mask],
        Mtot,
        Sphere(r),
    )
end

@doc raw"""
    spin_parameter(g::AbstractGalaxy, Mtot, r)

Returns the total spin parameter ``\lambda(r) = j / (\sqrt{2} r v_\mathrm{circ}(r))``, where ``j`` is
the [`specific_angular_momentum`](@ref). The spin parameter is computed for all particle components.
"""
function spin_parameter(g::AbstractGalaxy, Mtot::Number, r::Number)
    sph = Sphere(r)
    J = angular_momentum(g.stars, sph) .+ angular_momentum(g.gas, sph) .+ angular_momentum(g.dm, sph)
    return spin_parameter(J ./ Mtot, Mtot, r)
end


text:
swapext(f, new_ext) = "$(splitext(f)[1])$new_ext"

isarm() = Base.Sys.ARCH in (:arm, :arm32, :arm7l, :armv7l, :arm8l, :armv8l, :aarch64, :arm64)

@noinline function isblank(img)
    return all(c -> green(c) == 0, img) ||
           all(c -> blue(c) == 0, img) ||
           all(c -> red(c) == 0, img) ||
           maximum(rawview(channelview(img))) < 0xcf
end

function compare_colors(a::RGB, b::RGB, tol)
    ok = true
    for f in (red, green, blue)
        dev = abs(float(f(a)) - float(f(b)))
        ok &= dev <= tol
    end
    return ok
end

# Helper functions
function test_compare_frames(test_frame, ref_frame, tol = 0.05)
    if isarm()
        @test_skip test_frame == ref_frame
    else
        frames_similar = true
        for (a, b) in zip(test_frame, ref_frame)
            frames_similar &= compare_colors(a, b, tol)
        end
        @test frames_similar
    end
end

# uses read!
function read_frameno!(img, v, frameno)
    seekstart(v)
    i = 0
    while !eof(v) && i < frameno
        read!(v, img)
        i += 1
    end
end

function make_comparison_frame_png(vidpath::AbstractString, frameno::Integer, writedir = tempdir())
    vid_basename = first(splitext(basename(vidpath)))
    png_name = joinpath(writedir, "$(vid_basename)_$(frameno).png")
    FFMPEG.exe(
        `-y -v error -i $(vidpath) -vf "sws_flags=accurate_rnd+full_chroma_inp+full_chroma_int; select=eq(n\,$(frameno-1))" -vframes 1 $(png_name)`,
    )
    return png_name
end

function make_comparison_frame_png(f, args...)
    png_name = make_comparison_frame_png(args...)
    try
        f(png_name)
    finally
        rm(png_name, force = true)
    end
end

function get_video_extrema(v)
    img = read(v)
    raw_img = parent(img)
    # Test that the limited range of this video is converted to full range
    minp, maxp = extrema(raw_img)
    while !eof(v)
        read!(v, raw_img)
        this_minp, this_maxp = extrema(raw_img)
        minp = min(minp, this_minp)
        maxp = max(maxp, this_maxp)
    end
    return minp, maxp
end

function get_raw_luma_extrema(elt, vidpath, nw, nh)
    buff, align = VideoIO.openvideo(vidpath) do v
        return VideoIO.read_raw(v, 1)
    end
    luma_buff = view(buff, 1:nw*nh*sizeof(elt))
    luma_vals = reinterpret(elt, luma_buff)
    return reinterpret.(extrema(luma_vals))
end

using ColorTypes: RGB, HSV
using FixedPointNumbers: Normed, N6f10
using Base: ReinterpretArray

function test_tone!(a::AbstractMatrix{X}, offset = 0, minval = 0, maxval = reinterpret(one(X))) where {T,X<:Normed{T}}
    maxcodept = reinterpret(one(X))
    modsize = maxval - minval + 1
    @inbounds for i in eachindex(a)
        a[i] = reinterpret(X, T(minval + mod(i + offset - 1, modsize)))
    end
    return a
end

function test_tone!(
    a::AbstractMatrix{C},
    offset = 0,
    minval = 0,
    maxval = reinterpet(one(X)),
) where {T,X<:Normed{T},C<:RGB{X}}
    modsize = maxval - minval + 1
    @inbounds for i in eachindex(a)
        h = mod(i, 360)
        v = minval + mod(i + offset - 1, modsize) / maxcodept
        hsv = HSV(h, 1, v)
        a[i] = convert(RGB{X}, hsv)
    end
    return a
end

test_tone(::Type{T}, nx::Integer, ny, args...) where {T} = test_tone!(Matrix{T}(undef, nx, ny), args...)

test_tone(nx::Integer, ny, args...) = test_tone(N6f10, nx, ny, args...)

function make_test_tones(::Type{T}, nx, ny, nf, args...) where {T}
    imgs = Vector{Matrix{T}}(undef, nf)
    @inbounds for i in 1:nf
        imgs[i] = test_tone(T, nx, ny, i - 1, args...)
    end
    return imgs
end

sizeof_parent(buf::Array) = sizeof(buf)
sizeof_parent(buf::ReinterpretArray) = sizeof_parent(parent(buf))

function copy_imgbuf_to_buf!(
    buf::StridedArray{UInt8},
    bwidth::Integer,
    fheight::Integer,
    imgbufp::Ptr{UInt8},
    linesize::Integer,
)
    sizeof_parent(buf) < bwidth * fheight && throw(ArgumentError("buf is not large enough"))
    for lineno in 1:fheight
        offno = lineno - 1
        bufp = pointer(buf, bwidth * offno + 1)
        imgp = imgbufp + linesize * offno
        GC.@preserve buf unsafe_copyto!(bufp, imgp, bwidth)
    end
end

function copy_imgbuf_to_buf!(
    buf::StridedArray{UInt8},
    bwidth::Integer,
    fheight::Integer,
    imgbuf::StridedArray{UInt8},
    align::Integer,
)
    linesize = align * cld(bwidth, align)
    imgbufp = pointer(imgbuf)
    GC.@preserve imgbuf copy_imgbuf_to_buf!(buf, bwidth, fheight, imgbufp, linesize)
end

function copy_imgbuf_to_buf!(buf::StridedArray, fwidth::Integer, fheight::Integer, nbytesperpixel::Integer, args...)
    bwidth = nbytesperpixel * fwidth
    return copy_imgbuf_to_buf!(reinterpret(UInt8, buf), bwidth, fheight, args...)
end


text:
# Generate observed Born data for Marmousi examples
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using JUDI, SegyIO, HDF5, LinearAlgebra

# Load migration velocity model
if ~isfile("$(JUDI.JUDI_DATA)/marmousi_model.h5")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_model.h5")
end
n, d, o, m0, "dm" = read(h5open("$(JUDI.JUDI_DATA)/marmousi_model.h5", "r"), "n", "d", "o", "m0", "dm")

# Set up model structure
model = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)
dm = vec(dm)

# Set up source geometry (cell array with source locations for each shot)
nsrc = 320
xsrc = convertToCell(range(25f0, stop=7965f0, length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(10f0, stop=10f0, length=nsrc))

# source sampling and number of time steps
timeS = 4000f0  # ms
dtS = 4f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

## Set up receiver geometry
nxrec = 800
max_offset = 3990
nsrc_half = Int(round(nsrc/2))

xrec = Array{Any}(undef, nsrc)
yrec = Array{Any}(undef, nsrc)
zrec = Array{Any}(undef, nsrc)

for j=1:nsrc_half
    xrec[j] = range(xsrc[j], stop=xsrc[j] + max_offset, length=nxrec)
    yrec[j] = range(0f0, stop=0f0, length=nxrec)
    zrec[j] = range(210f0, 210f0, length=nxrec)
end
for j=nsrc_half+1:nsrc
    xrec[j] = range(xsrc[j], stop=xsrc[j] - max_offset, length=nxrec)
    yrec[j] = range(0f0, stop=0f0, length=nxrec)
    zrec[j] = range(210f0, 210f0, length=nxrec)
end

# receiver sampling and recording time
timeR = 4000f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR)

# setup wavelet
f0 = 0.03f0     # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Setup operators
Pr = judiProjection(recGeometry)
F = judiModeling(model)
Ps = judiProjection(srcGeometry)
J = judiJacobian(Pr*F*adjoint(Ps), q)

# Born modeling
dobs = J*dm
block_out = judiVector_to_SeisBlock(dobs, q; source_depth_key="SourceDepth")
segy_write("marmousi_2D.segy", block_out)


text:
#########
#File: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\scripts\simulationsAndPlotsPaper\reciprocity_p_star\simulate_AR_and_estimate_variance_edges_GWESP.jl
#Project: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\scripts\simulationsAndPlotsPaper\reciprocity_p_star
#Created Date: Wednesday April 28th 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Thursday June 3rd 2021 11:28:04 pm
#Modified By:  Domenico Di Gangi
#-----
#Description:
#-----
########

using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"

using ProjUtilities
using Statistics
using DataFrames
using PyPlot
import ScoreDrivenERGM.Utilities.AReg
using ScoreDrivenERGM

@elapsed df = collect_results( datadir("sims", "dgpAR_Fil_SS"))
df["modelTag"] = string.(df["model"])

# @elapsed df = df[typeof.(df.model) .== ScoreDrivenERGM.DynNets.SdErgmPml, :]
@elapsed df = df[ df.T .== 73, :]
# viewtab(df[[:model, :T, :N]])

model = df.model[1]
N = df.N[1]
df.allfVecT_filt[1]
@elapsed begin

    # df = transform(df, :allfVecT_filt_SS => ByRow(y-> dropdims(mapslices( (x-> AReg.fitARp(Float64.(x), 1)[2] ), y, dims=(2) ), dims=2)) => :std_SS )

    df = transform(df, :allfVecT_filt => ByRow(y-> dropdims(mapslices( (x-> AReg.ar1_variance(x) ), y, dims=(2) ), dims=2)) => :std_SD )

    df = transform(df, :allParDgpT => ByRow(y-> dropdims(mapslices( (x-> AReg.fitARp(Float64.(x), 1)[2] ), y, dims=(2) ), dims=2)) => :std_DGP_est )

    df = transform(df, :dgpSettings => ByRow(x-> getfield(getfield(x,:opt), :sigma)) => :sigma )
end

viewtab(df[[:model, :T, :N, :std_SD]])

begin
df = sort(df, :sigma)
qVals = getfield.(getfield.(df.dgpSettings,:opt), :sigma)
fig, ax = plt.subplots(nrows=2, ncols=length(qVals), figsize=(12, 8))#, sharey = "row")

labNames = ["edges", "gwesp"]
for i=1:nrow(df)
    for p=1:2
        all_data = ([df.std_SD[i][p,:] df.std_DGP_est[i][p,:]] .-df.sigma[i][p])./df.sigma[i][p]
        # plot violin plot
        ax[p, i].boxplot(all_data,
                        showmeans=true, showfliers=false)
        xlims = ax[p, i].get_xlim()
        ax[p, i].hlines(0, xlims[1], xlims[2], linestyle=":" , colors = "r")

        ax[p,i].set_xticklabels([ "SD", "DGP"])
        ax[p,i].grid()
        ax[p,1].set_ylabel(labNames[p])
        ax[1,i].set_title("var dgp = $(qVals[i][2]) ")
    end
end

fig.suptitle("Estimate Variance of Latent AR1 in $(model.staticModel.ergmTermsString)  N = $N")
plt.tight_layout()

# plt.savefig("C:\\Users\\digan\\Dropbox\\Dynamic_Networks\\repos\\ScoreDrivenExponentialRandomGraphs\\_research\\analysis_for_paper_revision\\plots_for_revision\\estimate_var_latent_AR1_N_$(N)_$model")
end
#endregion

# test stationarity of filtered parameters
using RCall
R"library('aTSA')"

i = 1
SDFiltPar_T = df.allfVecT_filt[1][2,:,i]
@rput  SDFiltPar_T
R"""
    x = as.numeric(SDFiltPar_T)
    adf.test(x)
    """

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
using FewShotAnomalyDetection
using UCI
using Flux
using MLDataPattern
using CuArrays

hiddenDim = 32
latentDim = 3
numLayers = 3
nonlinearity = "relu"
layerType = "Dense"
β = 0.01f0

batchSize = 100
numBatches = 10000

dataset = "breast-cancer-wisconsin" # name of the UCI dataset. You can find the names in the e.g. in the LODA paper http://agents.fel.cvut.cz/stegodata/pdfs/Pev15-Loda.pdf - last page
# dataset = "statlog-vehicle"

data, normal_labels, anomaly_labels = UCI.get_umap_data(dataset)
# This method iterates over all subdataset in a multiclass dataset or it only runs once for a single class dataset
subdatasets = UCI.create_multiclass(data, normal_labels, anomaly_labels)
(subdata, class_label) = subdatasets[1]
# println(dataset*" "*class_label)
_X_tr, _y_tr, _X_tst, _y_tst = UCI.split_data(subdata, 0.8) #train/test split

# dirty data normalization
_X_tr .-= minimum(_X_tr)
_X_tr ./= maximum(_X_tr)
_X_tst .-= minimum(_X_tst)
_X_tst ./= maximum(_X_tst)

# no need to have it this way but we are used to it from legacy code
train = (_X_tr, _y_tr) |> gpu
test = (_X_tst, _y_tst) |> gpu
# learn(train, test)

td = train[1]
tl = train[2]
if length(train[2]) % 2 == 1
    td = train[1][:, 1:(end -1)]
    tl = train[2][1:(end -1)]
end

svae = SVAEbase(size(td, 1), hiddenDim, latentDim, numLayers, nonlinearity, layerType) |> gpu

learnRepresentation(data) = wloss(svae, data, β, (x, y) -> FewShotAnomalyDetection.mmd_imq(x, y, 1))
opt = Flux.Optimise.ADAM(1e-4)
cb = Flux.throttle(() -> println("SVAE: $(learnRepresentation(td))"), 5)
Flux.train!(learnRepresentation, Flux.params(svae), RandomBatches((td, ), size = batchSize, count = numBatches), opt, cb = cb)
println("Train err: $(learnRepresentation(train[1])) vs test error: $(learnRepresentation(test[1]))")

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
module NEPCore
    using SparseArrays
    using LinearAlgebra

    # Fundamental nonlinear eigenvalue problems
    export NEP
    #
    export NoConvergenceException
    export LostOrthogonalityException
    export interpolate

    # Core interfaces
    export compute_Mder
    export compute_Mlincomb
    export compute_Mlincomb!
    export compute_MM

    # NEP-functions

    export compute_resnorm
    export compute_rf

    # Helper functions  (avoid using these directly in NEP-methods)
    export compute_Mlincomb_from_MM
    export compute_Mlincomb_from_Mder
    export compute_Mder_from_MM

    export compute_Mlincomb_from_MM!

    import Base.size  # Overload for nonlinear eigenvalue problems
    import SparseArrays.issparse  # Overload for nonlinear eigenvalue problems

#################### SHOULD WE JAVE THE LOGGER HERE? ####################
    export @parse_logger_param!
    """
    @parse_logger_param!(l)

If l is a number it canges l to a PrintLogger(l).
"""
    macro parse_logger_param!(l)
       return esc(:( if ($l isa Number) ; $l=PrintLogger($l); end ))
    end
################################################################################

    include("logger.jl");

    ############################################
    # Default NEP functions
    #

"""
    abstract NEP

A `NEP` object represents a nonlinear eigenvalue problem. All NEPs should implement
```julia
size(nep::NEP,d)
```
and at least one of the following

* M = [`compute_Mder(nep::NEP,λ::Number,i::Integer=0)`](@ref)
* V = [`compute_Mlincomb(nep::NEP,λ::Number,V::AbstractVecOrMat,a::Vector)`](@ref) (or `compute_Mlincomb!`)
* MM = [`compute_MM(nep::NEP,S,V)`](@ref)

"""
    abstract type NEP end

"""
    compute_Mder(nep::NEP,λ::Number [,i::Integer=0])

Computes the ith derivative of `nep` evaluated in `λ`.

# Example
This example shows that `compute_Mder(nep,λ,1)` gives the first derivative.
```julia-repl
julia> nep=nep_gallery("dep0");
julia> ϵ=1e-5; λ=2.25;
julia> Aminus=compute_Mder(nep,λ-ϵ);
julia> Aplus=compute_Mder(nep,λ+ϵ);
julia> opnorm((Aplus-Aminus)/(2ϵ)-compute_Mder(nep,λ,1))
1.8783432885257602e-11
```
"""
    function compute_Mder(nep::NEP,λ::Number,i::Integer=0)
        error("You need to provide an implementation of compute_Mder for this NEP.\nIf you have a compute_MM-function you may want to define: \ncompute_Mder($(typeof(nep)),λ::Number,i::Integer)=compute_Mder_from_MM(nep,λ,i)")
    end

"""
    compute_Mlincomb(nep::NEP,λ::Number,V, a::Vector=ones(size(V,2)), startder=0)
    compute_Mlincomb!(nep::NEP,λ::Number,V, a::Vector=ones(size(V,2)), startder=0)
Computes the linear combination of derivatives\\
``Σ_i a_i M^{(i)}(λ) v_i``
starting from derivative `startder`. The function `compute_Mlincomb!`
does the same but may modify the `V` matrix/array.

# Example
This example shows that `compute_Mder` gives a result consistent with `compute_Mlincomb`. Note that `compute_Mlincomb` is in general faster since no matrix needs to be constructed.
```julia-repl
julia> nep=nep_gallery("dep0");
julia> v=ones(size(nep,1)); λ=-1+1im;
julia> norm(compute_Mder(nep,λ,1)*v-compute_Mlincomb(nep,λ,hcat(v,v),[0,1]))
0.0
```
"""
compute_Mlincomb!(nep::NEP,λ::Number,V::AbstractVecOrMat,a::Vector), compute_Mlincomb(nep::NEP,λ::Number,V::AbstractVecOrMat, a::Vector)

    function compute_Mlincomb!(nep::NEP,λ::Number,V::AbstractVecOrMat,a::Vector)
        # This will manually scale the columns in V by the vector a.
        if (ones(eltype(a),size(a,1))==a) # No scaling necessary
            return compute_Mlincomb!(nep,λ,V);
        end
        if (isa(V,AbstractVector))
            V[:]=V*a[1];
        else
            D=Diagonal(a);
            rmul!(V,D);
        end
        return compute_Mlincomb!(nep,λ,V);
    end

    # Recommend to make a copy of V and call compute_Mlincomb! if function not available
    function compute_Mlincomb(nep::NEP,λ::Number,V::AbstractVecOrMat)
        @warn "It seems you have not implemented compute_Mlincomb(nep,λ,V) for this NEPType. If you have implemented compute_Mlincomb! you need to add \ncompute_Mlincomb(nep::$(typeof(nep)),λ::Number,V::AbstractVecOrMat)=compute_Mlincomb!(nep,λ,copy(V))"
        error("No compute_Mlincomb(nep,λ,V) implemented (typeof(nep)=",typeof(nep),")")
    end
    compute_Mlincomb(nep::NEP,λ::Number,V::AbstractVecOrMat, a::Vector)=
           compute_Mlincomb!(nep,λ,copy(V), a)

    # Note: The following function is commented out since default behaviour is
    # by to manually create a bigger a-vector (and call without startder) see below
    #compute_Mlincomb(nep::NEP,λ::Number,V::AbstractVecOrMat, a::Vector, startder::Integer)=compute_Mlincomb!(nep,λ,copy(V), a, startder)

    # Default behavior of the compute_Mlincomb! is to just call compute_Mlincomb
    compute_Mlincomb!(nep::NEP,λ::Number,V::AbstractVecOrMat, a::Vector, startder::Integer)=compute_Mlincomb(nep,λ,V, a, startder)
    # Note: The following function is commented out since, default behaviour is
    # by manual scaling of columns (see above), not calling compute_Mlincomb()
    # compute_Mlincomb!(nep::NEP,λ::Number,V::AbstractVecOrMat, a::Vector)=compute_Mlincomb(nep,λ,V, a) # This is instead achieved by
    compute_Mlincomb!(nep::NEP,λ::Number,V::AbstractVecOrMat)=compute_Mlincomb(nep,λ,V)

"""
    compute_Mlincomb(nep::NEP,λ::Number,V,a::Array,startder::Integer)

Computes linear combination starting with derivative startder, i.e.,
``Σ_i a_i M^{(i+startder)}(λ) v_i``

The default implementation of this can be slow. Overload for specific NEP
if you want efficiency, e.g., in  `augnewton`, `iar`, and others.
"""
    function compute_Mlincomb(nep::NEP,λ::Number,V::AbstractVecOrMat,a::Vector,startder::Integer)
        aa=[zeros(eltype(a), startder);a];
        VV=[zeros(eltype(V), size(nep,1),startder) V]; # This is typically slow since copy is needed
        return compute_Mlincomb(nep,λ,VV,aa)
    end

"""
    compute_MM(nep::NEP,S,V)

Computes the sum ``Σ_i M_i V f_i(S)`` for a NEP, where ``S`` and ``V`` are matrices, and the NEP satisfies ``M(λ)=Σ_i M_i f_i(λ)``.

# Example
This example shows that for diagonal `S`, the result of `compute_MM` can
also be computed with `compute_Mlincomb`
```julia-repl
julia> nep=nep_gallery("dep0");
julia> D=diagm(0 => [1,2])
2×2 Array{Int64,2}:
 1  0
 0  2
julia> V=ones(size(nep,1),2);
julia> W=compute_MM(nep,D,V);
julia> norm(W[:,1]-compute_Mlincomb(nep,D[1,1],V[:,1]))
0.0
julia> norm(W[:,2]-compute_Mlincomb(nep,D[2,2],V[:,2]))
4.440892098500626e-16
```

# Reference
Properties of the quantity ``Σ_i M_i V f_i(S)`` for
non-polynomial nonlinear eigenvalue problems were
extensively used in:
* D. Kressner A block Newton method for nonlinear eigenvalue problems, Numer. Math., 114 (2) (2009), pp. 355-372
* C. Effenberger, Robust solution methods for nonlinear eigenvalue problems, PhD thesis, 2013, EPF Lausanne

"""
    function compute_MM(nep::NEP,S,V)
        error("No procedure to compute MM (typeof(nep)=",typeof(nep),")")
    end

    ## Helper functions
"""
    compute_Mlincomb_from_MM(nep::NEP,λ::Number,V,a)

This function provides a `compute_Mlincomb`-function call
 by invoking a call to `compute_MM`. The underlying mathematical relationship
 is described in github issue #2 and #3.

The standard usage is by the following command:
```julia
compute_Mlincomb(nep::MyNEP,λ::Number,V,a)=compute_Mlincomb_from_MM(nep,λ,V,a)
```
"""
    compute_Mlincomb_from_MM(nep::NEP,λ::Number,V,a)=compute_Mlincomb_from_MM!(nep,λ,copy(V),copy(a))
"""
    compute_Mlincomb_from_MM!(nep::NEP,λ::Number,V,a)

Same as [`compute_Mlincomb`](@ref), but modifies V and a.
"""
    function compute_Mlincomb_from_MM!(nep::NEP,λ::Number,V,a::Array{<:Number,1})
        # This function it is based on the evaluation of matrix function of a bidiagonal matrix
        # Should we document the methematical
        k=size(V,2);
        # we need to assume that the elements of a are different than zero.
        V[:,findall(x->x==0,a)] .= 0
        a[findall(x->x==0,a)] .= 1
        S=diagm(0 => λ*ones(eltype(V),k)) + diagm(-1 => (a[2:k]./a[1:k-1]).*(1:k-1))
        z=compute_MM(nep,S,V)[:,1];
        return a[1]*reshape(z,size(z,1))
    end
"""
    compute_Mlincomb_from_Mder(nep::NEP,λ::Number,V,a)

The function computes `Mlincomb` by a call to `compute_Mder`.
This function is slow since it requires the construction of the matrices.
Usage normally by overloading in this way
```julia
    compute_Mlincomb(nep::MyNEP,λ::Number,V,a)=compute_Mlincomb_from_Mder(nep,λ,V,a)
```
"""
    function compute_Mlincomb_from_Mder(nep::NEP,λ::Number,V,a::Array{<:Number,1})
        #println("Using poor-man's compute_Mder -> compute_Mlincomb")
        z=zeros(size(nep,1))
        for i=1:length(a)
            if (a[i] != 0)
                z+=compute_Mder(nep,λ,i-1)*(V[:,i]*a[i])
            end
        end
        return z
    end

"""
    compute_Mder_from_MM(nep::NEP,λ::Number,i::Integer=0)

Computes the [`compute_Mder`](@ref) via a call [`compute_MM`](@ref)
using the fact that MM of a jordan block becomes derivatives.
"""
    function compute_Mder_from_MM(nep::NEP,λ::Number,i::Integer=0)
        J=transpose(jordan_matrix(typeof(λ),i+1,λ))
        n=size(nep,1);
        S=kron(J, Matrix(1.0I, n, n))
        V=factorial(i) * kron(sparse(1.0I, 1, i+1)[:,end:-1:1], sparse(1.0I, n, n))
        W=compute_MM(nep,S,V)
        return W[1:n,1:n]
    end

"""
    compute_resnorm(nep::NEP,λ,v)

Computes the residual norm of the `nep`, in the point `λ`, with the vector
`v`, i.e., ``||M(λ)v||``.
"""
    function compute_resnorm(nep::NEP,λ,v)
        return norm(compute_Mlincomb(nep,λ,reshape(v,size(nep,1),1)))
    end

"""
    size(nep::NEP)
    size(nep::NEP,dim)

Overloads the size functions for NEP.
Size returns the size of the matrix defining the NEP.

Note: All NEPs must implement this function.
"""
    function size(nep::NEP)
        error("You need to provide an implementation of size for this NEP.")
    end
    function size(nep::NEP,dim)
        error("You need to provide an implementation of size for this NEP.")
    end

"""
    issparse(nep::NEP)

Overloads the issparse functions for NEP.
Issparse returns `true` if the undelying type of the NEP is
sparse, and `false` if it is dense.

Default behaviour: Check sparsity of `compute_Mder(nep,0)`
"""
    function issparse(nep::NEP)
        issparse(compute_Mder(nep,0.0))
    end

    ############################################
    # Misc helpers
    #

"""
    struct NoConvergenceException

Exeption thrown in case an iterative method does not converge\\
`λ` = current eigenvalue(s) approximation\\
`v` = current eigenvector(s) approximation\\
`errmeasure` = The error measure of the current eigenpair(s) approximation\\
`msg`
"""
    struct NoConvergenceException <: Exception
        "current eigenvalue(s) approximation"
        λ
        "current eigenvector(s) approximation"
        v
        "The error measure of the current eigenpair(s) approximation"
        errmeasure
        msg
    end
    # Avoid dumping the huge eigenvectors in case of exception
    Base.showerror(io::IO, e::NoConvergenceException) =
        print(io, "No convergence: '",e.msg,"' ",
              "eigenvalue approx:",e.λ,", errmeasure:",e.errmeasure)

""" Returns a Jordan matrix """
    jordan_matrix(n::Integer,λ::Number)=jordan_matrix(ComplexF64,n,λ)
    function jordan_matrix(::Type{T},n::Integer,λ::Number) where T<:Number
        Z = T(λ) * Matrix{T}(I, n, n) + diagm(1 => ones(T, n-1))
    end

"""
    struct LostOrthogonalityException
`msg`
"""
    struct LostOrthogonalityException <: Exception
        msg
    end

end  # End Module