text:
using Profile

include("../src/MDDF.jl")

# Here we use the PDBTools package to read the pdb file (from http://github.com/m3g/PDBTools)
using PDBTools
atoms = PDBTools.readPDB("./structure.pdb")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
solute_indexes = [ atom.index for atom in filter( atom -> PDBTools.isprotein(atom), atoms ) ]
solute = MDDF.Selection( solute_indexes, nmols=1 )

# Input options for the calcualtion
options = MDDF.Options(binstep=0.2,lastframe=1,lcell=2)

Profile.clear()

for num in [ 10000000 ]

  println(" n = ", num)

  solvent_indexes = [ atom.index for atom in filter( atom -> (atom.resname == "TIP3" && atom.resnum <= num ), atoms ) ]
  solvent = MDDF.Selection( solvent_indexes, natomspermol=3 )

  # Initialize trajectroy data structure and open input stream

  # Run MDDF calculation, and get the resutls in the R structure
  println("linkedcells:")
  trajectory = MDDF.NamdDCD("./trajectory.dcd",solute,solvent)
  R = MDDF.mddf_linkedcells(trajectory,options)
  trajectory = MDDF.NamdDCD("./trajectory.dcd",solute,solvent)
  @time R = MDDF.mddf_linkedcells(trajectory,options)

  #println("naive:")
  trajectory = MDDF.NamdDCD("./trajectory.dcd",solute,solvent)
  R = MDDF.mddf_naive(trajectory,options)
  trajectory = MDDF.NamdDCD("./trajectory.dcd",solute,solvent)
  @time R = MDDF.mddf_naive(trajectory,options)

  println(" profiling linkedcells... ")
  trajectory = MDDF.NamdDCD("./trajectory.dcd",solute,solvent)
  @profile R = MDDF.mddf_linkedcells(trajectory,options)

  #trajectory = MDDF.NamdDCD("./trajectory.dcd",solute,solvent)
  #println(" profiling naive... ")
  #profile R = MDDF.mddf_naive(trajectory,options)

end


text:
"""
    ⊠(V₁::VectorSpace, V₂::VectorSpace)

Given two vector spaces `V₁` and `V₂` (`ElementarySpace` or `ProductSpace`), or thus,
objects of corresponding fusion categories ``C₁`` and ``C₂``, ``V₁ ⊠ V₂`` constructs the
Deligne tensor product, an object in ``C₁ ⊠ C₂`` which is the natural tensor product of
those categories. In particular, the corresponding type of sectors (simple objects) is
given by `sectortype(V₁ ⊠ V₂) == sectortype(V₁) ⊠ sectortype(V₂)` and can be thought of as
a tuple of the individual sectors.

The Deligne tensor product also works in the type domain and for sectors and tensors. For
group representations, we have `Rep[G₁] ⊠ Rep[G₂] == Rep[G₁ × G₂]`, i.e. these are the
natural representation spaces of the direct product of two groups.
"""
⊠(V1::VectorSpace, V2::VectorSpace) = (V1 ⊠ one(V2)) ⊗ (one(V1) ⊠ V2)

# define deligne products with empty tensor product: just add a trivial sector of the type of the empty space to each of the sectors in the non-empty space
function ⊠(V::GradedSpace, P0::ProductSpace{<:EuclideanSpace{ℂ},0})
    I1 = sectortype(V)
    I2 = sectortype(P0)
    return Vect[I1 ⊠ I2](ifelse(isdual(V), dual(c), c) ⊠ one(I2) => dim(V, c) for c in sectors(V); dual = isdual(V))
end

function ⊠(P0::ProductSpace{<:EuclideanSpace{ℂ},0}, V::GradedSpace)
    I1 = sectortype(P0)
    I2 = sectortype(V)
    return Vect[I1 ⊠ I2](one(I1) ⊠ ifelse(isdual(V), dual(c), c) => dim(V, c) for c in sectors(V); dual = isdual(V))
end

function ⊠(V::ComplexSpace, P0::ProductSpace{<:EuclideanSpace{ℂ},0})
    I2 = sectortype(P0)
    return Vect[I2](one(I2) => dim(V); dual = isdual(V))
end

function ⊠(P0::ProductSpace{<:EuclideanSpace{ℂ},0}, V::ComplexSpace)
    I1 = sectortype(P0)
    return Vect[I1](one(I1) => dim(V); dual = isdual(V))
end

function ⊠(P::ProductSpace{<:EuclideanSpace{ℂ},0}, P0::ProductSpace{<:EuclideanSpace{ℂ},0})
    I1 = sectortype(P)
    I2 = sectortype(P0)
    return one(Vect[I1 ⊠ I2])
end

function ⊠(P::ProductSpace{<:EuclideanSpace{ℂ}}, P0::ProductSpace{<:EuclideanSpace{ℂ},0})
    I1 = sectortype(P)
    I2 = sectortype(P0)
    S = Vect[I1 ⊠ I2]
    N = length(P)
    return ProductSpace{S,N}(map(V->V ⊠ P0, tuple(P...)))
end

function ⊠(P0::ProductSpace{<:EuclideanSpace{ℂ},0}, P::ProductSpace{<:EuclideanSpace{ℂ}})
    I1 = sectortype(P0)
    I2 = sectortype(P)
    S = Vect[I1 ⊠ I2]
    N = length(P)
    return ProductSpace{S,N}(map(V->P0 ⊠ V, tuple(P...)))
end

ext: jl


text:
#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

export iterlines, itertweets, loadtweets
import GZip
import JSON

"""
   iterlines(fun, filename; maxlines=typemax(Int))

Calls `fun` over each line of the given filename, skipping empty lines. It supports gzip compressed inputs.

"""
function iterlines(fun, filename::String; maxlines=typemax(Int))
    if endswith(filename, ".gz")
        f = GZip.open(filename)
        i = 0
        while !eof(f) && i < maxlines
            line = readline(f)
            i += 1
            if length(line) == 0
                continue
            end
            fun(line)
        end
        close(f)
    else
        open(filename) do f
            i = 0
            while !eof(f) && i < maxlines
                line = readline(f)
                i += 1
                fun(line)
            end
        end
    end
end

"""
   parsetweet(line) -> Dict{String,Any}

Parses a json-dictionary string, which is the assumed object format.
It supports two formats:

- a raw json-dictionaries
- a dictionary prefixed by some prefix key, that is, key<tab>json-dictionary

After parsing the object, the `key` is inserted into the dictionary with the keyword "key"
`dic["key"] = key`

"""
function parsetweet(line)
    if line[1] == '{'
        tweet = JSON.parse(line)
    else
        key, value = split(line, '\t', limit=2)
        tweet = JSON.parse(value)
        tweet["key"] = key
    end

    tweet
end

"""
   itertweets(fun, filename::String; maxlines=typemax(Int))
   itertweets(fun, file; maxlines=typemax(Int))

Calls `fun` for each object in the file. One object per line using
the format supported by `parsetweet`

"""
function itertweets(fun, filename::String; maxlines=typemax(Int))
    iterlines(filename, maxlines=maxlines) do line
        tweet = parsetweet(line)
        fun(tweet)
    end
end

function itertweets(fun, file; maxlines=typemax(Int))
    i = 0
    while !eof(file) && i < maxlines
        line = readline(file)
        i += 1
        try
            tweet = parsetweet(line)
            fun(tweet)
        catch
            continue
        end
    end
end

"""
    loadtweets(f; maxlines=typemax(Int))

loads a list of tweets from a file using `itertweets`
"""
function loadtweets(f; maxlines=typemax(Int))
    L = Dict[]
    itertweets(f, maxlines=maxlines) do t
        push!(L, t)
    end
    L
end


text:
# This file was generated, do not modify it. # hide
μ = randn(1000)
σ = rand(1000)
x = randn(1000)
y = randn(1000)

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
struct StateEstimationDC
    main::Array{Float64,2}
    flow::Array{Float64,2}
    estimate::Array{Float64,2}
    error::Array{Float64,1}
    baddata::Array{Float64,2}
    observability::Array{Array{Int64,1},1}
end

#### DC state estimation
function rundcse(system, measurements, num, numsys, settings, info)
    printstyled("Algorithm: DC state estimation\n"; bold = true)

    ########## Pre-processing ##########
    busi, type, Gshunt, Tini, branchi, reactance, transTap, transShift, branchOn,
    meanPij, variPij, onPij, meanPi, variPi, onPi,
    meanTi, variTi, onTi  = read_dcsystem(system, measurements)

    main = fill(0.0, numsys.Nbus, 3)
    flow = fill(0.0, numsys.Nbranch, 5)
    Va, Pinj, Pij, Pji, badsave, islands = write_dcsystem(system, numsys, main, flow)

  algtime = @elapsed begin
    ########## Convert in integers ##########
    fromi = convert(Array{Int64,1}, system.branch[:, 2])
    toi = convert(Array{Int64,1}, system.branch[:, 3])
    fromiPij = convert(Array{Int64,1}, measurements.legacyFlow[:, 2])
    toiPij = convert(Array{Int64,1}, measurements.legacyFlow[:, 3])
    busiPi = convert(Array{Int64,1}, measurements.legacyInjection[:, 1])
    busiTi = convert(Array{Int64,1}, measurements.pmuVoltage[:, 1])
    branchiPij = convert(Array{Int64,1}, measurements.legacyFlow[:, 1])

    ########## Numbering and slack bus ##########
    bus = collect(1:numsys.Nbus)
    branch = collect(1:numsys.Nbranch)
    newnumbering = false
    Tslack = 0.0
    slack = 0
    @inbounds for i = 1:numsys.Nbus
        if bus[i] != busi[i]
            newnumbering = true
        end
        if type[i] == 3
            Tslack = Tini[i]
            slack = bus[i]
        end
    end
    if newnumbering
        println("The new bus numbering is running.")
    end
    if slack == 0
        slack = 1
        println("The slack bus is not found. Slack bus is the first bus.")
    end

    numlabel = false
    @inbounds for i = 1:numsys.Nbranch
        if branch[i] != branchi[i]
            numlabel = true
            break
        end
    end
    if numlabel
        println("The new branch numbering is running.")
    end
    from = renumber(fromi, busi, bus, newnumbering)
    to = renumber(toi, busi, bus, newnumbering)
    fromPij = renumber(fromiPij, busi, bus, newnumbering)
    toPij = renumber(toiPij, busi, bus, newnumbering)
    busPi = renumber(busiPi, busi, bus, newnumbering)
    busTi = renumber(busiTi, busi, bus, newnumbering)
    branchPij = renumber(branchiPij, branchi, branch, newnumbering)

    ########## Ybus matrix ##########
    Ybus, admitance, Pshift = ybusdc(system, numsys, bus, from, to, branchOn, transTap, reactance, transShift)

    ########## Number of Jacobian elements and measurements ##########
    Nelement = 0; Nmeasure = 0; Nflow = 0; Nvol = 0; Ninj = 0
    @inbounds for i = 1:num.legacyNf
        if onPij[i] == 1
            Nelement += 2; Nmeasure += 1; Nflow += 1
        end
    end
    @inbounds for i = 1:num.legacyNi
        if onPi[i] == 1
            Nmeasure += 1; Ninj += 1
            for j in Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
                Nelement += 1
            end
        end
    end
    @inbounds for i = 1:num.pmuNv
        if onTi[i] == 1
            Nelement += 1; Nmeasure += 1; Nvol += 1
        end
    end

    ########## Jacobian and measurement data ##########
    row = fill(0, Nelement); col = similar(row); jac = fill(0.0, Nelement)
    mean = fill(0.0, Nmeasure); weight = similar(mean)
    index = 1; rowindex = 1
    @inbounds for i = 1:num.legacyNf
        if onPij[i] == 1
            row[index] = rowindex
            col[index] = fromPij[i]
            jac[index] = admitance[branchPij[i]]
            index += 1
            row[index] = rowindex
            col[index] = toPij[i]
            jac[index] = -admitance[branchPij[i]]

            if fromPij[i] == from[branchPij[i]]
                mean[rowindex] = meanPij[i] + transShift[branchPij[i]] * (pi / 180)  * admitance[branchPij[i]]
            else
                mean[rowindex] = meanPij[i] - transShift[branchPij[i]] * (pi / 180) * admitance[branchPij[i]]
            end
            weight[rowindex] = 1 / variPij[i]

            index += 1; rowindex += 1
        end
    end

    if settings.observe[:observe] == 1
        N = 2 * Nflow
        Jflow = sparse(row[1:N], col[1:N], jac[1:N], Nflow, numsys.Nbus)
    end

    cnt = 1
    @inbounds for i = 1:num.legacyNi
        if onPi[i] == 1
            for j in Ybus.colptr[busPi[i]]:(Ybus.colptr[busPi[i] + 1] - 1)
                row[index] = rowindex
                col[index] = Ybus.rowval[j]
                jac[index] = Ybus.nzval[j]
                index += 1
            end

            mean[rowindex] = meanPi[i] - Pshift[i] - Gshunt[i] / system.basePower
            weight[rowindex] = 1 / variPi[i]

            rowindex += 1; cnt += 1
        end
    end

    cnt = 1
    @inbounds for i = 1:num.pmuNv
        if onTi[i] == 1
            row[index] = rowindex
            col[index] = busTi[i]
            jac[index] = 1.0

            mean[rowindex] = meanTi[i] - Tslack * (pi / 180)
            weight[rowindex] = 1 / variTi[i]

            index += 1; rowindex += 1; cnt += 1
        end
    end
    J = sparse(row, col, jac, Nmeasure, numsys.Nbus)

    ########## Observability analysis ##########
    Npseudo = 0
    if settings.observe[:observe] == 1
        r = rank(J)
        if r < numsys.Nbus
            println("The power system is unobservable, the column rank of the Jacobian matrix is $r." )

            islands, Nislands = islandsflow(Jflow)
            if settings.observe[:islands] == 1
                islands, Nisland = islandstopological(settings, numsys, J, Jflow, Ybus, branch, from, to, busPi, onPi, islands)
            elseif settings.observe[:islands] == 2
                islands, Nisland = islandsbp(settings, numsys, J, bus, islands, Nislands)
            end

            if settings.observe[:restore] == 1
                J, mean, weight, Npseudo = restorationgram(system, settings, numsys, num, measurements, islands, Nisland, J, mean, weight,
                    Ybus, slack, Tslack, branch, from, to, fromPij, toPij, busPi, branchPij, busTi, onPi, onPij, onTi, Nvol, Nflow,
                    meanPij, meanPi, meanTi, variPi, variTi, variPij, admitance, Gshunt, transShift, Pshift)
            elseif settings.observe[:restore] == 2
                J, mean, weight, Npseudo = restorationbp(system, settings, numsys, num, measurements, islands, Nisland, J, mean, weight,
                    Ybus, branch, from, to, busPi, onPi, meanPi, variPi, Pshift, Gshunt)
            end
        else
            println("The power system is observable.")
        end
    end
    Nmeasure = Nmeasure + Npseudo

    ######### Remove column of the slack bus, but keep full column rank ##########
    for i in J.colptr[slack]:(J.colptr[slack + 1] - 1)
        J.nzval[i] = 0.0
    end
    J = [J; sparse([1], [slack], [1.0], 1, numsys.Nbus)]
    push!(mean, 0.0); push!(weight, 1.0);

    ######### LAV ##########
    if settings.lav[:lav] == 1
        Ti = Array{Float64}(undef, numsys.Nbus)
        x = lav(J, mean, numsys.Nbus, Nmeasure + 1, settings)
        for i = 1:numsys.Nbus
            Ti[i] = x[i] - x[numsys.Nbus + i]
        end
    end

    ######## WLS and bad data analysis ##########
    if settings.lav[:lav] == 0
        W = spdiagm(0 => sqrt.(weight))

        if settings.bad[:bad] == 1
            idx = findall(!iszero, J)
            v = 1 ./ weight
        end

        pass = 1; rbelow = true
        while rbelow
            if settings.bad[:pass] < pass || settings.bad[:bad] == 0
                rbelow = false
            end

            H = W * J
            G = transpose(H) * H
            b = transpose(H) * W * mean
            Ti = wls(G, b, settings.solve)

            if rbelow
                J, mean, W, badsave, rbelow = baddata(settings, numsys, Ti, mean, v, W, G, J, Nmeasure, idx, badsave, rbelow)
                pass += 1
            end
        end
    end

    ########## Post-processing ##########
    Ti =  (pi / 180) * Tslack .+ Ti

    @inbounds for i = 1:numsys.Nbranch
        if branchOn[i] == 1
            Pij[i] = admitance[i] * (Ti[from[i]] - Ti[to[i]] - (pi / 180) * transShift[i]) * system.basePower
            Pji[i] = -Pij[i]
        end
    end

    Pinj[:] = (Ybus * Ti + Pshift) * system.basePower + Gshunt
    for i = 1:numsys.Nbus
        Va[i] = (180 / pi) * Ti[i]
    end

    ePmuVol = size(measurements.pmuVoltage, 2) == 9
    eLegFlo = size(measurements.legacyFlow, 2) == 11
    eLegInj = size(measurements.legacyInjection, 2) == 9
    estimate_row = 9; error_row = 3; exact = false
    if ePmuVol && eLegFlo && eLegInj
        estimate_row = 11; error_row = 6; exact = true
    end
    estimate = zeros(Nmeasure, estimate_row)
    error = zeros(error_row)

    Nm = Nmeasure - size(badsave, 1)
    idx = 1; scaleTi = 180 / pi
    for (i, on) in enumerate(onPij)
        if on != 0
            estimate[idx, 1:5] = [idx 1.0 1.0 1.0 i]
            if idx in badsave[:, 1]
                estimate[idx, 2] = 2.0; onPij[i] = 0
            end
            if on == 2
                estimate[idx, 2] = 3.0; onPij[i] = 1
            end

            estimate[idx, 6] = measurements.legacyFlow[i, 4] * system.basePower
            estimate[idx, 7] = measurements.legacyFlow[i, 5] * system.basePower
            k = branchPij[i]
            if measurements.legacyFlow[i, 2] == system.branch[k, 2] && measurements.legacyFlow[i, 3] == system.branch[k, 3]
                estimate[idx, 8] = flow[k, 4]
            else
                estimate[idx, 8] = flow[k, 5]
            end
            estimate[idx, 9] = abs(estimate[idx, 6] - estimate[idx, 8])

            if estimate[idx, 1] != 2.0
                error[1] += estimate[idx, 9] / (Nm * system.basePower)
                error[2] += estimate[idx, 9]^2 / (Nm * system.basePower^2)
                error[3] += estimate[idx, 9]^2 / (measurements.legacyFlow[i, 5] * system.basePower^2)
            end

            if exact
                estimate[idx, 10] = measurements.legacyFlow[i, 10] * system.basePower
                estimate[idx, 11] = abs(estimate[idx, 8] - estimate[idx, 10])

                if estimate[idx, 1] != 2.0
                    error[4] += estimate[idx, 11] / (Nm * system.basePower)
                    error[5] += estimate[idx, 11]^2 / (Nm * system.basePower^2)
                    error[6] += estimate[idx, 11]^2 / (measurements.legacyFlow[i, 5] * system.basePower^2)
                end
            end
            idx += 1
        end
    end
    for (i, on) in enumerate(onPi)
        if on != 0
            estimate[idx, 1:5] = [idx 1.0 1.0 3.0 i]
            if idx in badsave[:, 1]
                estimate[idx, 2] = 2.0; onPi[i] = 0
            end
            if on == 2
                estimate[idx, 2] = 3.0; onPi[i] = 1
            end

            estimate[idx, 6] = measurements.legacyInjection[i, 2] * system.basePower
            estimate[idx, 7] = measurements.legacyInjection[i, 3] * system.basePower
            estimate[idx, 8] = main[busPi[i], 3]
            estimate[idx, 9] = abs(estimate[idx, 6] - estimate[idx, 8])

            if estimate[idx, 1] != 2.0
                error[1] += estimate[idx, 9] / (Nm * system.basePower)
                error[2] += estimate[idx, 9]^2 / (Nm * system.basePower^2)
                error[3] += estimate[idx, 9]^2 / (measurements.legacyInjection[i, 3] * system.basePower^2)
            end

            if exact
                estimate[idx, 10] = measurements.legacyInjection[i, 8] * system.basePower
                estimate[idx, 11] = abs(estimate[idx, 8] - estimate[idx, 10])

                if estimate[idx, 1] != 2.0
                    error[4] += estimate[idx, 11] / (Nm * system.basePower)
                    error[5] += estimate[idx, 11]^2 / (Nm * system.basePower^2)
                    error[6] += estimate[idx, 11]^2 / (measurements.legacyInjection[i, 3] * system.basePower^2)
                end
            end
            idx += 1
        end
    end
    for (i, on) in enumerate(onTi)
        if on != 0
            estimate[idx, 1:5] = [idx 1.0 2.0 8.0 i]
            if idx in badsave[:, 1]
                estimate[idx, 2] = 2.0; onTi[i] = 0
            end
            if on == 2
                estimate[idx, 2] = 3.0; onTi[i] = 1
            end

            estimate[idx, 6] = measurements.pmuVoltage[i, 5] * scaleTi
            estimate[idx, 7] = measurements.pmuVoltage[i, 6] * scaleTi
            estimate[idx, 8] = main[busTi[i], 2]
            estimate[idx, 9] = abs(estimate[idx, 6] - estimate[idx, 8])

            if estimate[idx, 1] != 2.0
                error[1] += estimate[idx, 9] / (Nm * scaleTi)
                error[2] += estimate[idx, 9]^2 / (Nm * scaleTi^2)
                error[3] += estimate[idx, 9]^2 / (measurements.pmuVoltage[i, 6] * scaleTi^2)
            end

            if exact
                estimate[idx, 10] = measurements.pmuVoltage[i, 9] * scaleTi
                estimate[idx, 11] = abs(estimate[idx, 8] - estimate[idx, 10])

                if estimate[idx, 1] != 2.0
                    error[4] += estimate[idx, 11] / (Nm * scaleTi)
                    error[5] += estimate[idx, 11]^2 / (Nm * scaleTi^2)
                    error[6] += estimate[idx, 11]^2 / (measurements.pmuVoltage[i, 6] * scaleTi^2)
                end
            end
            idx += 1
        end
    end
    error[2] = sqrt(error[2])
    if exact
        error[5] = sqrt(error[5])
    end

    pass = collect(1:size(badsave, 1))
    idxbad = trunc.(Int, badsave[:, 1])
    bad = [pass estimate[idxbad, 3:5] badsave[:, 2] estimate[idxbad, 2]]

    if settings.observe[:observe] == 1 && Npseudo != 0 && newnumbering
        Nisland = size(islands, 1)
        for k = 1:Nisland
            for (n, i) in enumerate(k)
                islands[k][n] = system.bus[i, 1]
            end
        end
    end
 end # algtime

    ######### Results ##########
    results = StateEstimationDC(main, flow, estimate, error, bad, islands)
    header, group = results_estimatedc(system, numsys, measurements, num, settings, results, idxbad, Npseudo, algtime, slack)
    if !isempty(settings.save)
        savedata(results, measurements, system; info = info, group = group, header = header, path = settings.save)
    end

    return results
end

### Read system
function read_dcsystem(system, measurements)
    busi = @view(system.bus[:, 1])
    type = @view(system.bus[:, 2])
    Gshunt = @view(system.bus[:, 5])
    Tini = @view(system.bus[:, 9])

    branchi = @view(system.branch[:, 1])
    reactance = @view(system.branch[:, 5])
    transTap = @view(system.branch[:, 10])
    transShift = @view(system.branch[:, 11])
    branchOn = @view(system.branch[:, 12])

    meanPij = @view(measurements.legacyFlow[:, 4])
    variPij = @view(measurements.legacyFlow[:, 5])
    onPij = @view(measurements.legacyFlow[:, 6])

    meanPi = @view(measurements.legacyInjection[:, 2])
    variPi = @view(measurements.legacyInjection[:, 3])
    onPi = @view(measurements.legacyInjection[:, 4])

    meanTi = @view(measurements.pmuVoltage[:, 5])
    variTi = @view(measurements.pmuVoltage[:, 6])
    onTi = @view(measurements.pmuVoltage[:, 7])

    return busi, type, Gshunt, Tini,
           branchi, reactance, transTap, transShift, branchOn,
           meanPij, variPij, onPij,
           meanPi, variPi, onPi,
           meanTi, variTi, onTi
end

### Write data
function write_dcsystem(system, numsys, main, flow)
    for i = 1:numsys.Nbus
        main[i, 1] = system.bus[i, 1]
    end
    Va = @view(main[:, 2])
    Pinj = @view(main[:, 3])

    for i = 1:numsys.Nbranch
        flow[i, 1] = system.branch[i, 1]
        flow[i, 2] = system.branch[i, 2]
        flow[i, 3] = system.branch[i, 3]
    end
    Pij = @view(flow[:, 4])
    Pji = @view(flow[:, 5])

    badsave = reshape(Float64[], 0, 2)
    islands = [Int64[]]

    return Va, Pinj, Pij, Pji, badsave, islands
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
using Statistics, SparseArrays, LinearAlgebra

# mutable struct structureFactor
# 	fields
# end
#=
#Correlation Functions
=#
function longitudinalDensityCurrent(data, N; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates Jₓ(k,t) = m/N * Σ vₓⁱ(t) exp(-𝚤 k ⋅ rⁱ(t))
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then Jₓ(0,t) = m/2
    =#
    x, y, z, vx, vy, vz = data
    #Longidutinal to x
    Jxkx = hcat([m/√N*sum(vx .* cos.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)
    JxkxIm = hcat([m/√N*sum(vx .* sin.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)

    #Transversal to y
    Jyky = hcat([m/√N*sum(vy .* cos.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)
    JykyIm = hcat([m/√N*sum(vy .* sin.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)

    #Transversal to z
    Jzkz = hcat([m/√N*sum(vz .* cos.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)
    JzkzIm = hcat([m/√N*sum(vz .* sin.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)

    Jxkx, Jyky, Jzkz, JxkxIm, JykyIm, JzkzIm
end

function transversalDensityCurrent(data, N; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates Jₓ(k,t) = m/N * Σ vₓⁱ(t) exp(-𝚤 k ⋅ rⁱ(t))
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then Jₓ(0,t) = m/2
    =#
    x, y, z, vx, vy, vz = data
    #Transversal to z
    Jxkz = hcat([m/√N*sum(vx .* cos.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)
    JxkzIm = hcat([m/√N*sum(vx .* sin.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)

    Jykz = hcat([m/√N*sum(vy .* cos.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)
    JykzIm = hcat([m/√N*sum(vy .* sin.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)
    #Transversal to x
    Jzkx = hcat([m/√N*sum(vz .* cos.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)
    JzkxIm = hcat([m/√N*sum(vz .* sin.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)

    Jykx = hcat([m/√N*sum(vy .* cos.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)
    JykxIm = hcat([m/√N*sum(vy .* sin.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)
    #Transversal to y
    Jxky = hcat([m/√N*sum(vx .* cos.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)
    JxkyIm = hcat([m/√N*sum(vx .* sin.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)

    Jzky = hcat([m/√N*sum(vz .* cos.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)
    JzkyIm = hcat([m/√N*sum(vz .* sin.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)

    Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx, JxkzIm, JykzIm, JxkyIm, JzkyIm, JykxIm, JzkxIm
end

function longitudinalDensityCorrelation(data, N, tmax; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates ⟨ JxJx ⟩(k,τ) =  Σₜ Jₓ(k,t+τ)Jₓ(k,t)/Z where Z = tmax-τ
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then ⟨ JxJx ⟩(k,τ) = m²/4
    =#

    Jxkx, Jyky, Jzkz, JxkxIm, JykyIm, JzkzIm = longitudinalDensityCurrent(data, N; kmin = kmin, kmax = kmax, L = L, m=m)

    JxJxkx = hcat([sum([Jxkx[τ+t,:] .* Jxkx[t,:] .+ JxkxIm[τ+t,:] .* JxkxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    # JxJxky = hcat([sum([Jxky[τ+t,:] .* Jxky[t,:] .+ JxkyIm[τ+t,:] .* JxkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkxIm = hcat([sum([Jxkx[τ+t,:] .* JxkxIm[t,:] .+ JxkxIm[τ+t,:] .* Jxkx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    # JxJxkyIm = hcat([sum([Jxky[τ+t,:] .* JxkyIm[t,:] .+ JxkyIm[τ+t,:] .* Jxky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JyJyky = hcat([sum([Jyky[τ+t,:] .* Jyky[t,:] .+ JykyIm[τ+t,:] .* JykyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    # JyJykx = hcat([sum([Jykx[τ+t,:] .* Jykx[t,:] .+ JykxIm[τ+t,:] .* JykxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JyJykyIm = hcat([sum([Jyky[τ+t,:] .* JykyIm[t,:] .+ JykyIm[τ+t,:] .* Jyky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    # JyJykxIm = hcat([sum([Jykx[τ+t,:] .* JykxIm[t,:] .+ JykxIm[τ+t,:] .* Jykx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JzJzkz = hcat([sum([Jzkz[τ+t,:] .* Jzkz[t,:] .+ JzkzIm[τ+t,:] .* JzkzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    # JzJzky = hcat([sum([Jzky[τ+t,:] .* Jzky[t,:] .+ JzkyIm[τ+t,:] .* JzkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JzJzkzIm = hcat([sum([Jzkz[τ+t,:] .* JzkzIm[t,:] .+ JzkzIm[τ+t,:] .* Jzkz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    # JzJzkyIm = hcat([sum([Jzky[τ+t,:] .* JzkyIm[t,:] .+ JzkyIm[τ+t,:] .* Jzky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkx, JyJyky, JzJzkz, JxJxkxIm, JyJykyIm, JzJzkzIm
end

function longitudinalDensityCorrelation(data, N, tmax, β::Int; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates ⟨ JxJx ⟩(k,τ) =  Σₜ Jₓ(k,t+τ)Jₓ(k,t)/Z where Z = tmax-τ
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then ⟨ JxJx ⟩(k,τ) = m²/4
    =#

    Jxkx, Jyky, Jzkz, JxkxIm, JykyIm, JzkzIm = longitudinalDensityCurrent(data, N; kmin = kmin, kmax = kmax, L = L, m=m)

	JxJxkx = hcat([diag((Jxkx' * sparse_with_diag(i,tmax) * Jxkx .+ JxkxIm' * sparse_with_diag(i,tmax) * JxkxIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkz = hcat([sum([Jxkz[τ+t,:] .* Jxkz[t,:] .+ JxkzIm[τ+t,:] .* JxkzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	# JxJxky = hcat([diag((Jxky' * sparse_with_diag(i,tmax) * Jxky .+ JxkyIm' * sparse_with_diag(i,tmax) * JxkyIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxky = hcat([sum([Jxky[τ+t,:] .* Jxky[t,:] .+ JxkyIm[τ+t,:] .* JxkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JxJxkxIm = hcat([diag((Jxkx' * sparse_with_diag(i,tmax) * JxkxIm .+ JxkxIm' * sparse_with_diag(i,tmax) * Jxkx)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkzIm = hcat([sum([Jxkz[τ+t,:] .* JxkzIm[t,:] .+ JxkzIm[τ+t,:] .* Jxkz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	# JxJxkyIm = hcat([diag((Jxky' * sparse_with_diag(i,tmax) * JxkyIm .+ JxkyIm' * sparse_with_diag(i,tmax) * Jxky)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkyIm = hcat([sum([Jxky[τ+t,:] .* JxkyIm[t,:] .+ JxkyIm[τ+t,:] .* Jxky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JyJyky = hcat([diag((Jyky' * sparse_with_diag(i,tmax) * Jyky .+ JykyIm' * sparse_with_diag(i,tmax) * JykyIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykz = hcat([sum([Jykz[τ+t,:] .* Jykz[t,:] .+ JykzIm[τ+t,:] .* JykzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	# JyJykx = hcat([diag((Jykx' * sparse_with_diag(i,tmax) * Jykx .+ JykxIm' * sparse_with_diag(i,tmax) * JykxIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykx = hcat([sum([Jykx[τ+t,:] .* Jykx[t,:] .+ JykxIm[τ+t,:] .* JykxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JyJykyIm = hcat([diag((Jyky' * sparse_with_diag(i,tmax) * JykyIm .+ JykyIm' * sparse_with_diag(i,tmax) * Jyky)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykzIm = hcat([sum([Jykz[τ+t,:] .* JykzIm[t,:] .+ JykzIm[τ+t,:] .* Jykz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	# JyJykxIm = hcat([diag((Jykx' * sparse_with_diag(i,tmax) * JykxIm .+ JykxIm' * sparse_with_diag(i,tmax) * Jykx)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykxIm = hcat([sum([Jykx[τ+t,:] .* JykxIm[t,:] .+ JykxIm[τ+t,:] .* Jykx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JzJzkz = hcat([diag((Jzkz' * sparse_with_diag(i,tmax) * Jzkz .+ JzkzIm' * sparse_with_diag(i,tmax) * JzkzIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkx = hcat([sum([Jzkx[τ+t,:] .* Jzkx[t,:] .+ JzkxIm[τ+t,:] .* JzkxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	# JzJzky = hcat([diag((Jzky' * sparse_with_diag(i,tmax) * Jzky .+ JzkyIm' * sparse_with_diag(i,tmax) * JzkyIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzky = hcat([sum([Jzky[τ+t,:] .* Jzky[t,:] .+ JzkyIm[τ+t,:] .* JzkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JzJzkzIm = hcat([diag((Jzkz' * sparse_with_diag(i,tmax) * JzkzIm .+ JzkzIm' * sparse_with_diag(i,tmax) * Jzkz)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkxIm = hcat([sum([Jzkx[τ+t,:] .* JzkxIm[t,:] .+ JzkxIm[τ+t,:] .* Jzkx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	# JzJzkyIm = hcat([diag((Jzky' * sparse_with_diag(i,tmax) * JzkyIm .+ JzkyIm' * sparse_with_diag(i,tmax) * Jzky)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkyIm = hcat([sum([Jzky[τ+t,:] .* JzkyIm[t,:] .+ JzkyIm[τ+t,:] .* Jzky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkx, JyJyky, JzJzkz, JxJxkxIm, JyJykyIm, JzJzkzIm
end

function transversalDensityCorrelation(data, N, tmax; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates ⟨ JxJx ⟩(k,τ) =  Σₜ Jₓ(k,t+τ)Jₓ(k,t)/Z where Z = tmax-τ
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then ⟨ JxJx ⟩(k,τ) = m²/4
    =#

    Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx, JxkzIm, JykzIm, JxkyIm, JzkyIm,
        JykxIm, JzkxIm = transversalDensityCurrent(data, N; kmin = kmin, kmax = kmax, L = L, m=m)

    JxJxkz = hcat([sum([Jxkz[τ+t,:] .* Jxkz[t,:] .+ JxkzIm[τ+t,:] .* JxkzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JxJxky = hcat([sum([Jxky[τ+t,:] .* Jxky[t,:] .+ JxkyIm[τ+t,:] .* JxkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkzIm = hcat([sum([Jxkz[τ+t,:] .* JxkzIm[t,:] .+ JxkzIm[τ+t,:] .* Jxkz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JxJxkyIm = hcat([sum([Jxky[τ+t,:] .* JxkyIm[t,:] .+ JxkyIm[τ+t,:] .* Jxky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JyJykz = hcat([sum([Jykz[τ+t,:] .* Jykz[t,:] .+ JykzIm[τ+t,:] .* JykzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JyJykx = hcat([sum([Jykx[τ+t,:] .* Jykx[t,:] .+ JykxIm[τ+t,:] .* JykxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JyJykzIm = hcat([sum([Jykz[τ+t,:] .* JykzIm[t,:] .+ JykzIm[τ+t,:] .* Jykz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JyJykxIm = hcat([sum([Jykx[τ+t,:] .* JykxIm[t,:] .+ JykxIm[τ+t,:] .* Jykx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JzJzkx = hcat([sum([Jzkx[τ+t,:] .* Jzkx[t,:] .+ JzkxIm[τ+t,:] .* JzkxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JzJzky = hcat([sum([Jzky[τ+t,:] .* Jzky[t,:] .+ JzkyIm[τ+t,:] .* JzkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JzJzkxIm = hcat([sum([Jzkx[τ+t,:] .* JzkxIm[t,:] .+ JzkxIm[τ+t,:] .* Jzkx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JzJzkyIm = hcat([sum([Jzky[τ+t,:] .* JzkyIm[t,:] .+ JzkyIm[τ+t,:] .* Jzky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkz, JxJxky, JyJykz, JyJykx, JzJzkx, JzJzky, JxJxkzIm, JxJxkyIm, JyJykzIm, JyJykxIm, JzJzkxIm, JzJzkyIm
end

function transversalDensityCorrelation(data, N, tmax, β::Int; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates ⟨ JxJx ⟩(k,τ) =  Σₜ Jₓ(k,t+τ)Jₓ(k,t)/Z where Z = tmax-τ
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then ⟨ JxJx ⟩(k,τ) = m²/4
    =#

    Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx, JxkzIm, JykzIm, JxkyIm, JzkyIm,
        JykxIm, JzkxIm = transversalDensityCurrent(data, N; kmin = kmin, kmax = kmax, L = L, m=m)

	JxJxkz = hcat([diag((Jxkz' * sparse_with_diag(i,tmax) * Jxkz .+ JxkzIm' * sparse_with_diag(i,tmax) * JxkzIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkz = hcat([sum([Jxkz[τ+t,:] .* Jxkz[t,:] .+ JxkzIm[τ+t,:] .* JxkzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JxJxky = hcat([diag((Jxky' * sparse_with_diag(i,tmax) * Jxky .+ JxkyIm' * sparse_with_diag(i,tmax) * JxkyIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxky = hcat([sum([Jxky[τ+t,:] .* Jxky[t,:] .+ JxkyIm[τ+t,:] .* JxkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JxJxkzIm = hcat([diag((Jxkz' * sparse_with_diag(i,tmax) * JxkzIm .+ JxkzIm' * sparse_with_diag(i,tmax) * Jxkz)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkzIm = hcat([sum([Jxkz[τ+t,:] .* JxkzIm[t,:] .+ JxkzIm[τ+t,:] .* Jxkz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JxJxkyIm = hcat([diag((Jxky' * sparse_with_diag(i,tmax) * JxkyIm .+ JxkyIm' * sparse_with_diag(i,tmax) * Jxky)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkyIm = hcat([sum([Jxky[τ+t,:] .* JxkyIm[t,:] .+ JxkyIm[τ+t,:] .* Jxky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JyJykz = hcat([diag((Jykz' * sparse_with_diag(i,tmax) * Jykz .+ JykzIm' * sparse_with_diag(i,tmax) * JykzIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykz = hcat([sum([Jykz[τ+t,:] .* Jykz[t,:] .+ JykzIm[τ+t,:] .* JykzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JyJykx = hcat([diag((Jykx' * sparse_with_diag(i,tmax) * Jykx .+ JykxIm' * sparse_with_diag(i,tmax) * JykxIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykx = hcat([sum([Jykx[τ+t,:] .* Jykx[t,:] .+ JykxIm[τ+t,:] .* JykxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JyJykzIm = hcat([diag((Jykz' * sparse_with_diag(i,tmax) * JykzIm .+ JykzIm' * sparse_with_diag(i,tmax) * Jykz)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykzIm = hcat([sum([Jykz[τ+t,:] .* JykzIm[t,:] .+ JykzIm[τ+t,:] .* Jykz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JyJykxIm = hcat([diag((Jykx' * sparse_with_diag(i,tmax) * JykxIm .+ JykxIm' * sparse_with_diag(i,tmax) * Jykx)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykxIm = hcat([sum([Jykx[τ+t,:] .* JykxIm[t,:] .+ JykxIm[τ+t,:] .* Jykx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JzJzkx = hcat([diag((Jzkx' * sparse_with_diag(i,tmax) * Jzkx .+ JzkxIm' * sparse_with_diag(i,tmax) * JzkxIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkx = hcat([sum([Jzkx[τ+t,:] .* Jzkx[t,:] .+ JzkxIm[τ+t,:] .* JzkxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JzJzky = hcat([diag((Jzky' * sparse_with_diag(i,tmax) * Jzky .+ JzkyIm' * sparse_with_diag(i,tmax) * JzkyIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzky = hcat([sum([Jzky[τ+t,:] .* Jzky[t,:] .+ JzkyIm[τ+t,:] .* JzkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JzJzkxIm = hcat([diag((Jzkx' * sparse_with_diag(i,tmax) * JzkxIm .+ JzkxIm' * sparse_with_diag(i,tmax) * Jzkx)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkxIm = hcat([sum([Jzkx[τ+t,:] .* JzkxIm[t,:] .+ JzkxIm[τ+t,:] .* Jzkx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JzJzkyIm = hcat([diag((Jzky' * sparse_with_diag(i,tmax) * JzkyIm .+ JzkyIm' * sparse_with_diag(i,tmax) * Jzky)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkyIm = hcat([sum([Jzky[τ+t,:] .* JzkyIm[t,:] .+ JzkyIm[τ+t,:] .* Jzky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkz, JxJxky, JyJykz, JyJykx, JzJzkx, JzJzky, JxJxkzIm, JxJxkyIm, JyJykzIm, JyJykxIm, JzJzkxIm, JzJzkyIm
end

function sparse_with_diag(idx, tmax)
    A = blockdiag(sparse([i for i in 1+idx:tmax],[i for i in 1:tmax-idx],[1 for i in 1:tmax-idx]), sparse(0I, 0, idx))
    A
end

function getTransversalCorrelationFunction(data, N, tmax; kmin = 1, kmax = 10, L = 1, m=1)
	use_sparse = 1
	if use_sparse == 0
    	CT = transversalDensityCorrelation(data, N, tmax; kmin = kmin, kmax = kmax, L = L, m=m)
		return mean(CT[1:6]), mean(CT[7:12])
	elseif use_sparse == 1
		CT = transversalDensityCorrelation(data, N, tmax,1; kmin = kmin, kmax = kmax, L = L, m=m)
		return  mean(CT[1:6]), mean(CT[7:12])
	end
end

function getLongitudinalCorrelationFunction(data, N, tmax; kmin = 1, kmax = 10, L = 1, m=1)
	use_sparse = 1
	if use_sparse == 0
    	CL = longitudinalDensityCorrelation(data, N, tmax; kmin = kmin, kmax = kmax, L = L, m=m)
		return mean(CL[1:3]), mean(CL[4:6])
	elseif use_sparse == 1
		CL = longitudinalDensityCorrelation(data, N, tmax,1; kmin = kmin, kmax = kmax, L = L, m=m)
		return  mean(CL[1:3]), mean(CL[4:6])
	end
end

max_stars_repo_name: jquetzalcoatl/DynamicStructureFactor

max_issues_repo_name: jquetzalcoatl/DynamicStructureFactor

max_forks_repo_name: jquetzalcoatl/DynamicStructureFactor


text: include("krylov.jl")


text:
using Distributed, DistributedArrays
using LinearAlgebra

# [CO] This can't be called here, otherwise can't call ACEfit from another
# package that doesn't have Distributed, DistributedArrays in its dependencies
# note sure yet how to resolve this
# @everywhere using Distributed, DistributedArrays

using ACEbase: ACEBasis

# DO WE EVEN NEED THIS?
#
# struct LLSQ
#    basis
#    configs
#    solver
#    meta::Dict{String, Any}
# end

# LLSQ(; basis = nothing, configs = nothing,
#        solver = QR(), meta = Dict{String, Any}()) =
#                         LLSQ(basis, configs, solver, meta)

# todo - move this to ACEbase I think
# wcw: fix these hacks
function get_basis(m::Any)
    return m
end
function set_params!(m::Any, x::Any)
    return m
end

# TODO: Maybe we can provide more fine-tuned execution models
#    struct Distributed parameters end  struct Serial ... end and so forth
#    rather than those symbols :serial, :mt, :dist

function llsq!(model, data::AbstractVector, par = :serial; solver = QR())
   basis = get_basis(model) # should return an ACEBasis
   θ, errors = llsq(basis, data, par; solver = solver)
   set_params!(model, θ)
   return θ, errors
end

function llsq(basis, data::AbstractVector, par = :serial; solver = QR())
   if par == :serial
      _iterate = siterate
      A, y, w = asm_llsq(basis, data, _iterate)
   elseif par == :mt
      _iterate = titerate
      A, y = asm_llsq(basis, data, _iterate)
   elseif par == :dist
      _iterate = siterate
      A, y, w = asm_llsq_dist(basis, data, _iterate)
   else
      error("unknown assembly type")
   end

   coef = solve_llsq(solver, Diagonal(w)*A, w.*y)

   config_errors = error_llsq(data, A*coef, y)

   return coef, config_errors

end

function get_lsq_indices(data)
   # count the number of observations and assign indices in the lsq matrix
   # we do this always in serial since it should take essentially no time
   # (but what if the data lives distributed???)
   Nobs = 0
   firstidx = zeros(Int, length(data))
   function count_Nobs(i, dat)
      firstidx[i] = Nobs + 1
      for o in ACEfit.observations(dat)
         Nobs += length(ACEfit.vec_obs(o))
      end
   end
   ACEfit.siterate(count_Nobs, data)
   return firstidx, Nobs
end

function asm_llsq(basis, data, _iterate)
   _, Nobs = get_lsq_indices(data)

   # allocate - maybe we need to check somewhere what the
   #            eltypes of A, y should be? real, complex, precision?
   A = zeros(Nobs, length(basis))
   Y = zeros(Nobs)
   W = zeros(Nobs)

   # inner assembly (this knows about A and Y)
   idx = 1
   # TODO: i0 is not used anymore, should revisit
   function asm_lsq_inner(i0, dat)
      for o in observations(dat)
         # TODO: this isn't type stable; for very cheap models, this inner
         #       loop could be a bottleneck, can it be fixed?
         oB = basis_obs(typeof(o), basis, dat.config)
         y = vec_obs(o)
         w = get_weight(o)
         # TODO: make this an input parameter eventually
         if hasproperty(o, :E) || hasproperty(o, :V)
            w = w ./ sqrt(length(dat.config))
         end
         inds = idx:idx+length(y)-1
         Y[inds] .= y[:]
         W[inds] .= w*ones(length(y))
         for ib = 1:length(basis)
            ovec = vec_obs(oB[ib])
            A[inds, ib] .= ovec[:]
         end
         idx += length(y)
      end
      return nothing
   end

   _iterate(asm_lsq_inner, data)

   return A, Y, W
end

function asm_llsq_dist(basis, data, _iterate)

   # for now, all data is sent to all workers, should revisit
   _, Nobs = get_lsq_indices(data)

   A = dzeros(Nobs, length(basis))
   Y = dzeros(Nobs)
   W = dzeros(Nobs)

   idx = 0
   # TODO: i0 is not used anymore, should revisit
   function asm_lsq_inner(i0, dat)
      for o in observations(dat)
         oB = basis_obs(typeof(o), basis, dat.config)
         y = vec_obs(o)
         w = get_weight(o)
         # TODO: make this an input parameter eventually
         if hasproperty(o, :E) || hasproperty(o, :V)
            w = w ./ sqrt(length(dat.config))
         end
         localrows = localindices(Y)[1]
         # fill rows of Y and A
         for i in 1:length(y)
            if (idx+i)<localrows[1] || (idx+i)>localrows[end]
               continue
            end
            localpart(Y)[idx+i-localrows[1]+1] = y[i]
            localpart(W)[idx+i-localrows[1]+1] = w
            for ib = 1:length(basis)
               ovec = vec_obs(oB[ib])
               localpart(A)[idx+i-localrows[1]+1, ib] = ovec[i]
            end
         end
         idx += length(y)
      end
      return nothing
   end

   @sync [@spawnat w _iterate(asm_lsq_inner, data) for w in workers()]

   # send data to main process for use with any solver
   #A = convert(Array, A)
   #Y = convert(Vector, Y)

   return A, Y, W
end

function error_llsq(data, approx, exact)

   errors = approx - exact
   config_types = String[]
   config_counts = Dict("set"=>Dict("E"=>0, "F"=>0, "V"=>0))
   config_errors = Dict("set"=>Dict("E"=>0.0, "F"=>0.0, "V"=>0.0))
   config_norms = Dict("set"=>Dict("E"=>0.0, "F"=>0.0, "V"=>0.0))
   for dat in data
       if !(dat.configtype in config_types)
          push!(config_types, dat.configtype)
          merge!(config_counts, Dict(dat.configtype=>Dict("E"=>0,   "F"=>0,   "V"=>0)))
          merge!(config_errors, Dict(dat.configtype=>Dict("E"=>0.0, "F"=>0.0, "V"=>0.0)))
          merge!(config_norms, Dict(dat.configtype=>Dict("E"=>0.0, "F"=>0.0, "V"=>0.0)))
       end
   end

   i = 1
   for dat in data
      for o in observations(dat)
         obs_len = length(vec_obs(o))
         obs_errors = errors[i:i+obs_len-1]
         obs_values = exact[i:i+obs_len-1]
         # TODO: we store the ref energy because it is used for the relrmse
         #       calculation ... but does it make sense to use the total energy?
         if hasproperty(o, :E)
            obs_values = obs_values .+ o.E_ref
         end
         if hasproperty(o, :E) || hasproperty(o, :V)
            obs_errors = obs_errors ./ length(dat.config)
            obs_values = obs_values ./ length(dat.config)
         end
         obs_error = sum(obs_errors.^2)
         obs_norm = sum(obs_values.^2)
         if hasproperty(o, :E)
            config_counts["set"]["E"] += obs_len
            config_errors["set"]["E"] += obs_error
            config_norms["set"]["E"] += obs_norm
            config_counts[dat.configtype]["E"] += obs_len
            config_errors[dat.configtype]["E"] += obs_error
            config_norms[dat.configtype]["E"] += obs_norm
         elseif hasproperty(o, :F)
            config_counts["set"]["F"] += obs_len
            config_errors["set"]["F"] += obs_error
            config_norms["set"]["F"] += obs_norm
            config_counts[dat.configtype]["F"] += obs_len
            config_errors[dat.configtype]["F"] += obs_error
            config_norms[dat.configtype]["F"] += obs_norm
         elseif hasproperty(o, :V)
            config_counts["set"]["V"] += obs_len
            config_errors["set"]["V"] += obs_error
            config_norms["set"]["V"] += obs_norm
            config_counts[dat.configtype]["V"] += obs_len
            config_errors[dat.configtype]["V"] += obs_error
            config_norms[dat.configtype]["V"] += obs_norm
         else
            println("something is wrong")
         end
         i += obs_len
      end
   end

   for i in keys(config_errors)
      for j in keys(config_errors[i])
         config_errors[i][j] = sqrt(config_errors[i][j] / config_counts[i][j])
         config_norms[i][j] = sqrt(config_norms[i][j] / config_counts[i][j])
         config_norms[i][j] = config_errors[i][j] / config_norms[i][j]
      end
   end

   return Dict("rmse"=>config_errors, "relrmse"=>config_norms)

end


text:
export brainscan_poincaresos

"""
    brainscan_poincaresos(A::Dataset, j::Int; kwargs...)
    brainscan_poincaresos(As::Vector{Dataset}, j::Int; kwargs...)
Launch an interactive application for scanning a Poincare surface of section of `A`
like a "brain scan", where the plane that defines the section can be arbitrarily
moved around via a slider. Return `figure, ax3D, ax2D`.

The input dataset must be 3 dimensional, and here the crossing plane is always
chosen to be when the `j`-th variable of the dataset crosses a predefined value.
The slider automatically gets all possible values the `j`-th variable can obtain.

If given multiple datasets, the keyword `colors` attributes a color to each one, e.g.
`colors = [JULIADYNAMICS_COLORS[mod1(i, 6)] for i in 1:length(As)]`.

The keywords `linekw, scatterkw` are named tuples or dictionaries that are propagated to
as keyword arguments to the line
and scatter plot respectively, while the keyword `direction = -1` is propagated
to the function `DyamicalSystems.poincaresos`.
"""
function brainscan_poincaresos(A::DynamicalSystems.AbstractDataset, j::Int; kwargs...)
    brainscan_poincaresos([A], j; kwargs...)
end

function brainscan_poincaresos(
    As::Vector{<:DynamicalSystems.AbstractDataset}, j::Int;
    linekw = (), scatterkw = (), direction = -1,
    colors = [CYCLIC_COLORS[i] for i in 1:length(As)]
)

for A in As; @assert size(A, 2) == 3; end
@assert j ∈ 1:3
mi, ma = total_minmaxima(As)
ms = 25maximum(abs(ma[i] - mi[i]) for i in 1:3)

otheridxs = DynamicalSystems.SVector(setdiff(1:3, j)...)

figure = Figure(resolution = (2000, 800))
display(figure)
ax = figure[1, 1] = Axis3(figure)
axp = figure[1, 2] = Axis(figure)
sll = labelslider!(
    figure, "$(('x':'z')[j]) =", range(mi[j], ma[j]; length = 100);
    sliderkw = Dict(:startvalue => (ma[j]+mi[j])/2)
)
figure[2, :] = sll.layout
y = sll.slider.value

for i in 1:length(As)
    A = As[i]
    # plot 3D trajectory
    lines!(ax, A.data; color = colors[i], transparency = false, linekw...)

    # Poincare sos
    psos = lift(y) do y
        DynamicalSystems.poincaresos(A, (j, y); direction, warning = false)
    end
    psos2d = lift(p -> p[:, otheridxs].data, psos)
    psos3d = lift(p -> p.data, psos)

    Makie.scatter!(axp, psos2d; color = colors[i], scatterkw...)
    Makie.scatter!(ax, psos3d; color = colors[i], markersize = ms, scatterkw...)
end

xlims!(axp, mi[otheridxs[1]], ma[otheridxs[1]])
ylims!(axp, mi[otheridxs[2]], ma[otheridxs[2]])

# plot transparent plane
ss = [mi...]
ws = [(ma - mi)...]
ws[j] = 0

p = lift(y) do y
    ss[j] = y
    o = Makie.Point3f0(ss...)
    w = Makie.Point3f0(ws...)
    Makie.FRect3D(o, w)
end

a = Makie.RGBAf0(0,0,0,0)
c = Makie.RGBAf0(0.2, 0.2, 0.25, 1.0)
img = Makie.ImagePattern([c a; a c]);
Makie.mesh!(ax, p; color = img);

return figure, ax, axp
end

function total_minmaxima(As::Vector{<:DynamicalSystems.AbstractDataset})
    mi, ma = DynamicalSystems.minmaxima(As[1])
    mi = Vector(mi); ma = Vector(ma)
    for j in 2:length(As)
        mi2, ma2 = DynamicalSystems.minmaxima(As[j])
        for i in 1:size(As[j], 2)
            mi[i] = min(mi2[i], mi[i])
            ma[i] = max(ma2[i], ma[i])
        end
    end
    return mi, ma
end


text:
@testset "Visualizations" begin
    ω = [0.1, 0.1]
    A = [0.5 0; 0 0.5]
    B = [0.5 0; 0 0.5]
    simulation = simulate_GAS_1_1(Normal, 0.0, ω, A, B, 1)
    # LBFGS
    gas = Model(1, 1, Normal, 0.0)
    f = fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f)
    @test length(rec) == 7
    ω = [0.1, 0.1]
    A = [0.5 0; 0 0.5]
    B = [0.5 0; 0 0.5]
    simulation = simulate_GAS_1_1(Normal, 0.0, ω, A, B, 1)
    gas = ScoreDrivenModels.Model(1, 1, Normal, 0.0)
    bac = cross_validation(gas, simulation, 10, 4985)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), bac, "name")
    @test length(rec) == 2

    simulation = simulate_GAS_1_1(Normal, 0.0, ω, A, B, 1)
    gas = ScoreDrivenModels.Model(1, 1, Normal, 0.0)
    fit!(gas, simulation)
    forec = forecast(simulation, gas, 50)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), simulation, forec)
    @test length(rec) == 4
end


text:
module LinearVortexSource

## Package imports
#==========================================================================================#

using LinearAlgebra
using StaticArrays
using Base.Iterators

import ..MathTools: rotation, inverse_rotation, midpair_map

import ..NonDimensional: pressure_coefficient

import ..PanelGeometry: AbstractPanel2D, Panel2D, WakePanel2D, collocation_point, p1, p2, transform_panel, affine_2D, panel_length, panel_angle, panel_tangent, panel_normal, distance, wake_panel, wake_panels, panel_points, panel_vector, panel_velocity

## Singularities
#==========================================================================================#

include("singularities.jl")

# export constant_source_velocity, linear_source_velocity_a, linear_source_velocity_b, linear_vortex_velocity_a, linear_vortex_velocity_b

## Matrix setups
#==========================================================================================#

include("matrix.jl")

# export total_velocity, source_velocity, vortex_velocity, vortex_influence_matrix, source_influence_matrix, neumann_boundary_condition, kutta, two_point_matrix, source_matrix, vortex_matrix, constant_source_matrix, constant_source_boundary_condition

end


text:
import GraphIdx.LinA: IncMat
import GraphIdx: enumerate_edges, WeightedGraph
import GraphIdx.Tree: WeightedTree
import GraphIdx.Grid: GridGraph

"""
    gap_vec!(γ, x, α, g [, c = 1.0])

Compute the duality gap: for each edge `e` we set `γ[e]` to the
difference of the edge cost ``\\lambda_{ij} |x_i - x_j|`` and the
linear approximation via the dual ``α_{ij}(x_i - x_j)``.

The graph `g` needs to implement a method `enumerate_edges(::Function, g)`.

If provided, the result is multiplied by `c`.
"""
function gap_vec!(
    γ::Vector{F},
    x::Array{F},
    α::Vector{F},
    g::Union{WeightedGraph, WeightedTree},
    c::F = F(1.0),
) where {F<:Real,Graph}
    enumerate_edges(g) do ei::Int, u::Int, v::Int, lam::Float64
        let diff = x[u] - x[v]
            γ[ei] = c * (lam*abs(diff) + α[ei] * diff)
        end
    end
end

function gap_vec!(
    γ::Vector{F},
    x::Array{F},
    α::Vector{F},
    g::GridGraph,
    c::F = F(1.0),
) where {F<:Real,Graph}
    enumerate_edges(g, weighted=true) do ei::Int, u::Int, v::Int, lam::Float64
        let diff = x[u] - x[v]
            γ[ei] = c * (lam*abs(diff) + α[ei] * diff)
        end
    end
end

"""
    gap_vec!(γ, dif, x, y, D, Dt, α, c)

Compute the gap vector (modifies `dif` and `x`) and stores it in `γ`.

!!! warning
    This is the old method using Julia's `SparseArrays.SparseMatrixCSC`.

"""
function gap_vec!(
    γ::Vector{Float64},
    dif::Vector{Float64},
    x::Vector{Float64},
    y::Vector{Float64},
    D::IncMat,
    Dt::IncMat,
    α::Vector{Float64},
    c::Float64 = Float64(1.0),
)
    mul!(x, Dt, α)
    x .-= y
    mul!(dif, D, x)
    γ .= c .* ((α .* dif) .+ abs.(dif))
    return γ
end

"""
    wolfe_gap_step(x0, α0, x1, α1)

Optimal gap step width according to Wolfe's condition, i.e.
compute the optimal step width ``θ ∈ [0, 1]`` to walk
from `(x0, α0)` to `(x1, α1)` such that the `sum(gap_vec)` will be minimized.
"""
function wolfe_gap_step(
    x0::Array{F},
    α0::Vector{F},
    x1::Array{F},
    α1::Vector{F},
    g::Graph
)::F where {F<:Real,Graph}
    θ::F = 0
    local m = num_edges(g)
    d0 = Vector{F}(undef, m)
    d1 = Vector{F}(undef, m)
    enumerate_edges(g) do ei::Int, u::Int, v::Int, lam::Float64
        d0[ei] = x0[u] - x0[v]
        d1[ei] = x1[u] - x1[v]
    end
    return θ
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
include("utils.jl")
filename = "data/small_earth_snapshots.jld2"
 #filename = "data/earth_snapshots.jld2"
jl_file = jldopen(filename, "r+")

# anomaly
# state_val = jl_file["T"]["1"][:,:,:,1] # ./ jl_file["ρ"]["1"][:,:,:,1]
# state_val = state_val .- mean(state_val, dims = 1)
# velocity
state_val = jl_file["ρw"]["1"][:,:,:,1] ./ jl_file["ρ"]["1"][:,:,:,1]
title_string = "Vertical Velocity"
height_index_start = 1
height_index_end = 31

clims = quantile(state_val[:,:,height_index_start:height_index_end][:], 0.99)
clims = (-clims,clims)

fig = Figure(resolution =(1520, 980))
ax = LScene(fig, scenekw = (camera = cam3d!, show_axis = true))
ax_text = Label(fig, title_string,
        textsize = 30, color = (:black, 0.85))

cmap = :balance # :Blues_9
cmapa = RGBAf0.(to_colormap(cmap), 1);
cmap = vcat(cmapa[1:15], fill(RGBAf0(0,0,0,0), 10), cmapa[25:end])

v1 = volume!(ax, 0..20, 0..10, 0..5, state_val[:,30:180-30,height_index_start:height_index_end],
                 colorrange=clims, algorithm=:absorption, absorption=10.0f0,
                 colormap=cmap)
axis = ax.scene[OldAxis]
axis[:names, :axisnames] = ("longitude [ᵒ]", "latitude [ᵒ]", "height [km]")
tstyle = axis[:names] #  get the nested attributes and work directly with them

tstyle[:textsize] = 15
tstyle[:textcolor] = (:black, :black, :black)
tstyle[:font] = "helvetica"
tstyle[:gap] = 10
axis[:ticks][:textcolor] = :black
axis[:ticks][:textsize] = 10
cbar1 = Colorbar(fig, v1, label = L" $w$ [m/s]",width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize=25, tickalign = 1, height = Relative(3/4)
)

axis[:ticks][:ranges] = ([0.0, 5.0, 10.0, 15.0, 20.0], [0.0, 2.5, 5.0, 7.5, 10.0], [0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
axis[:ticks][:labels] = (["180W", "90W", "0", "90E", "180E"], ["60S", "30S", "0", "30N", "60N"], ["0", "6", "12", "18", "24", "30"])

fig[2:10, 1:10] = ax
fig[3:8, 11] = cbar1
fig[1, 5:6] = ax_text

zoom!(ax.scene, 0.8)
update!(fig.scene)

#=
# use this to record video
# bring up the figure and then run this piece of code
seconds = 10
fps = 30
frames = round(Int, fps * seconds )
record(fig, pwd() * "/zonal_velocity.mp4"; framerate = fps) do io
    for i = 1:frames
        sleep(1/fps)
        recordframe!(io)
    end
end
=#

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Clearances_Main
	Clearances_Main()=begin
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Sealing Strips pairs",
				:Lower=>1
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Height Under Shell Inlet Nozzle",
				:Lower=>1E-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Height Under Shell Outlet Nozzle",
				:Lower=>1E-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Bundle-to-Shell Clearance",
				:Symbol=>"L_{cf}",
				:Lower=>1E-8
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Baffle-to-Shell Clearance",
				:Symbol=>"L_{cd}",
				:Lower=>1E-8
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tube-to-Baffle Clearance",
				:Symbol=>"L_{td}",
				:Lower=>1E-8
			)),
			[:SealStrip,:Hinozzle_Shell,:Honozzle_Shell,:BundleToShell,:BaffleToShell,:TubeToBaffle,],
		)
	end
	SealStrip::DanaInteger
	Hinozzle_Shell::length
	Honozzle_Shell::length
	BundleToShell::length
	BaffleToShell::length
	TubeToBaffle::length
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Clearances_Main
function atributes(in::Clearances_Main,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Main parameters for diametral clearances in a shell and tube heat exchanger."
	drive!(fields,_)
	return fields
end
Clearances_Main(_::Dict{Symbol,Any})=begin
	newModel=Clearances_Main()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
import Base: promote_rule

# SparseMatrixCSC
promote_rule(::Type{SparseMatrixCSC{Tv, Ti}}, ::Type{Matrix{T}}) where {Tv, Ti, T} = Matrix{promote_type(T, Tv)}

# IMatrix
promote_rule(::Type{IMatrix{N, T}}, ::Type{PermMatrix{Tv, Ti, Vv, Vi}}) where {N, T, Tv, Ti, Vv, Vi} = (TT=promote_type(T, Tv); PermMatrix{TT, Ti, Vector{TT}, Vi})
promote_rule(::Type{IMatrix{N, T}}, ::Type{SparseMatrixCSC{Tv, Ti}}) where {N, T, Tv, Ti} = SparseMatrixCSC{promote_type(T, Tv), Ti}
promote_rule(::Type{IMatrix{M, TA}}, ::Type{Matrix{TB}}) where {M, TA, TB} = Array{TB, 2}

# PermMatrix
promote_rule(::Type{PermMatrix{TvA, TiA}}, ::Type{SparseMatrixCSC{TvB, TiB}}) where {TvA, TiA, TvB, TiB} =
    SparseMatrixCSC{promote_type(TvA, TvB), promote_type(TiA, TiB)}
promote_rule(::Type{PermMatrix{Tv, Ti}}, ::Type{Matrix{T}}) where {Tv, Ti, T} =
    Array{promote_type(Tv, T), 2}


text:
# Basic manipulation of points on an elliptic curve

mutable struct EllipticCurve
    a1
    a2
    a3
    a4
    a6
end

mutable struct Point
    x
    y
    E
end

function Base.:(-)(a::Point)
    Point(a.x,-a.E.a1*a.x-a.E.a3-a.y,a.E)
end

function Base.:(+)(a::Point, b::Point)
    if !(a.E == b.E)
        error("error: elliptic curve must be the same")
    end
    if (a == b)
        L = (3*a.x^2+2*a.E.a2*a.x-a.E.a1*a.y+a.E.a4)//(2*a.y+a.E.a1*a.x+a.E.a3)
        x = L^2+L*a.E.a1-a.E.a2-2*a.x
        y = L*x-L*a.x+a.y
        Point(x,-y,a.E)
    else
        L = (b.y-a.y)//(b.x-a.x)
        x = L^2+L*a.E.a1-a.E.a2-a.x-b.x
        y = L*x-L*a.x+a.y
        Point(x,-y,a.E)
    end
end

function Base.:(*)(n::Integer, P::Point)
    Q = nothing
    F = P
    m = n
    while m >= 1
        if (m % 2==1)
            if Q==nothing
                Q = F
            else
                Q += F
            end
        end
        m = floor(m/2)
        F += F
    end
    Q
end

E = EllipticCurve(BigInt(0),BigInt(109),BigInt(0),BigInt(224),BigInt(0))

P = Point(BigInt(-100),BigInt(260),E)

println(P)
println(-P)
println(P+P)
Q = 9*P

sa = (56-Q.x+Q.y)//(56-14*Q.x)
sb = (56-Q.x-Q.y)//(56-14*Q.x)
sc = (-28-6*Q.x)//(28-7*Q.x)
l = lcm(denominator(sa),denominator(sb),denominator(sc))
println(numerator(sa*l))
println(numerator(sb*l))
println(numerator(sc*l))

max_stars_repo_name: octonion/elliptic-curves

max_issues_repo_name: octonion/elliptic-curves

max_forks_repo_name: octonion/elliptic-curves


text:
using DataFrames
using DataConvenience

sample(DataFrame(a=1:100), 10)

sample(DataFrame(a=1:100), 0.1)

sample(DataFrame(a=1:100), 1//7)

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
"""

```
sphereradiusfromshellvolume(volume,step)
```

Computes the radius that corresponds to a spherical shell of
a given volume.

"""
function sphereradiusfromshellvolume(volume, step)

    fourthirdsofpi = 4 * pi / 3
    if 3 * step * volume - pi * step^4 <= 0.0
        return 0.0
    end
    rmin =
        (sqrt(3 * pi) * sqrt(3 * step * volume - pi * step^4) - 3 * pi * step^2) /
        (6 * pi * step)
    return (0.5 * (volume / fourthirdsofpi + 2 * rmin^3))^(1 / 3)

end


text:
# Test `getModel` and `solveModel!`

function test_milp()
    c = hcat([3; 5])
    a = [1 0; 0 2; 3 2]
    b = hcat([4; 12; 18])

    model = ModelMilp(vec_c, a, b)
    solve!(model)

    return model.solution.obj
end

@test test_milp() == 42


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

using Test

function recompile_azureclusterlesshpc(; runtime=false)
    azureclusterlesshpc_pkgid = Base.PkgId(Base.UUID("0090f474-c28d-4577-9e1b-9155814dbffa"), "AzureClusterlessHPC")
    if runtime
        ENV["AZ_BATCH_TASK_WORKING_DIR"] = pwd()
    else
        if haskey(ENV, "AZ_BATCH_TASK_WORKING_DIR")
            delete!(ENV, "AZ_BATCH_TASK_WORKING_DIR")
        end
    end
    Base.compilecache(azureclusterlesshpc_pkgid)
end

# Run tests
if length(ARGS) == 0

    # Base tests
    recompile_azureclusterlesshpc(runtime=false)
    include("TestCore.jl")
    using Main.TestCore
    TestCore.runtests()

    # Runtime tests
    recompile_azureclusterlesshpc(runtime=true)
    include("TestRuntime.jl")
    using Main.TestRuntime
    TestRuntime.runtests()
    recompile_azureclusterlesshpc(runtime=false)
else
    if ARGS[1] == "core"

        recompile_azureclusterlesshpc(runtime=false)
        include("TestCore.jl")
        using Main.TestCore
        TestCore.runtests()

    elseif ARGS[1] == "runtime"

        recompile_azureclusterlesshpc(runtime=true)
        include("TestRuntime.jl")
        using Main.TestRuntime
        TestRuntime.runtests()

    else
        print("Pass \"core\" or \"runtime\".\n")
    end
    recompile_azureclusterlesshpc(runtime=false)
end


text:
##########################
#### HELPER FUNCTIONS ####
##########################

confining_pressure(s::NumericalSetup) = s.pc
confining_pressure(m::Model) = confining_pressure(m.setup)

plastic_yield_stress(plas::Plasticity,p) = plas.σy(p)
plastic_yield_stress(m::Model,p) = plastic_yield_stress(m.cm.plasticity,p)

ϵ2ϵdev(ϵax,E,ν,pc) = 2/3 * (1+ν)/E * (E*ϵax + (1 - 2ν)*pc)
ϵ̇2ϵ̇dev(ϵ̇,ν) = 2/3 * (1+ν) * ϵ̇

is3D(g::Geom) = g isa Geom3D ? true : false
is3D(setup::TriaxialSetup) = is3D(setup.geom)
is3D(m::Model) = is3D(m.setup)
# variables initializers : [s_ax,D] or [ϵ_ax, D] depending on the control variable type (ConstantStress or ConstantStrainRate)
init_variables(::SD) = [0.0,0.0]
init_variables(dam::MD) = [0.0,dam.r.D₀]
init_variables(dam::MD, Dᵢ) = [0.0,max(Dᵢ,dam.r.D₀)]
init_variables(m::Model) = init_variables(m.cm.damage)
init_variables(m::Model, Dᵢ) = init_variables(m.cm.damage, Dᵢ)

Dmax_default = 0.95
init_params(::Model{<:CM{tW,tD,tE,Nothing}}    ; Dmax=Dmax_default) where {tW,tD,tE} = (;Dmax)
init_params(::Model{<:CM{tW,tD,tE,<:P{<:DT}}}  ; Dmax=Dmax_default) where {tW,tD,tE} = (;Dmax)
init_params(::Model{<:CM{tW,tD,tE,<:P{<:MVT}}} ; Dmax=Dmax_default) where {tW,tD,tE} = (;Dmax)

# get shear and normal stress on a defect oriented ψ degree wrt σ1, negative in compression by convention
"shear stress from σ₁, σ₃ and ψ"
flaw_shear_stress(σ₁,σ₃,ψ) = 0.5 * (σ₃ - σ₁) * sind(2*ψ)
"normal stress from σ₁, σ₃ and ψ"
flaw_normal_stress(σ₁,σ₃,ψ) =  0.5 *((σ₃ + σ₁) + (σ₃ - σ₁)*cosd(2*ψ))

####################
#### INVARIANTS ####
####################
# geometric factor linking second invariants of the deviatoric stress or strain tensor and s_ax or ϵ'_ax in triaxial setup
geom_τ_factor(setup::TriaxialSetup) = is3D(setup) ? sqrt(3)/2 : 1.0
geom_τ_factor(model::M) = geom_τ_factor(model.setup)

# stress invariants for a triaxial setup,
"stress invariants τ and σₘ, sₐₓ negative in compression"
function stress_invariants(sₐₓ, setup::TriaxialSetup)
    pc = confining_pressure(setup)
	τ_fact = geom_τ_factor(setup)
	if is3D(setup.geom) # all 3D cases
		τ = τ_fact * abs(sₐₓ)
    	σₘ = 0.5*sₐₓ - pc
	else # 2D plane strain incompressible or inplane compressibility only.
		τ = τ_fact * abs(sₐₓ)
    	σₘ = sₐₓ - pc
	end
	return τ, σₘ
end
stress_invariants(s, model::M{<:CM,<:TS}) = stress_invariants(s, model.setup)

# conversions between axial
sₐₓ2σₐₓ(setup::TriaxialSetup,s) = (Δσ=sₐₓ2Δσ(setup,s) ; Δσ - setup.pc)
sₐₓ2σₐₓ(m::Model,s) = sₐₓ2σₐₓ(m.setup,s)

σₐₓ2sₐₓ(setup::TriaxialSetup{Geom3D},σₐₓ) = (2/3)*(σₐₓ + setup.pc)
σₐₓ2sₐₓ(setup::TriaxialSetup{Geom2D},σₐₓ) = 0.5*(σₐₓ + setup.pc)
σₐₓ2sₐₓ(m::Model,σₐₓ) = σₐₓ2sₐₓ(m.setup,σₐₓ)

sₐₓ2Δσ(geom::Geom,s) = is3D(geom) ? (3/2)*s : 2.0*s
sₐₓ2Δσ(setup::TriaxialSetup,s) = sₐₓ2Δσ(setup.geom,s)
sₐₓ2Δσ(model::M,s) = sₐₓ2Δσ(model.setup,s)

Δσ2sₐₓ(geom::Geom,Δσ) = is3D(geom) ? (2/3)*Δσ : 0.5*Δσ
Δσ2sₐₓ(setup::TriaxialSetup,Δσ) = Δσ2sₐₓ(setup.geom,Δσ)
Δσ2sₐₓ(model::M,Δσ) = Δσ2sₐₓ(model.setup,Δσ)

# ϵII for incompressible triaxial experiment as defined in SiStER
ϵII_func(ϵ, model::M{<:CM{tW,tD,<:IE},<:TS}) where {tW,tD} = geom_τ_factor(model)*abs(ϵ)
ϵ̇II_func(ϵ̇, model::M{<:CM{tW,tD,<:IE},<:TS}) where {tW,tD} = ϵII_func(ϵ̇, model) # same for the time derivative

#############################
#### WEAKENING FUNCTIONS ####
#############################
#cm.weakening.max is the residual fraction of G0 remaining at D=1
# Linear simple damage
weakening_func(D,cm::CM{<:LW,<:SD}) = (1-(1-cm.weakening.γ)*D)
weakening_derivative(D,cm::CM{<:LW,<:SD}) = -(1-cm.weakening.γ*one(D))
# Linear simple micromechanical damage
weakening_func(D,cm::CM{<:LW,<:MD}) = ((cm.weakening.γ-1)*D + 1 - cm.damage.r.D₀*cm.weakening.γ) / (1-cm.damage.r.D₀)
weakening_derivative(D,cm::CM{<:LW,<:MD}) = (cm.weakening.γ-1)/(1-cm.damage.r.D₀)
# asymptotic simple damage
weakening_func(D,cm::CM{<:AW,<:SD}) = 1/(1+D)^cm.weakening.γ #TODO: TESTS
weakening_derivative(D,cm::CM{<:AW,<:SD}) = -cm.weakening.γ/(1+D)^(cm.weakening.γ+1) #TODO: TESTS
# asymptotic micromechanical damage untested
weakening_func(D,cm::CM{<:AW,<:MD}) = ((D<cm.damage.r.D₀)&&(D=cm.damage.r.D₀) ; ((1+cm.damage.r.D₀)/(1+D))^cm.weakening.γ) #TODO: TESTS
weakening_derivative(D,cm::CM{<:AW,<:MD}) = -cm.weakening.γ*(1+cm.damage.r.D₀)^(cm.weakening.γ)/(1+D)^(cm.weakening.γ+1) #TODO: TESTS

weakening_func(D,m::Model) = weakening_func(D,m.cm)
weakening_derivative(D,m::Model) = weakening_derivative(D,m.cm)

function get_energy_based_G(D, τ, σ, s₁₂, model::M{<:CM{<:EBW,<:IKI}})
	r = model.cm.damage.r
	G = model.cm.elasticity.G
	A, B = DSB.compute_AB(r,D)
	A1, B1 = DSB.compute_A1B1(r,A,B)
	#return 1 / ( 1/(4G) * ( 1 + (B1/2) * (B1 + A1*σ/τ - A1*σ*s₁₂^2/τ^3) ) )
	return 1/(1 + B1^2/2 + A1*B1*σ /(2*τ))
end

#############
#### KI  ####
#############
"KI with yield stress from stress invariants (sign convention agnostic)"
function get_KI_from_s(s, D, model::M{<:CM{<:W,<:IKI}})
	r = model.cm.damage.r
    #D₀ = r.D₀ #τ_heal = 60
	τ, σₘ = stress_invariants(s, model)
	D = IfElse.ifelse(D < r.D₀, r.D₀, D)
	KI = DSB.compute_KI(r,σₘ,τ,D)
	return KI
end

"KI with yield stress from principal stresses (sign convention agnostic)"
function get_KI_from_s(s, D, model::M{<:CM{<:W,<:PKI}})
	r = model.cm.damage.r
	pc = confining_pressure(model)
	D = IfElse.ifelse(D < r.D₀, r.D₀, D)
	KI = compute_KI(r, sₐₓ2σₐₓ(model,s), -pc, D)
	return KI
end

#######################
#### DAMAGE GROWTH ####
#######################

"damage growth rate when damage type is SimpleCharlesLaw"
function Ḋ_func(s, D, p, model::M{<:CM{tW,<:SD}}) where {tW}
    d = model.cm.damage
    #τ_heal = 0
	τ, σₘ = stress_invariants(s, model)
	pin = -σₘ
    σy = model.cm.damage.σy(pin)
	Ḋ = IfElse.ifelse(((τ-σy)>0) & ((D < p.Dmax) | (tW<:AW)),
		1/(d.Td*(d.S)^d.m) * (τ-σy)^d.m,
		zero(D))
	return Ḋ #- (D-D₀)/τ_heal
end

"damage growth rate when damage type is MicroMechanicalCharlesLaw"
# function Ḋ_func(s, D, p, model::M{<:CM{tW,<:MD}}, τ_heal=Inf) where {tW}
# 	r = model.cm.damage.r
# 	(D >= p.Dmax) && (return -(p.Dmax-r.D₀)/τ_heal)
# 	(D <= 0) && (return 0.0)

# 	KI = get_KI_from_s(s, D, model)
# 	(KI <= 0) && (return -(D-r.D₀)/τ_heal)

# 	#(KI >= 50000) && @show D KI

#   	dDdl = DSB.compute_dDdl(r,D) # damage derivative wrt crack length
# 	#(KI >= 50000) && @show dDdl

#   	dldtmax = 2000.0
#   	dldt = min(r.l̇₀*(KI/r.K₁c)^(r.n),dldtmax)  #Vr cracks growth rate
#   	#@assert dDdl * dldt >= 0
# 	#(KI >= 50000) && @show dDdl*dldt ; println()
#   	return dDdl * dldt - (D-r.D₀)/τ_heal
# end

"damage growth rate when damage type is MicroMechanicalCharlesLaw"
function Ḋ_func(s, D, p, model::M{<:CM{tW,<:MD}}, τ_heal=Inf) where {tW}
	r = model.cm.damage.r
	D = IfElse.ifelse(D > p.Dmax, p.Dmax, D)
	D = IfElse.ifelse(D < r.D₀, r.D₀, D)

	KI = get_KI_from_s(s, D, model)
	KI = IfElse.ifelse(KI <= 0, zero(D), KI)
	#(KI >= 50000) && @show D KI

  	dDdl = DSB.compute_dDdl(r,D) # damage derivative wrt crack length
	#(KI >= 50000) && @show dDdl

  	dldtmax = 2000.0
  	dldt = min(r.l̇₀*(KI/r.K₁c)^(r.n),dldtmax)  #Vr cracks growth rate
	dldt = IfElse.ifelse(D == p.Dmax, zero(D), dldt)

  	#@assert dDdl * dldt >= 0
	#(KI >= 50000) && @show dDdl*dldt ; println()
  	return dDdl * dldt - (D-r.D₀)/τ_heal
end
#############################
#### VISCOSITY FUNCTIONS ####
#############################

function η_dam_func(model,D,Ḋ)
	G = model.cm.elasticity.G
	f = weakening_func(D,model)
	df = weakening_derivative(D,model)
	return - G * f^2 / (Ḋ * df)
end

# plastic viscosity for constant strain rate setup
function η_plas_func(model::M{tCM,<:TS{tG,<:CSR}}, s, ϵ̇) where {tCM,tG}
	_, σₘ = stress_invariants(s, model)
	pin = -σₘ
    σy = plastic_yield_stress(model,pin)
	τ_fact = geom_τ_factor(model)
	#ϵ̇II = ϵ̇II_func(ϵ̇, model)
	return abs(σy/(2*τ_fact*ϵ̇))
end

# plastic viscosity for constant stress setup
function η_plas_func(model::M{tCM,<:TS{tG,<:CS}}, s, ϵ̇) where {tCM,tG}
	_, σₘ = stress_invariants(s, model)
	pin = -σₘ
	σy = plastic_yield_stress(model,pin)
	τ_fact = geom_τ_factor(model)
	#ϵ̇II = ϵ̇II_func(ϵ̇, model)
	return abs(σy/(2*τ_fact*ϵ̇))
end

# when no plasticity, viscosity is always linked to damage
get_viscosity(m::Model{<:CM{tW,tD,tE,Nothing}}, s, ϵ̇, D, Ḋ, p) where {tW,tD,tE} = η_dam_func(m,D,Ḋ)

# when plastic depends on the switch condition from damaged-elastic to elastic-plastic
function get_viscosity(m::Model{<:CM{tW,tD,tE,<:Plasticity{<:MinViscosityThreshold{Nothing}}}}, s, ϵ̇, D, Ḋ, p) where {tW,tD,tE}
	#D₀ = m.cm.damage.r.D₀
	ηp =  η_plas_func(m,s,ϵ̇)

	#τ, σₘ = stress_invariants(s, m)
	#pin = -σₘ
	#σy = plastic_yield_stress(m,pin)
	ηd = η_dam_func(m,D,Ḋ)

	#isaboveyield_τ = IfElse.ifelse(τ >= σy, true, false)
	#Dmincond = IfElse.ifelse(D >= D₀ + 0.5*(p.Dmax-D₀), true, false)
	η = IfElse.ifelse(ηd <= ηp, ηp, ηd)
	η = IfElse.ifelse(D >= p.Dmax, ηp, ηd)

	return η
end

function get_viscosity(m::Model{<:CM{tW,tD,tE,<:Plasticity{<:DamageThreshold}}}, s, ϵ̇, D, Ḋ, p) where {tW,tD,tE}
	return IfElse.ifelse(D >= m.cm.plasticity.threshold.value, η_plas_func(m,s,ϵ̇), η_dam_func(m,D,Ḋ))
end

#is_plastic(m::Model{<:CM{tW,tD,tE,<:Plasticity{DamageThreshold}}},s,D) = (D >= m.cm.plasticity.threshold.value)
#is_plastic(m::Model{<:CM{tW,tD,tE,<:Plasticity{MinViscosityThreshold}}},s,D) =

##########################
#### UPDATE FUNCTIONS ####
##########################

## INPLACE functions ##
### Constant strain rate :
"inplace derivative update function for constant strain rate test, variables are the deviatoric axial stress and the damage"
function update_derivatives!(du,u,p,t,model::M{tCM,<:TS{tG,<:CSR}}) where {tCM,tG}
    G = model.cm.elasticity.G
    ϵ̇ = model.setup.control.ϵ̇
	s, D = u
	Ḋ = Ḋ_func(s, D, p, model)
	η = get_viscosity(model, s, ϵ̇, D, Ḋ, p)
	fw = weakening_func(D,model)
	du[1] = 2*G*fw * (ϵ̇ - s/(2η))
	du[2] = Ḋ
    du
end

### Constant stress :
"inplace derivative update function for constant strain rate test, variables are the axial strain and the damage"
function update_derivatives!(du,u,p,t,model::M{tCM,<:TS{tG,<:CS}}) where {tCM,tG}
    s = model.setup.control.s
	ϵ, D = u
	ϵ̇_prev = du[1]
	Ḋ = Ḋ_func(s, D, p, model)
	η = get_viscosity(model, s, ϵ̇_prev, D, Ḋ, p)
	du[1] = s/(2η)
	du[2] = Ḋ
    du
end

## OUT OF PLACE functions ##
"derivative update function for constant strain rate test, variables are the deviatoric axial stress and the damage"
function update_derivatives(u,p,t,model::M{tCM,<:TS{tG,<:CSR}}) where {tCM,tG}
    G = model.cm.elasticity.G
    ϵ̇ = model.setup.control.ϵ̇
	s, D = u
	Ḋ = Ḋ_func(s, D, p, model)
	η = get_viscosity(model, s, ϵ̇, D, Ḋ, p)
	fw = weakening_func(D,model)
	ṡ = 2*G*fw * (ϵ̇ - s/(2η))
    return SA[ṡ, Ḋ]
end

"derivative update function for constant strain rate test, variables are the axial strain and the damage"
function update_derivatives(u,p,t,model::M{tCM,<:TS{tG,<:CS}}) where {tCM,tG}
    s = model.setup.control.s
	ϵ, D = u
	ϵ̇_prev = du[1]
	Ḋ = Ḋ_func(s, D, p, model)
	η = get_viscosity(model, s, ϵ̇_prev, D, Ḋ, p)
	ϵ̇ = s/(2η)
    return SA[ϵ̇, Ḋ]
end

max_stars_repo_name: Leooop/SimpleDamage

max_issues_repo_name: Leooop/SimpleDamage

max_forks_repo_name: Leooop/SimpleDamage


text:
using Test
using ScottishTaxBenefitModel
using .EquivalenceScales:
    EQScales,
    EQ_Person,
    get_equivalence_scales,
    EQ_P_Type,
    eq_spouse_of_head,
    eq_other_adult,
    eq_head

using .ExampleHelpers

#
# Replicates cases based on tables 1,2 of:
# Chanfreau, Jenny, and Tania Burchard. 2008.
# ‘Scottish Government - Income and Poverty Statistics - Equivalence Scales Paper’. September 2008. https://www2.gov.scot/Topics/Statistics/Browse/Social-Welfare/IncomePoverty/equivalence-scales-paper.
# https://www2.gov.scot/Topics/Statistics/Browse/Social-Welfare/IncomePoverty/equivalence-scales-paper
#
@testset "Eq Scales" begin
    for (key,hh) in get_all_examples()
        eqs :: EQScales = get_equivalence_scales(
            Float64,
            collect(values(hh.people)))
        println( "hh $key $eqs" )

        if key == single_hh
            @test eqs.oxford == 1
            @test eqs.oecd_bhc == 1
            @test eqs.oecd_ahc == 1
            @test eqs.mcclements_bhc == 1
            @test eqs.mcclements_ahc == 1
            @test eqs.square_root == 1
            @test eqs.per_capita == 1
        elseif key == single_parent_hh
            # 1 adult 2 kids 18 and 8
            @test eqs.oxford ≈ 1+0.7+0.5
            @test eqs.oecd_bhc ≈ 1+0.3+0.5
            @test eqs.oecd_ahc ≈ 1+0.34+0.72
            @test eqs.mcclements_bhc ≈ 1+0.377+0.59
            @test eqs.mcclements_ahc ≈ 1+0.42+0.69
            @test eqs.square_root ≈ sqrt(3)
            @test eqs.per_capita ≈ 3
        elseif key == childless_couple_hh
            @test eqs.oxford ≈ 1+0.7
            @test eqs.oecd_bhc ≈ 1+0.5
            @test eqs.oecd_ahc ≈ 1+0.72
            @test eqs.mcclements_bhc ≈ 1+0.64
            @test eqs.mcclements_ahc ≈ 1+0.82
            @test eqs.square_root == sqrt(2)
            @test eqs.per_capita == 2
        elseif key == cpl_w_2_children_hh
            # chidren 2 5
            @test eqs.oxford ≈ 1+0.7+0.5+0.5
            @test eqs.oecd_bhc ≈ 1+0.5+0.3+0.3
            @test eqs.oecd_ahc ≈ 1+0.72+0.34+0.34
            @test eqs.mcclements_bhc ≈ 1+0.64+0.295+0.344
            @test eqs.mcclements_ahc ≈ 1+0.82+0.33+0.38
            @test eqs.square_root ≈  sqrt(4)
            @test eqs.per_capita ≈ 4
        end
        @test eqs == hh.equivalence_scales
    end
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
struct PkgId
    uuid::Union{UUID,Nothing}
    name::String

    PkgId(u::UUID, name::AbstractString) = new(UInt128(u) == 0 ? nothing : u, name)
    PkgId(::Nothing, name::AbstractString) = new(nothing, name)
end
PkgId(name::AbstractString) = PkgId(nothing, name)

function PkgId(m::Module, name::String = String(nameof(moduleroot(m))))
    uuid = UUID(ccall(:jl_module_uuid, NTuple{2, UInt64}, (Any,), m))
    UInt128(uuid) == 0 ? PkgId(name) : PkgId(uuid, name)
end

==(a::PkgId, b::PkgId) = a.uuid == b.uuid && a.name == b.name

function hash(pkg::PkgId, h::UInt)
    h += 0xc9f248583a0ca36c % UInt
    h = hash(pkg.uuid, h)
    h = hash(pkg.name, h)
    return h
end

show(io::IO, pkg::PkgId) =
    print(io, pkg.name, " [", pkg.uuid === nothing ? "top-level" : pkg.uuid, "]")

function binpack(pkg::PkgId)
    io = IOBuffer()
    write(io, UInt8(0))
    uuid = pkg.uuid
    write(io, uuid === nothing ? UInt128(0) : UInt128(uuid))
    write(io, pkg.name)
    return String(take!(io))
end

function binunpack(s::String)
    io = IOBuffer(s)
    @assert read(io, UInt8) === 0x00
    uuid = read(io, UInt128)
    name = read(io, String)
    return PkgId(UUID(uuid), name)
end


text:
function calc_P_W(reactor, R_0, B_0)
  nu_term = ( 1 + reactor.nu_n ) ^ 2 * ( 1 + reactor.nu_T ) ^ 2
  nu_term /= ( 1 + 2 * ( reactor.nu_n + reactor.nu_T ) )

  K_W = (
    2.824e-3 * nu_term * (reactor.f_D^2) * reactor.epsilon * reactor.kappa ^0.4
  )

  tmp_n_20 = calc_n_20(reactor, R_0, B_0)
  tmp_T_k = calc_T_k(reactor, R_0, B_0)

  K_W * tmp_n_20^2 * tmp_T_k^2 * R_0
end


text:
using PowerModelsReliability
using PowerModels
using InfrastructureModels
using Ipopt
using Mosek
using Juniper
using Cbc
using CPLEX

mosek = MosekSolver()
ipopt = IpoptSolver()
juniper = JuniperSolver(IpoptSolver(print_level=0); mip_solver=CplexSolver())

data_tf = PowerModels.parse_file("./test/data/case5_tf.m")
data_tf["multinetwork"] = false
data_tf["load"] = data_tf["tf_load"]
resultDCunittf = run_unittfopf(data_tf, DCPPowerModel, mosek; multinetwork=false, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => true))

data_tf = PowerModels.parse_file("./test/data/case5_tf.m")
data_tf["multinetwork"] = false
data_tf["load"] = data_tf["tf_load"]
resultACunittf = run_unittfopf(data_tf, ACPPowerModel, juniper; multinetwork=false, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => true))

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
using BenchmarkTools, SymbolicUtils
using SymbolicUtils: isnumber

using Random

SUITE = BenchmarkGroup()

@syms a b c d; Random.seed!(123);

let r = @rule(~x => ~x), rs = RuleSet([r]),
    acr = @rule(~x::isnumber + ~y => ~y)

    overhead = SUITE["overhead"]  = BenchmarkGroup()
    overhead["rule"]  = BenchmarkGroup()

    overhead["rule"]["noop:Int"]  = @benchmarkable $r(1)
    overhead["rule"]["noop:Sym"]  = @benchmarkable $r($a)
    overhead["rule"]["noop:Term"] = @benchmarkable $r($(a+2))

    overhead["acrule"]  = BenchmarkGroup()
    overhead["acrule"]["noop:Int"]  = @benchmarkable $acr(1)
    overhead["acrule"]["noop:Sym"]  = @benchmarkable $acr($a)
    overhead["acrule"]["a+2"] = @benchmarkable $acr($(a+2))
    overhead["acrule"]["a+b"] = @benchmarkable $acr($(a+b))
    overhead["acrule"]["a+2+b"] = @benchmarkable $acr($(a+2+b))

    overhead["ruleset"]  = BenchmarkGroup()
    overhead["ruleset"]["noop:Int"]  = @benchmarkable $rs(1)
    overhead["ruleset"]["noop:Sym"]  = @benchmarkable $rs($a)
    overhead["ruleset"]["noop:Term"] = @benchmarkable $rs($(a+2))

    overhead["simplify"]  = BenchmarkGroup()
    overhead["simplify"]["noop:Int"]  = @benchmarkable simplify(1)
    overhead["simplify"]["noop:Sym"]  = @benchmarkable simplify($a)
    overhead["simplify"]["noop:Term"] = @benchmarkable simplify($(a+2))

    function random_term(len; atoms, funs, fallback_atom=1)
        xs = rand(atoms, len)
        while length(xs) > 1
            xs = map(Iterators.partition(xs, 2)) do xy
                x = xy[1]; y = get(xy, 2, fallback_atom)
                rand(funs)(x, y)
            end
        end
        xs[]
    end
    ex1 = random_term(1000, atoms=[a, b, c, d, a^(-1), b^(-1), 1, 2.0], funs=[+, *])
    ex2 = random_term(1000, atoms=[a, b, c, d, a^(-1), b^(-1), 1, 2.0], funs=[/, *])

    overhead["simplify"]["randterm (+, *):serial"] = @benchmarkable simplify($ex1, threaded=false)
    overhead["simplify"]["randterm (/, *):serial"] = @benchmarkable simplify($ex2, threaded=false)
    overhead["simplify"]["randterm (+, *):thread"] = @benchmarkable simplify($ex1, threaded=true)
    overhead["simplify"]["randterm (/, *):thread"] = @benchmarkable simplify($ex2, threaded=true)
end


text:
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Flux
using Flux: Optimise.update!

function act(model, s, ϵ, nA)
    rand() <= ϵ && return rand(1:nA)
    return Flux.onecold(model(s))
end

function run!(model, opt, env;
        epochs=100, steps=300, maxt=100, batchsize=128,
        trainiters=100, ϵ=0.03, γ=99f-2, τ=1f-2)
    rewards = zeros(Float32, epochs)
    target = deepcopy(model)
    memory = ReplayMemory{Float32, Int64, Float32}(
        length(env.observationspace), 1, 10*steps, batchsize
    )
    nA = length(env.actionspace)
    ps = params(model)

    for i in 1:epochs
        action = let ϵ = ϵ; s -> act(model, s, ϵ, nA) end
        episodes!(action, memory, env, steps, maxt)

        for _ in 1:trainiters
            s1, a, r, s2, d = sample(memory)
            Qtar = r .+ (1 .- d) .* γ .* vec(maximum(target(s2), dims=1))
            gs = gradient(ps) do
                Qs1a = gather(model(s1), a)
                return Flux.mse(Qs1a, Qtar)
            end
            update!(opt, ps, gs)
            softupdate!(target, model, τ)
        end

        ϵ = max(ϵ * 0.99, 0.01)
        rewards[i] = test(s -> Flux.onecold(model(s)), env)
    end
    return rewards
end

env = CartPole{Float32}()

model = Chain(
    Dense(length(env.observationspace), 32, relu),
    Dense(32, length(env.actionspace))
)

r = run!(model, ADAM(), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "dqn")


text:
using FSM

cn = Constants{Float64}()

ebm = EBM{Float64}(
    am=1,
    cm=1,
    dm=1,
    em=1,
    hm=1,
)

Sf = 1.0

@code_warntype surf_props(ebm, cn, Sf)

time_surf_props(ebm, cn, Sf) = @time surf_props(ebm, cn, Sf)

time_surf_props(ebm, cn, Sf)
time_surf_props(ebm, cn, Sf)


text:
using EnsembleInference
using EnsembleInference: E1, E2, E3, representation_block
using LinearAlgebra
using Manifolds: Manifolds
using Rotations
using Serialization
using Test

isskewhermitian(A) = ishermitian(im * A)

isunitary(A) = A * A' ≈ I && abs(det(A)) ≈ 1

# generate and sparsify all blocks up to maximum degree L
generate_blocks(M, p, L) = [representation_block(M, p, ℓ) for ℓ in 0:L]

# test against reference file. if reference file doesn't exist, create it.
# used to ensure that future optimizations don't break correctness
function test_reference(fn, A)
    dir = Base.source_dir()
    path = abspath(joinpath(dir, fn))
    dirname, filename = splitdir(path)
    if !isfile(path)
        mkpath(dirname)
        open(io -> serialize(io, A), path, "w")
    end
    Aref = open(deserialize, path, "r")
    @test A ≈ Aref
end

# explicit little-d matrix entries, from Tables 4.3-4.6 of
# Varshalovich. Quantum Theory of Angular Momentum. 1998.
# NOTE: the tables order the indices in decreasing order, where we order them in increasing
# order.
function dmat(β, ℓ)
    sinβ, cosβ = sincos(β)
    α = β / 2
    sinα, cosα = sincos(α)
    return dℓ = if ℓ == 0
        [one(cosβ)]
    elseif ℓ == 1//2
        [cosα sinα; -sinα cosα]
    elseif ℓ == 1
        [
            (1 + cosβ)/2 sinβ/√2 (1 - cosβ)/2
            -sinβ/√2 cosβ sinβ/√2
            (1 - cosβ)/2 -sinβ/√2 (1 + cosβ)/2
        ]
    elseif ℓ == 3//2
        [
            cosα^3 √3*sinα*cosα^2 √3*sinα^2*cosα sinα^3
            -√3*sinα*cosα^2 cosα*(3cosα^2 - 2) -sinα*(3sinα^2 - 2) √3*sinα^2*cosα
            √3*sinα^2*cosα sinα*(3sinα^2 - 2) cosα*(3cosα^2 - 2) √3*sinα*cosα^2
            -sinα^3 √3*sinα^2*cosα -√3*sinα*cosα^2 cosα^3
        ]
    else
        throw(DomainError(ℓ, "ℓ must be less than or equal to 3/2"))
    end
end

@testset "Representations" begin
    @testset "so(3) representations" begin
        e = Matrix{Float64}(I, 3, 3)
        basis = Manifolds.DefaultOrthogonalBasis()
        @testset "basis vector E$i" for (i, Ei) in enumerate((E1, E2, E3))
            Xⁱ = Float64.(i .== (1:3))
            Eimat = Manifolds.get_vector(SO3, e, Xⁱ, basis)
            @test @inferred(representation_block(so3, Ei, 0)) ≈ zeros(1, 1)
            @testset "ℓ=$ℓ" for ℓ in 1:0.5:100
                u_Ei = @inferred representation_block(so3, Ei, ℓ)
                u_Eimat = @inferred representation_block(so3, Eimat, ℓ)
                if Ei === E3
                    @test u_Ei isa Diagonal
                else
                    @test u_Ei isa Tridiagonal
                end
                @test u_Eimat isa Tridiagonal
                @test isskewhermitian(u_Ei)
                @test isskewhermitian(u_Eimat)
                @test representation_block(so3, Ei, ℓ) ≈ representation_block(so3, Eimat, ℓ)
            end
        end
        # arrays serialized with v"1.5", deserialization with old Julia versions may fail
        VERSION ≥ v"1.5.3" && @testset "reference tests" begin
            test_reference("representations/so3_E1_30.jls", generate_blocks(so3, E1, 30))
            test_reference("representations/so3_E2_30.jls", generate_blocks(so3, E2, 30))
            test_reference("representations/so3_E3_30.jls", generate_blocks(so3, E3, 30))
        end
        @testset "non-basis vector" begin
            @testset "ℓ=$ℓ" for ℓ in 1:0.5:100
                Xⁱ = randn(3)
                X = Manifolds.get_vector(SO3, e, Xⁱ, basis)
                u1 = @inferred representation_block(so3, Xⁱ, ℓ)
                u2 = @inferred representation_block(so3, X, ℓ)
                u3 = (
                    Xⁱ[1] * representation_block(so3, E1, ℓ) +
                    Xⁱ[2] * representation_block(so3, E2, ℓ) +
                    Xⁱ[3] * representation_block(so3, E3, ℓ)
                )
                @test isskewhermitian(u3)
                @test u1 ≈ u3
                @test u2 ≈ u3
            end
        end
    end

    @testset "SO(3) representations" begin
        # R = Rotations.UnitQuaternion(normalize(randn(4)))
        @testset "basic group properties" begin
            @testset "ℓ=$ℓ" for ℓ in 0:20
                R1 = Rotations.UnitQuaternion(normalize(randn(4)))
                R2 = Rotations.UnitQuaternion(normalize(randn(4)))
                R12 = R1 * R2
                # identity
                Ue = @inferred representation_block(SO3, Matrix{Float64}(I, 3, 3), ℓ)
                @test Ue ≈ I
                # composition
                U1 = @inferred representation_block(SO3, R1, ℓ)
                @test isunitary(U1)
                U2 = @inferred representation_block(SO3, R2, ℓ)
                @test isunitary(U2)
                U12 = @inferred representation_block(SO3, R12, ℓ)
                @test isunitary(U12)
                @test U1 * U2 ≈ U12
                # invertible
                @test inv(U1) ≈ representation_block(SO3, inv(R1), ℓ)
                # don't check associativity, because implied by matrix
            end
        end

        @testset "little d matrix" begin
            @testset "ℓ=$ℓ" for ℓ in 0:0.5:1.5
                βs = randn(20)
                for β in βs
                    # Varshalovich uses the ZYZ Euler angle convention, see §1.4.1
                    R = Matrix(Rotations.RotY(β))
                    dℓ = representation_block(SO3, R, ℓ)
                    @test isunitary(dℓ)
                    @test dℓ ≈ dmat(β, ℓ)
                end
            end
        end

        # TODO: test to ensure conventions used are consistent with Varshalovich
        # TODO: use FiniteDifferences to test connection of algebra and group representations
    end
end


text:
using Zygote, OptimKit
using ..ITensorAutoHOOT
using ..ITensorNetworks
using ..ITensorAutoHOOT: batch_tensor_contraction
using ..ITensorNetworks: PEPS, inner_network, inner_networks, flatten, rayleigh_quotient

function loss_grad_wrap(peps::PEPS, Hs::Array)
  function loss(peps::PEPS)
    peps_bra = addtags(linkinds, peps, "bra")
    peps_ket = addtags(linkinds, peps, "ket")
    sites = commoninds(peps_bra, peps_ket)
    peps_ket_ham = prime(sites, peps_ket)
    network_H = inner_networks(peps_bra, peps_ket, peps_ket_ham, Hs)
    network_inner = inner_network(peps_bra, peps_ket)
    network_list = vcat(network_H, [network_inner])
    variables = flatten([peps_bra, peps_ket, peps_ket_ham])
    inners = batch_tensor_contraction(network_list, variables...)
    return rayleigh_quotient(inners)
  end
  loss_w_grad(peps::PEPS) = loss(peps), gradient(loss, peps)[1]
  return loss_w_grad
end

"""Update PEPS based on gradient descent
Parameters
----------
peps: a peps network with datatype PEPS
Hs: An array of MPO operators with datatype LocalMPO or LineMPO
stepsize: step size used in the gradient descent
num_sweeps: number of gradient descent sweeps/iterations
Returns
-------
An array containing Rayleigh quotient losses after each iteration.
"""
function gradient_descent(peps::PEPS, Hs::Array; stepsize::Float64, num_sweeps::Int)
  loss_w_grad = loss_grad_wrap(peps, Hs)
  return gradient_descent(peps, loss_w_grad; stepsize=stepsize, num_sweeps=num_sweeps)
end

function OptimKit.optimize(peps::PEPS, Hs::Array; num_sweeps::Int, method="GD")
  loss_w_grad = loss_grad_wrap(peps, Hs)
  return optimize(peps, loss_w_grad; num_sweeps=num_sweeps, method=method)
end


text:
# ---
# title: 208. Implement Trie (Prefix Tree)
# id: problem208
# author: zhwang
# date: 2022-03-17
# difficulty: Medium
# categories: Design, Trie
# link: <https://leetcode.com/problems/implement-trie-prefix-tree/description/>
# hidden: true
# ---
#
# Implement a trie with `insert`, `search`, and `startsWith` methods.
#
# **Example:**
#
#
#
#     Trie trie = new Trie();
#
#     trie.insert("apple");
#     trie.search("apple");   // returns true
#     trie.search("app");     // returns false
#     trie.startsWith("app"); // returns true
#     trie.insert("app");
#     trie.search("app");     // returns true
#
#
# **Note:**
#
#   * You may assume that all inputs are consist of lowercase letters `a-z`.
#   * All inputs are guaranteed to be non-empty strings.
#
#
## @lc code=start
using LeetCode

Base.@kwdef mutable struct PrefixNode
    isend::Bool = false
    children = Dict{Char, PrefixNode}()
end

function insert_node!(node::PrefixNode, word::String)::Nothing
    for c in word
        children = node.children
        haskey(children, c) || (children[c] = PrefixNode())
        node = children[c]
    end
    node.isend = true
    return
end

function search_prefix_node(node::PrefixNode, prefix::String)::Union{Nothing, PrefixNode}
    for c in prefix
        children = node.children
        haskey(children, c) || return
        node = children[c]
    end
    node
end

function search_word(node::PrefixNode, word::String)::Bool
    node = search_prefix_node(node, word)
    return !isnothing(node) && node.isend
end

starts_with(node::PrefixNode, prefix::String) = !isnothing(search_prefix_node(node, prefix))

## @lc code=end

## @lc test=start
@testset "208.implement-trie-prefix-tree.jl" begin
    @testset "208 case 1" begin
        node = PrefixNode()
        insert_node!(node, "apple")
        @test search_word(node, "apple")
        @test !search_word(node, "app")
        @test starts_with(node, "app")
        insert_node!(node, "app")
        @test search_word(node, "app")
    end
end
## @lc test=end

## @lc info=start
# link: [solution to 208]()
## @lc info=end


text:
function solve()
    n, C = [parse(Int, x) for x in split(readline())]
    event = []
    for i=1:n
        a, b, c = [parse(Int, x) for x in split(readline())]
        push!(event, (a-1, c))
        push!(event, (b, -c))
    end
    sort!(event)

    ret = 0
    fee = 0
    t = 0
    for (x, y) in event
        if x ≠ t
            ret += min(C, fee) * (x - t)
            t = x
        end
        fee += y
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia