text:
##################################################
using HDF5 # To have access to .hf5 files
########################################
include("../sources/Main.jl") # Loading the main code
########################################
# z=1.3
qminMeasure, qmaxMeasure = 0.01,1.0 # Range in q where the GR are computed
xminMeasure, xmaxMeasure = 0.0,0.4 # Range in x where the GR are computed

nbqMeasure = 10 # Number of q for which the GR are computed
nbxMeasure = 10 # Number of x for which the GR are computed
nbqxGrid = nbqMeasure*nbxMeasure
tabqMeasure = exp.(range(log(qminMeasure),length=nbqMeasure,log(qmaxMeasure)))
tabxMeasure = collect(range(xminMeasure,length=nbxMeasure,xmaxMeasure))

const tabqxGrid = zeros(Float64,2,nbqxGrid) # Location (q,x) of the grid points where the diffusion coefficients are computed
const tabGRGrid = zeros(Float64,nbqxGrid)
const tabRotGrid = zeros(Float64,nbqxGrid)

function tabqxGrid!()
    index = 1
    for iq=1:nbqMeasure
    qMeasure = tabqMeasure[iq]
        for ix=1:nbxMeasure
            xMeasure = tabxMeasure[ix]
            tabqxGrid[1,index], tabqxGrid[2,index] = qMeasure, xMeasure
            index += 1
        end
    end
end

function tabGRRotGrid!()

    if (PARALLEL == "yes") # Computation is made in parallel
        println("parallel")
        Threads.@threads for iGrid=1:nbqxGrid
            q, x = tabqxGrid[1,iGrid], tabqxGrid[2,iGrid]

            tabOmega_parallel = zeros(Float64,nbK)
            tabKappaSqOverTwoOmega_parallel = zeros(Float64,nbK)

            tabAln_parallel = zeros(Float64,N+1,N+1)
            tabBln_parallel = zeros(Float64,N+1,N+1)
            tabCln_parallel = zeros(Float64,N+1,N+1)
            tabDln_parallel = zeros(Float64,N+1,N+1)
            tabFln_parallel = zeros(Float64,N+1,N+1)
            tabGln_parallel = zeros(Float64,N+1,N+1)
            tabHln_parallel = zeros(Float64,N+1,N+1)

            tabTruncMln_parallel = Vector{Array{Float64,2}}([zeros(Float64,3*k,3*k) for k=1:N+1])
            tabEigValsMln_parallel = Vector{Vector{Complex{Float64}}}([zeros(Float64,k) for k=1:N+1])

            re, im = growth_rate_with_rotation(x,q,tabOmega_parallel,tabKappaSqOverTwoOmega_parallel,
                            tabAln_parallel,tabBln_parallel,tabCln_parallel,tabDln_parallel,
                            tabFln_parallel,tabGln_parallel,tabHln_parallel,tabTruncMln_parallel,
                            tabEigValsMln_parallel)

            if (im > threshold_GR)
                tabRotGrid[iGrid], tabGRGrid[iGrid] = re, im
            else
                tabRotGrid[iGrid], tabGRGrid[iGrid] = 0.0, 0.0
            end
        end
    else # Computation is not made in parallel
        println("not parallel")

        for iGrid=1:nbqxGrid # Loop over the elements of the (a,j)-grid
            q, x = tabqxGrid[1,iGrid], tabqxGrid[2,iGrid] # Current (a,j) location

            re, im = growth_rate_with_rotation(x,q)

            if (im > threshold_GR)
                tabRotGrid[iGrid], tabGRGrid[iGrid] = re, im
            else
                tabRotGrid[iGrid], tabGRGrid[iGrid] = 0.0, 0.0
            end
        end
    end
end

########################################
namefile = "../data/Dump_Growth_Rate_Precession.hf5"
########################################
# Function that writes to .hf5 files
########################################
function writedump!(namefile)
    file = h5open(namefile,"w") # Opening the file
    write(file,"a",a) # Characteristic length of disk
    write(file,"b",b) # Characteristic length of halo
    write(file,"c",c) # Characteristic length of bulge
    write(file,"eps",eps)
    write(file,"nbq",nbqMeasure) # Dumping the grid of (q,x)
    write(file,"nbx",nbxMeasure) # Dumping the grid of (q,x)
    write(file,"tabqx",tabqxGrid) # Dumping the grid of (q,x)
    write(file,"tabq",tabqMeasure) # Dumping the grid of q
    write(file,"tabx",tabxMeasure) # Dumping the grid of x
    write(file,"tabGR",tabGRGrid) # Dumping the total GR(q,x) for the current used value of lmax
    write(file,"tabRot",tabRotGrid) # Dumping the total Rot(q,x) for the current used value of lmax
    close(file) # Closing the file
end

########################################

tabqxGrid!()

println("Computing the maximum growth rate...")
@time tabGRRotGrid!()

########################################
println("Saving data...")
writedump!(namefile) # Dumping the computed table

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
"""
    root_to_leaves_ordering(mechanism; exclude_origin, exclude_loop_joints)

    ordered list of ids from root to leaves, all nodes are visited a single time
        excluding: origin, joints forming a loop which are not visited

    mechanism: Mechanism
    exclude_origin: flag to exclude origin from traversal
    exclude_loop_joints: flag to exclude loop joints from traversal
"""
function root_to_leaves_ordering(mechanism::Mechanism{T};
        exclude_origin::Bool=true,
        exclude_loop_joints::Bool=true) where T
    nodes = [mechanism.origin; mechanism.bodies; mechanism.joints; mechanism.contacts]
    loop_joints = get_loop_joints(mechanism.bodies, mechanism.joints)
    return root_to_leaves_ordering(nodes, loop_joints,
        exclude_origin=exclude_origin, exclude_loop_joints=exclude_loop_joints)
end

function root_to_leaves_ordering(nodes::Vector{Node{T}}, loop_joints;
        exclude_origin::Bool=true,
        exclude_loop_joints::Bool=true) where T
    ids = Vector{Int64}()
    stack = [0]
    while length(stack) > 0
        ids, stack = explore(ids, stack, nodes, loop_joints)
    end
    exclude_origin && (ids = ids[2:end]) # remove origin
    !exclude_loop_joints && push!(ids, getfield.(loop_joints, :id)...) # add loop_joints
    return ids
end

function explore(ids::Vector{Int}, stack::Vector{Int}, nodes::Vector{Node{T}}, loop_joints) where T
    loop_joints_ids = getfield.(loop_joints, :id)
    id = pop!(stack)
    push!(ids, id)
    node = nodes[findfirst(n -> n.id == id, nodes)]
    child_ids = get_child_ids(node, nodes)
    setdiff!(child_ids, loop_joints_ids)
    push!(stack, child_ids...)
    return ids, stack
end

function get_child_ids(joint::JointConstraint, nodes::Vector{Node{T}}) where T
    [joint.child_id]
end

function get_child_ids(node::Node, nodes::Vector{Node{T}}) where T
    child_ids = Vector{Int64}()
    for cnode in nodes
        if !(cnode isa Origin) && !(cnode isa Body)
            (cnode.parent_id == node.id) && push!(child_ids, cnode.id)
        end
    end
    return child_ids
end

function get_child_ids(body::Contact, nodes::Vector{Node{T}}) where T
    Vector{Int64}()
end

function get_loop_joints(bodies::Vector{<:Body}, joints::Vector{<:JointConstraint})
    # Assumes that the Origin cannot be the child of a joint.
    # the results is dependent on the ordering of the joints provided to the method.
    loop_joints = []
    body_ids = getfield.(bodies, :id)
    visited_bodies = zeros(Int,0)
    for joint in joints
        if joint.child_id in visited_bodies
            push!(loop_joints, joint)
        else
            push!(visited_bodies, joint.child_id)
        end
    end
    return [loop_joints...]
end

ext: jl


text:
using Spectra
using Test
using Random

# Dummy data
x = collect(0.:1.:100.)

# Our peak is located at freq_th with a half-width of hwhm_th
int_th = 10.
freq_th = 40.
hwhm_th = 10.

# we now generate a perfect y as well as a noisy y
y = int_th.*exp.(-log.(2) .*((x.-freq_th)./hwhm_th).^2)
Random.seed!(42); # we fix the seed
y_noise = y+randn((size(y,1),1))

# we measured the width and frequency with the peakhw function
int_mea1, freq_meas1, hwhm_meas1, centroid_mea1 = peakmeas(x,y,smoothing="no")
#int_mea2, freq_meas2, hwhm_meas2, centroid_mea2 = peakmeas(x,vec(y_noise),method="savgol",window_length=5,polyorder=2) # with applying a Savitsky Golay filter

# test of the perfect version
@test isapprox(int_th,int_mea1,atol=1e-5)
@test isapprox(freq_th,freq_meas1,atol=1e-5)
@test isapprox(hwhm_th,hwhm_meas1,atol=1e-5)

# test of the noisy version
#@test isapprox(int_th,int_mea2,atol=1.)
#@test isapprox(freq_th,freq_meas2,atol=1.)
#@test isapprox(hwhm_th,hwhm_meas2,atol=1.)

# Now we test 2 gaussians to test the centroid
int2_th = 20.00000 # the two intensities cannot be fully equal, as this will raise an error.
freq2_th = 60.
hwhm2_th = 10.

y2 = int2_th.*exp.(-log(2) .*((x.-freq2_th)./hwhm2_th).^2)
y_tot = y + y2

~, ~, ~, centroid_mea3 = peakmeas(x,y_tot,smoothing="no")

int_tot = int2_th + int_th

centroid_th = int_th./int_tot.*freq_th .+ int2_th./int_tot.*freq2_th # this works because the peaks have the same hwhm
@test isapprox(centroid_th,centroid_mea3,atol=1e-4)


text:
function RK4(previous)
    k1 = ode(previous)
    k2 = ode(previous + 0.5*dt*k1)
    k3 = ode(previous + 0.5*dt*k2)
    k4 = ode(previous + dt*k3)

    return previous + (dt/6)*(k1 + 2*k2 + 2*k3 + k4)
end

function ode(x)
    return (circshift(x, -1) - circshift(x, 2)) .* circshift(x, 1) - x .+ F
end

const F = 8.0
const dt = 0.005

nx = 50
state = zeros(nx)
state[1] = 8.0

for i in 1:20000
    state = RK4(state)
end

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
#Julia implementation of "load_network" function
export NNet, evaluate_network, evaluate_network_multiple, num_inputs
type NNet
    file::AbstractString
    weights::Array{Any,1}
    biases::Array{Any,1}
    symmetric::Int32
    numLayers::Int32
    inputSize::Int32
    outputSize::Int32
    maxLayerSize::Int32

    layerSizes::Array{Int32,1}
    mins::Array{Float64,1}
    maxes::Array{Float64,1}
    means::Array{Float64,1}
    ranges::Array{Float64,1}

    numNetworks::Int32
    numTau::Int32
    numPra::Int32
    tauCut::Array{Float64,1}
    praCut::Array{Float64,1}
    networkArray::Bool

    function NNet(file::AbstractString)
        this  = new()
        this.file = file
        f = open(this.file)

        ###############
        # Begin reading .nnet file
        ###############

        #Read header line
        line = readline(f)
        this.networkArray = false
        if line[1] == 'A'
            this.networkArray = true
        end

        if this.networkArray
            #Read number of neural networks
            line = readline(f)
            record = split(line,[',','\n'])
            this.numNetworks = parse(Int32,record[1])

            #Read number of previous RA cutpoints
            line = readline(f)
            record = split(line,[',','\n'])
            this.numPra = parse(Int32,record[1])

            #Read each previous RA cutpoint
            line = readline(f)
            record = split(line,[',','\n'])
            this.praCut = zeros(this.numPra)
            for i=1:(this.numPra)
                this.praCut[i] = parse(Float64,record[i])
            end

            #Read the number of tau cutpoints
            line = readline(f)
            record = split(line,[',','\n'])
            this.numTau = parse(Int32,record[1])

            #Read each tau cutpoint
            line = readline(f)
            record = split(line,[',','\n'])
            this.tauCut = zeros(this.numTau)
            for i=1:(this.numTau)
                this.tauCut[i] = parse(Float64,record[i])
            end

        else
            this.numNetworks = 1
        end

        #Read information about the neural network
        line = readline(f)
        record = split(line,[',','\n'])
        this.numLayers = parse(Int32,record[1])
        this.inputSize = parse(Int32,record[2])
        this.outputSize = parse(Int32,record[3])
        this.maxLayerSize=parse(Int32,record[4])

        line = readline(f)
        record = split(line,[',','\n'])
        this.layerSizes = zeros(this.numLayers+1)
        for i=1:(this.numLayers+1)
            this.layerSizes[i]=parse(Int32,record[i])
        end

        line = readline(f)
        record = split(line,[',','\n'])
        this.symmetric = parse(Int32,record[1])

        line = readline(f)
        record = split(line,[',','\n'])
        this.mins = zeros(this.inputSize)
        for i=1:(this.inputSize)
            this.mins[i]=parse(Float64,record[i])
        end

        line = readline(f)
        record = split(line,[',','\n'])
        this.maxes = zeros(this.inputSize)
        for i=1:(this.inputSize)
            this.maxes[i]=parse(Float64,record[i])
        end

        line = readline(f)
        record = split(line,[',','\n'])
        this.means = zeros(this.inputSize+1)
        for i=1:(this.inputSize+1)
            this.means[i]=parse(Float64,record[i])
        end

        line = readline(f)
        record = split(line,[',','\n'])
        this.ranges = zeros(this.inputSize+1)
        for i=1:(this.inputSize+1)
            this.ranges[i]=parse(Float64,record[i])
        end

        this.weights = Any[zeros(this.layerSizes[2],this.layerSizes[1])]
        this.biases  = Any[zeros(this.layerSizes[2])]
        for i=2:this.numLayers
            this.weights = [this.weights;Any[zeros(this.layerSizes[i+1],this.layerSizes[i])]]
            this.biases  = [this.biases;Any[zeros(this.layerSizes[i+1])]]
        end
        this.weights = Any[this.weights]
        this.biases  = Any[this.biases]
        for i=2:this.numNetworks
            w2 = Any[zeros(this.layerSizes[2],this.layerSizes[1])]
            b2 = Any[zeros(this.layerSizes[2])]
            for i=2:this.numLayers
                w2 = [w2;Any[zeros(this.layerSizes[i+1],this.layerSizes[i])]]
                b2 = [b2;Any[zeros(this.layerSizes[i+1])]]
            end
            this.weights = [this.weights;Any[w2]]
            this.biases  = [this.biases;Any[b2]]
        end

        nnetInd = 1
        layer=1
        i=1
        j=1
        line = readline(f)
        record = split(line,[',','\n'])
        while !eof(f)
            if layer>this.numLayers
                layer=1
                nnetInd=nnetInd+1
            end
            while i<=this.layerSizes[layer+1]
                while record[j]!=""
                    this.weights[nnetInd][layer][i,j] = parse(Float64,record[j])
                    j=j+1
                end
                j=1
                i=i+1
                line = readline(f)
                record = split(line,[',','\n'])
            end
            i=1
            while i<=this.layerSizes[layer+1]
                this.biases[nnetInd][layer][i] = parse(Float64,record[1])
                i=i+1
                line = readline(f)
                record = split(line,[',','\n'])
            end
            layer=layer+1
            i=1
            j=1
        end
        close(f)

        return this
    end
end

#Evaluates one set of inputs
function evaluate_network(nnet::NNet,input::Array{Float64,1})
    numLayers = nnet.numLayers
    inputSize = nnet.inputSize
    outputSize = nnet.outputSize
    symmetric = nnet.symmetric
    biases = nnet.biases
    weights = nnet.weights

    nnetInd = 1
    if nnet.numNetworks>1
        tau = input[6]
        pra = input[7]
        nnetInd += round(Int32,pra) - 2
        nnetInd *= nnet.numTau

        upper = nnet.numTau
        lower = 1
        while lower < upper
            middle = floor(Int32,(upper+lower)/2)
            if tau - nnet.tauCut[middle] > nnet.tauCut[middle+1]-tau
                lower = middle + 1
            else
                upper = middle
            end
        end
        nnetInd += lower
    end

    inputs = zeros(inputSize)
    for i = 1:inputSize
        if input[i]<nnet.mins[i]
            inputs[i] = (nnet.mins[i]-nnet.means[i])/nnet.ranges[i]
        elseif input[i] > nnet.maxes[i]
            inputs[i] = (nnet.maxes[i]-nnet.means[i])/nnet.ranges[i]
        else
            inputs[i] = (input[i]-nnet.means[i])/nnet.ranges[i]
        end
    end

    for layer = 1:numLayers-1
        temp = max.(weights[nnetInd][layer]*inputs[1:nnet.layerSizes[layer]]+biases[nnetInd][layer],0)
        inputs = temp
    end
    outputs = weights[nnetInd][end]*inputs[1:nnet.layerSizes[end-1]]+biases[nnetInd][end]
    for i=1:outputSize
        outputs[i] = outputs[i]*nnet.ranges[end]+nnet.means[end]
    end
    return outputs
end

#Evaluates multiple inputs at once. Each set of inputs should be a column in the input array
#Returns a column of output Q values for each input set
function evaluate_network_multiple(nnet::NNet,input::Array{Float64,2})
    numLayers = nnet.numLayers
    inputSize = nnet.inputSize
    outputSize = nnet.outputSize
    symmetric = nnet.symmetric
    biases = nnet.biases
    weights = nnet.weights

    _,numInputs = size(input)
    symmetryVec = zeros(numInputs)

    nnetInd = 1
    if nnet.numNetworks>1
        tau = input[6]
        pra = input[7]
        nnetInd += round(Int32,pra) - 2
        nnetInd *= nnet.numTau

        upper = nnet.numTau
        lower = 1
        while lower < upper
            middle = floor(Int32,(upper+lower)/2)
            if tau - nnet.tauCut[middle] > nnet.tauCut[middle+1]-tau
                lower = middle + 1
            else
                upper = middle
            end
        end
        nnetInd += lower
    end

    inputs = zeros(inputSize,numInputs)
    for i = 1:inputSize
        for j = 1:numInputs
            if input[i,j]<nnet.mins[i]
                inputs[i,j] = (nnet.mins[i]-nnet.means[i])/nnet.ranges[i]
            elseif input[i,j] > nnet.maxes[i]
                inputs[i,j] = (nnet.maxes[i]-nnet.means[i])/nnet.ranges[i]
            else
                inputs[i,j] = (input[i,j]-nnet.means[i])/nnet.ranges[i]
            end
        end
    end

    for layer = 1:numLayers-1
        inputs = max.(weights[nnetInd][layer]*inputs[1:nnet.layerSizes[layer],:]+biases[nnetInd][layer]*ones(1,numInputs),0)
    end
    outputs = weights[nnetInd][end]*inputs[1:nnet.layerSizes[end-1],:]+biases[nnetInd][end]*ones(1,numInputs)
    for i=1:outputSize
        for j=1:numInputs
            outputs[i,j] = outputs[i,j]*nnet.ranges[end]+nnet.means[end]
        end
    end
    return outputs
end

function num_inputs(nnet::NNet)
    if nnet.numNetworks==1
        return nnet.inputSize
    else
        return nnet.inputSize+2
    end
end
function num_outputs(nnet::NNet)
    return nnet.outputSize
end

max_stars_repo_name: alisamiii/VerticalCAS

max_issues_repo_name: alisamiii/VerticalCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
function create_drawing(img::Array{UInt32,2})
  return Cairo.CairoRGBSurface(img)
end

function create_contex(surface::Cairo.CairoSurface)
  return Cairo.CairoContext(surface)
end

function draw_vectors(ctx::Cairo.CairoContext, vectors, color, factor=1)
  Cairo.save(ctx)
  Cairo.set_source_rgb(ctx, color...)
  for line in vectors
    diff = line[2] - line[1]
    Cairo.move_to(ctx, line[1]...)
    Cairo.line_to(ctx, line[1]+diff*factor...)
  end
  Cairo.set_line_width(ctx, 3.0)
  Cairo.stroke(ctx)
  return ctx
end

function draw_line(ctx::Cairo.CairoContext, line, color=(0,0,0), thickness=2.0)
  Cairo.save(ctx)
  Cairo.set_source_rgb(ctx, color...)
  Cairo.move_to(ctx, line[1]...)
  Cairo.line_to(ctx, line[2]...)
  Cairo.set_line_width(ctx, thickness)
  Cairo.stroke(ctx)
  return ctx
end

function draw_points(ctx::Cairo.CairoContext, points, color=(0,0,0), factor=1)
  d = 5
  Cairo.save(ctx)
  Cairo.set_source_rgb(ctx, color...)
  for point in points
    Cairo.move_to(ctx, point+[-d,-d]...)
    Cairo.line_to(ctx, point+[d,d]...)
    Cairo.move_to(ctx, point+[-d,d]...)
    Cairo.line_to(ctx, point+[d,-d]...)
  end
  Cairo.set_line_width(ctx, 3.0)
  Cairo.stroke(ctx)
  return ctx
end

function draw_point(ctx::Cairo.CairoContext, point, color, factor=1)
  d = 5
  Cairo.save(ctx)
  Cairo.set_source_rgb(ctx, color...)
  Cairo.move_to(ctx, point+[-d,-d]...)
  Cairo.line_to(ctx, point+[d,d]...)
  Cairo.move_to(ctx, point+[-d,d]...)
  Cairo.line_to(ctx, point+[d,-d]...)
  Cairo.set_line_width(ctx, 2.0)
  Cairo.stroke(ctx)
  return ctx
end

function draw_text(ctx::Cairo.CairoContext, text, point, offset, fontsize, color)
  fill_color = color
  line_color = [1,1,1] - color
  Cairo.save(ctx)
  Cairo.select_font_face(ctx, "Sans", Cairo.FONT_SLANT_NORMAL,
                                          Cairo.FONT_WEIGHT_BOLD)
  Cairo.set_font_size(ctx, fontsize)
  Cairo.move_to(ctx, point+offset...)
  Cairo.rotate(ctx, -pi/2)
  Cairo.scale(ctx, -1, 1)
  Cairo.text_path(ctx, text)
  Cairo.set_source_rgb(ctx, fill_color...)
  Cairo.fill_preserve(ctx)
  Cairo.set_source_rgb(ctx, line_color...)
  Cairo.set_line_width(ctx, 1.0)
  Cairo.scale(ctx, -1, 1)
  Cairo.rotate(ctx, pi/2)
  Cairo.stroke(ctx)
  return ctx
end

function draw_indices(ctx::Cairo.CairoContext, points, offset, fontsize, color)
  fill_color = color
  line_color = [1,1,1] - color
  Cairo.save(ctx)
  Cairo.select_font_face(ctx, "Sans", Cairo.FONT_SLANT_NORMAL,
                                          Cairo.FONT_WEIGHT_BOLD)
  Cairo.set_font_size(ctx, fontsize)
  for (k, point) in enumerate(points)
    Cairo.move_to(ctx, point+offset...)
    Cairo.rotate(ctx, -pi/2)
    Cairo.scale(ctx, -1, 1)
    Cairo.text_path(ctx, "$k")
    Cairo.set_source_rgb(ctx, fill_color...)
    Cairo.fill_preserve(ctx)
    Cairo.set_source_rgb(ctx, line_color...)
    Cairo.set_line_width(ctx, 1.0)
    Cairo.scale(ctx, -1, 1)
    Cairo.rotate(ctx, pi/2)
  end
  Cairo.stroke(ctx)
  return ctx
end

function get_drawing(surface::Cairo.CairoSurface)
  return convert(Array{RGB24}, surface.data)
end

function draw_reference(ctx, h, w, factor)
  reference = 1
  margin = 10
  a = [h-margin, w-margin-reference*factor]
  b = [h-margin, w-margin]
  color = [0.9,0.9,0.9]
  draw_line(ctx, (a,b), color, 3.0)
  draw_text(ctx, "1", a+[-20, reference*factor/2], [0,0], 24.0, color)
  return ctx
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    laplacematrix(mesh; weights=:uniform, normalize=true)

The Laplace-Beltrami (a.k.a. Laplacian) matrix of the `mesh`.
Optionally specify the discretization `weights` as either
`:uniform` or `:cotangent` and `normalize` the rows by the
diagonal value.

## References

* Vallet, B & Lévy, B. 2008. [Spectral Geometry Processing with Manifold
  Harmonics](https://onlinelibrary.wiley.com/doi/10.1111/j.1467-8659.2008.01122.x)

* Zhang et al. 2007. [Spectral Methods for Mesh Processing and Analysis]
  (https://diglib.eg.org/handle/10.2312/egst.20071052.001-022)
"""
function laplacematrix(mesh; weights=:uniform, normalize=true)
  # convert to half-edge topology
  m = topoconvert(HalfEdgeTopology, mesh)

  # retrieve adjacency relation
  t = topology(m)
  𝒩 = Adjacency{0}(t)

  # initialize matrix
  n = nvertices(t)
  L = spzeros(n, n)

  # fill matrix with weights
  if weights == :uniform
    for i in 1:n
      js = 𝒩(i)
      for j in js
        L[i,j] = 1.0
      end
      L[i,i] = -1.0*length(js)
      if normalize
        for j in js
          L[i,j] /= -L[i,i]
        end
        L[i,i] /= -L[i,i]
      end
    end
  elseif weights == :cotangent
    v = vertices(m)
    for i in 1:n
      js = CircularVector(𝒩(i))
      for k in 1:length(js)
        j₋, j, j₊ = js[k-1], js[k], js[k+1]
        vᵢ, vⱼ =  v[i],  v[j]
        v₋, v₊ = v[j₋], v[j₊]
        αᵢⱼ = ∠(vⱼ, v₋, vᵢ)
        βᵢⱼ = ∠(vᵢ, v₊, vⱼ)
        L[i,j] = cot(αᵢⱼ) + cot(βᵢⱼ)
      end
      L[i,i] = -sum(L[i,js])
      if normalize
        for j in js
          L[i,j] /= -L[i,i]
        end
        L[i,i] /= -L[i,i]
      end
    end
  else
    throw(ArgumentError("invalid discretization weights"))
  end

  L
end


text:
using BenchmarkTools

SUITE = BenchmarkGroup()
SUITE["ISS"] = BenchmarkGroup()

include("iss.jl")

# ==============================
# Setup options
# ==============================

# general options
𝑂_iss = Options(:T=>time_horizon, :mode=>"check")

# algorithm-specific options
𝑂_dense_ISS01 = Options(:δ=>6e-4, :vars=>136:270, :assume_sparse=>true,
                        :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])
𝑂_discrete_ISS01 = Options(:discretization=>"nobloating", :δ=>0.01,
                           :vars=>136:270, :assume_sparse=>true,
                           :lazy_inputs_interval=>-1,
                           :partition=>[1:135, 136:270])

# ==============================
# ISU01 and ISS01
# ==============================

# specification options
𝑂_ISU01 = merge(𝑂_iss, Options(:property=>ISU01))
𝑂_ISS01 = merge(𝑂_iss, Options(:property=>ISS01))

# single run
sol_ISU01_dense = solve(iss_TV, 𝑂_ISU01, op=BFFPSV18(𝑂_dense_ISS01))
sol_ISS01_dense = solve(iss_TV, 𝑂_ISS01, op=BFFPSV18(𝑂_dense_ISS01))
sol_ISU01_discrete = solve(iss_TV, 𝑂_ISU01, op=BFFPSV18(𝑂_discrete_ISS01))
sol_ISS01_discrete = solve(iss_TV, 𝑂_ISS01, op=BFFPSV18(𝑂_discrete_ISS01))

# verify that specifications hold
@assert !sol_ISU01_dense.satisfied
@assert sol_ISS01_dense.satisfied
@assert !sol_ISU01_discrete.satisfied
@assert sol_ISS01_discrete.satisfied

# benchmark
SUITE["ISS"]["ISU01", "dense"] =
    @benchmarkable solve($iss_TV, $𝑂_ISU01, op=BFFPSV18($𝑂_dense_ISS01))
SUITE["ISS"]["ISS01", "dense"] =
    @benchmarkable solve($iss_TV, $𝑂_ISS01, op=BFFPSV18($𝑂_dense_ISS01))
SUITE["ISS"]["ISU01", "discrete"] =
    @benchmarkable solve($iss_TV, $𝑂_ISU01, op=BFFPSV18($𝑂_discrete_ISS01))
SUITE["ISS"]["ISS01", "discrete"] =
    @benchmarkable solve($iss_TV, $𝑂_ISS01, op=BFFPSV18($𝑂_discrete_ISS01))

# ==============================
# ISU02 and ISS02
# ==============================

# specification options
𝑂_ISU02 = merge(𝑂_iss, Options(:property=>ISU02))
𝑂_ISS02 = merge(𝑂_iss, Options(:property=>ISS02))

𝑂_dense_ISS02 = merge(𝑂_dense_ISS01,
                      Options(:δ=>5e-3, :partition=>[1:135, 136:270, 271:273]))
𝑂_discrete_ISS02 = merge(𝑂_discrete_ISS01,
                         Options(:partition=>[1:135, 136:270, 271:273]))

# single run
sol_ISU02_dense = solve(iss_CONST, 𝑂_ISU02, op=BFFPSV18(𝑂_dense_ISS02))
sol_ISS02_dense = solve(iss_CONST, 𝑂_ISS02, op=BFFPSV18(𝑂_dense_ISS02))
sol_ISU02_discrete = solve(iss_CONST, 𝑂_ISU02, op=BFFPSV18(𝑂_discrete_ISS02))
sol_ISS02_discrete = solve(iss_CONST, 𝑂_ISS02, op=BFFPSV18(𝑂_discrete_ISS02))

# verify that specifications hold
@assert !sol_ISU02_dense.satisfied
@assert sol_ISS02_dense.satisfied
@assert !sol_ISU02_discrete.satisfied
@assert sol_ISS02_discrete.satisfied

# benchmark
SUITE["ISS"]["ISU02", "dense"] =
    @benchmarkable solve($iss_CONST, $𝑂_ISU02, op=BFFPSV18($𝑂_dense_ISS02))
SUITE["ISS"]["ISS02", "dense"] =
    @benchmarkable solve($iss_CONST, $𝑂_ISS02, op=BFFPSV18($𝑂_dense_ISS02))
SUITE["ISS"]["ISU02", "discrete"] =
    @benchmarkable solve($iss_CONST, $𝑂_ISU02, op=BFFPSV18($𝑂_discrete_ISS02))
SUITE["ISS"]["ISS02", "discrete"] =
    @benchmarkable solve($iss_CONST, $𝑂_ISS02, op=BFFPSV18($𝑂_discrete_ISS02))

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
"""
    PoissonLikelihood(l::AbstractLink=ExpLink())

Poisson likelihood with rate defined as `l(f)`.

```math
    p(y|f) = Poisson(y | θ=l(f))
```

This is to be used if  we assume that the uncertainity associated
with the data follows a Poisson distribution.
"""
struct PoissonLikelihood{L<:AbstractLink}
    invlink::L
end

PoissonLikelihood() = PoissonLikelihood(ExpLink())

(l::PoissonLikelihood)(f::Real) = Poisson(l.invlink(f))

(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(Poisson.(l.invlink.(fs)))


text:
function trace_back(geometry::Geometry{T,S,A,E}, shaders::Vector{Function}, ray::Ray{T}, nBands::S, maxNumberOfCycles::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon for rendering
    count_while_loop = 0
    continue_iterating::Bool = true
    radiance::Vector{T} = ones(nBands) .* GREAT_NUM
    mtl::Int = 0
    hit::Bool = false
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfCycles
            continue_iterating = false
            continue
        end
        hit_, In, eidx, nrml, rmin = intersect(geometry, ray)
        if hit_ === true
	    hit = true
	    mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
	    if scatter === false
		continue_iterating = false
	    end
        else
	    radiance .*= 0.0
            continue_iterating = false
        end
    end
    return hit, radiance
end

function trace_back(geometry::Geometry{T,S,A,E}, geometry_bvh::Bvh{T,S,A,E}, shaders::Vector{Function}, ray::Ray{T}, nBands::S, maxNumberOfCycles::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon for rendering
    count_while_loop = 0
    continue_iterating::Bool = true
    radiance::Vector{T} = ones(nBands) .* GREAT_NUM
    mtl::Int = 0
    hit::Bool = false
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfCycles
            continue_iterating = false
            continue
        end
        hit_, In, eidx, nrml, rmin = intersect(geometry, geometry_bvh, ray)
        if hit_ === true
	    hit = true
	    mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
	    if scatter === false
		continue_iterating = false
	    end
        else
	    radiance .*= 0.0
            continue_iterating = false
        end
    end
    return hit, radiance
end

function trace_back(geometry::Geometry{T,S,A,E}, shaders::Vector{Function}, skymap::B, ray::Ray{T}, nBands::S, maxNumberOfCycles::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, B<:AbstractArray{T,2}, E<:AbstractArray{S,1}}
    # traces ray/photon for rendering
    count_while_loop = 0
    continue_iterating::Bool = true
    radiance::Vector{T} = ones(nBands) .* GREAT_NUM
    mtl::Int = 0
    hit::Bool = false
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfCycles
            continue_iterating = false
            continue
        end
        hit_, In, eidx, nrml, rmin = intersect(mesh, ray)
        if hit_mesh === true
	    hit = true
	    mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
	    if scatter === false
		continue_iterating = false
	    end
        else
	    radiance .*= sample_skymap(skymap, ray.dx, ray.dy, ray.dz)
            continue_iterating = false
        end
    end
    return hit, radiance
end

function trace_back(geometry::Geometry{T,S,A,E}, geometry_bvh::Bvh{T,S,A,E}, shaders::Vector{Function}, skymap::B, ray::Ray{T}, nBands::S, maxNumberOfCycles::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, B<:AbstractArray{T,2}, E<:AbstractArray{S,1}}
    # traces ray/photon for rendering
    count_while_loop = 0
    continue_iterating::Bool = true
    radiance::Vector{T} = ones(nBands) .* GREAT_NUM
    mtl::Int = 0
    hit::Bool = false
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfCycles
            continue_iterating = false
            continue
        end
        hit_, In, eidx, nrml, rmin = intersect(geometry, geometry_bvh, ray)
        if hit_ === true
	    hit = true
	    mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
	    if scatter === false
		continue_iterating = false
	    end
        else
	    radiance .*= sample_skymap(skymap, ray.dx, ray.dy, ray.dz)
            continue_iterating = false
        end
    end
    return hit, radiance
end

function trace_back(assets::AbstractArray{Asset{T,S},1}, geometries::AbstractArray{Geometry{T,S,A,E},1}, palettes::AbstractArray{Vector{Function},1}, geometry_bvhs::AbstractArray{Bvh{T,S,A,E},1}, scene_bvh::Bvh{T,S,A,E}, ray::Ray{T}, skymap::B, nBands::S, maxNumberOfCycles::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, B<:AbstractArray{T,2}, E<:AbstractArray{S,1}}
    # traces ray/photon for rendering
    count_while_loop::S = 0
    continue_iterating::Bool = true
    radiance::Vector{T} = ones(nBands)
    hit::Bool = false
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfCycles
            continue_iterating = false
            continue
        end
        hit_scene, I_scene, tray, nrml, range_scene, aidx, oidx, eidx = intersect_scene(assets, geometries, geometry_bvhs, scene_bvh, ray)
        if hit_scene === true
	    hit = true
            mtl = geometries[oidx].mtl[eidx]
            pidx = assets[aidx].pidx
            shader = palettes[pidx][mtl]
            scatter, ray, r = shader(I_scene.x, I_scene.y, I_scene.z, -tray.dx, -tray.dy, -tray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
	    if scatter === false
		continue_iterating = false
	    end
        else
	    radiance = radiance .* sample_skymap(skymap, ray.dx, ray.dy, ray.dz)
            continue_iterating = false
        end
    end
    return radiance
end

# following traces in forward direction with cyclic Walls.
# These functions record the "faterecs" and don't use skymaps.
# These function can be used to compute e.g. fAPAR or surface BDRF/BRF/BHR etc.
# but some are relatively slow because of the instantiation of small arrays...

function record(geometry::Geometry{T,S,A,E}, shaders::Vector{Function}, walls::Walls{T}, ray::Ray{T}, maxNumberOfBounces::Int) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon until it is absorbed or bounces through the reference plane
    #
    # TERMINAL CODES:
    #
    # -1    the maximum number of cycles was reached
    #  1    photon escapes through the reference plane
    #  2    photon escapes through the bottom plane
    #  3    photon cyclicly re-enters the scene (through walls)'
    #  4    photon is absorbed
    #  5    photon is scattered
    #
    path::AbstractArray{T, 2} = zeros(Float64, maxNumberOfBounces+1, 3)
    path[1,:] = [ray.x, ray.y, ray.z]
    codes::AbstractArray{Int, 1} = zeros(Int, maxNumberOfBounces+1) # code[1] === 0
    index::Int = 1
    count_while_loop::Int = 0
    continue_iterating::Bool = true
    mtl::Int = -1
    while continue_iterating
        count_while_loop += 1
        index += 1
        if count_while_loop >= maxNumberOfBounces
            path[index, :] = [ray.x, ray.y, ray.z]
            codes[index] = -1
            continue_iterating = false
            continue
        end
        I_wall, side, r_wall = intersect_walls(walls, ray)
        hit_, In, eidx, nrml, rmin = intersect(geometry, ray)
        wall = r_wall <= rmin
        if wall === true
            if side === 't'
                # photon escapes through reference plane
                continue_iterating = false
                codes[index] = 1
                path[index, :] = [I_wall.x, I_wall.y, I_wall.z]
            elseif (side === 'b')
                # photon is absorbed by a black surface underneath the scene
                continue_iterating = false
                codes[index] = 2
                path[index, :] = [I_wall.x, I_wall.y, I_wall.z]
            else # thus the condition: ' !(side ==='t') && !(side === 'b') ' is true
                # photon cyclicly re-enters from the opposite wall
                ray = cycle_ray(walls, ray)
                count_while_loop -= 1 # don't count these events toward the number of cycles
                index -= 1
            end
        else # implies hit_ is true, because the ray either hits the scene or a wall
            hit = true
            mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
            if scatter === false
                continue_iterating = false
                codes[index] = 4
                path[index, :] = [I_scene.x, I_scene.y, I_scene.z]
            else
                codes[index] = 5
                path[index, :] = [I_scene.x, I_scene.y, I_scene.z]
            end
	end
    end
    return path, codes
end

function record(geometry::Geometry{T,S,A,E}, geometry_bvh::Bvh{T,S,A,E}, shaders::Vector{Function}, walls::Walls{T}, ray::Ray{T}, maxNumberOfBounces::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon until it is absorbed or bounces through the reference plane
    #
    # TERMINAL CODES:
    #
    # -1    the maximum number of cycles was reached
    #  1    photon escapes through the reference plane
    #  2    photon escapes through the bottom plane
    #  3    photon cyclicly re-enters the scene (through walls)'
    #  4    photon is absorbed
    #  5    photon is scattered
    #
    path::AbstractArray{T, 2} = zeros(Float64, maxNumberOfBounces+1, 3)
    path[1,:] = [ray.x, ray.y, ray.z]
    codes::AbstractArray{Int, 1} = zeros(Int, maxNumberOfBounces+1) # code[1] === 0
    index::Int = 1
    count_while_loop::Int = 0
    continue_iterating::Bool = true
    mtl::Int = -1
    while continue_iterating
        count_while_loop += 1
        index += 1
        if count_while_loop >= maxNumberOfBounces
            path[index, :] = [ray.x, ray.y, ray.z]
            codes[index] = -1
            continue_iterating = false
            continue
        end
        I_wall, side, r_wall = intersect_walls(walls, ray)
        hit_, In, eidx, nrml, rmin = intersect(geometry, geometry_bvh, ray)
        wall = r_wall <= rmin
        if wall === true
            if side === 't'
                # photon escapes through reference plane
                continue_iterating = false
                codes[index] = 1
                path[index, :] = [I_wall.x, I_wall.y, I_wall.z]
            elseif (side === 'b')
                # photon is absorbed by a black surface underneath the scene
                continue_iterating = false
                codes[index] = 2
                path[index, :] = [I_wall.x, I_wall.y, I_wall.z]
            else # thus the condition: ' !(side ==='t') && !(side === 'b') ' is true
                # photon cyclicly re-enters from the opposite wall
                ray = cycle_ray(walls, ray)
                count_while_loop -= 1 # don't count these events toward the number of cycles
                index -= 1
            end
        else # implies hit_ is true, because the ray either hits the scene or a wall
            hit = true
            mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
            if scatter === false
                continue_iterating = false
                codes[index] = 4
                path[index, :] = [I_scene.x, I_scene.y, I_scene.z]
            else
                codes[index] = 5
                path[index, :] = [I_scene.x, I_scene.y, I_scene.z]
            end
	end
    end
    return path, codes
end

function record(assets::AbstractArray{Asset{T,S},1}, geometries::AbstractArray{Geometry{T,S,A,E},1}, palettes::AbstractArray{Vector{Function},1}, geometry_bvhs::AbstractArray{Bvh{T,S,A,E},1}, scene_bvh::Bvh{T,S,A,E}, walls::Walls{T}, ray::Ray{T}, maxNumberOfBounces::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon until it is absorbed or bounces through the reference plane
    #
    # TERMINAL CODES:
    #
    # -1    the maximum number of cycles was reached
    #  1    photon escapes through the reference plane
    #  2    photon escapes through the bottom plane
    #  3    photon cyclicly re-enters the scene (through walls)'
    #  4    photon is absorbed
    #  5    photon is scattered
    #
    path::AbstractArray{T, 2} = zeros(Float64, maxNumberOfBounces+1, 3)
    path[1,:] = [ray.x, ray.y, ray.z]
    codes::AbstractArray{Int, 1} = zeros(Int, maxNumberOfBounces+1) # code[1] === 0
    index::Int = 1
    count_while_loop::Int = 0
    continue_iterating::Bool = true
    mtl::Int = -1
    while continue_iterating
        count_while_loop += 1
	index += 1
        if count_while_loop >= maxNumberOfBounces
            path[index, :] = [ray.x, ray.y, ray.z]
            codes[index] = -1
            continue_iterating = false
            continue
        end
        I_wall, side, r_wall = intersect_walls(walls, ray)
        hit_scene, I_scene, tray, nrml, range_scene, aidx, oidx, eidx = intersect_scene(assets, geometries, geometry_bvhs, scene_bvh, ray)
        wall = r_wall <= range_scene
        if wall === true
            if side === 't'
                # photon escapes through reference plane
                continue_iterating = false
                codes[index] = 1
                path[index, :] = [I_wall.x, I_wall.y, I_wall.z]
            elseif (side === 'b')
                # photon is absorbed by a black surface underneath the scene
                continue_iterating = false
                codes[index] = 2
                path[index, :] = [I_wall.x, I_wall.y, I_wall.z]
            else # thus the condition: ' !(side ==='t') && !(side === 'b') ' is true
                # photon cyclicly re-enters from the opposite wall
                ray = cycle_ray(walls, ray)
                count_while_loop -= 1 # don't count these events toward the number of cycles
                index -= 1
            end
        else # wall is false, hence ray must have hit the (objects in the) scene
            hit = true
            @assert(hit_scene)
	    geometry = geometries[oidx]
            mtl = geometry.mtl[eidx]
            pidx = assets[aidx].pidx
            shader = palettes[pidx][mtl]
            scatter, ray, r = shader(I_scene.x, I_scene.y, I_scene.z, -tray.dx, -tray.dy, -tray.dz, nrml.x, nrml.y, nrml.z)
            if scatter === false
                continue_iterating = false
                codes[index] = 4
                path[index, :] = [I_scene.x, I_scene.y, I_scene.z]
            else
                codes[index] = 5
                path[index, :] = [I_scene.x, I_scene.y, I_scene.z]
            end
	end
    end
    return path, codes
end

function trace_forward(geometry::Geometry{T,S,A,E}, shaders::Vector{Function}, walls::Walls{T}, ray::Ray{T}, maxNumberOfBounces::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon until it is absorbed or bounces through the reference plane
    count_while_loop::Int = 0
    continue_iterating::Bool = true
    mtl::Int = -1
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfBounces
            continue_iterating = false
            continue
        end
        I_wall, side, r_wall = intersect_walls(walls, ray)
        hit_, In, eidx, nrml, rmin = intersect(geometry, ray)
        wall = r_wall <= rmin
        if wall === true
            if side === 't'
                # photon escapes through reference plane
                continue_iterating = false
            elseif (side === 'b')
                # photon is absorbed by a black surface underneath the scene
                continue_iterating = false
            else # thus the condition: ' !(side ==='t') && !(side === 'b') ' is true
                # photon cyclicly re-enters from the opposite wall
                ray = cycle_ray(walls, ray)
                count_while_loop -= 1 # don't count these events toward the number of cycles
            end
        else # implies hit_ is true, because the ray either hits the scene or a wall
            hit = true
            mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
            if scatter === false
                continue_iterating = false
            end
	end
    end
    return ray
end

function trace_forward(geometry::Geometry{T,S,A,E}, geometry_bvh::Bvh{T,S,A,E}, shaders::Vector{Function}, walls::Walls{T}, ray::Ray{T}, maxNumberOfBounces::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon until it is absorbed or bounces through the reference plane
    count_while_loop::Int = 0
    continue_iterating::Bool = true
    mtl::Int = -1
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfBounces
            continue_iterating = false
            continue
        end
        I_wall, side, r_wall = intersect_walls(walls, ray)
        hit_, In, eidx, nrml, rmin = intersect(geometry, geometry_bvh, ray)
        wall = r_wall <= rmin
        if wall === true
            if side === 't'
                # photon escapes through reference plane
                continue_iterating = false
            elseif (side === 'b')
                # photon is absorbed by a black surface underneath the scene
                continue_iterating = false
            else # thus the condition: ' !(side ==='t') && !(side === 'b') ' is true
                # photon cyclicly re-enters from the opposite wall
                ray = cycle_ray(walls, ray)
                count_while_loop -= 1 # don't count these events toward the number of cycles
            end
        else # implies hit_ is true, because the ray either hits the scene or a wall
            hit = true
            mtl = geometry.mtl[eidx]
            scatter, ray, r = shaders[mtl](In.x, In.y, In.z, -ray.dx, -ray.dy, -ray.dz, nrml.x, nrml.y, nrml.z)
            radiance .*= r
            if scatter === false
                continue_iterating = false
            end
	end
    end
    return ray
end

function trace_forward(assets::AbstractArray{Asset{T,S},1}, geometries::AbstractArray{Geometry{T,S,A,E},1}, palettes::AbstractArray{Vector{Function},1}, geometry_bvhs::AbstractArray{Bvh{T,S,A,E},1}, scene_bvh::Bvh{T,S,A,E}, walls::Walls{T}, ray::Ray{T}, maxNumberOfBounces::S) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    # traces ray/photon until it is absorbed or bounces through the reference plane
    count_while_loop::Int = 0
    continue_iterating::Bool = true
    mtl::Int = -1
    while continue_iterating
        count_while_loop += 1
        if count_while_loop >= maxNumberOfBounces
            push!(codes, -1)
            push!(path, Line(ray.x, ray.y, ray.z, ray.x, ray.y, ray.z)) # add a zero-length line, basically indicating that we stop tracing here...
            continue_iterating = false
            continue
        end
        I_wall, side, r_wall = intersect_walls(walls, ray)
        hit_scene, I_scene, tray, nrml, range_scene, aidx, oidx, eidx = intersect_scene(assets, geometries, geometry_bvhs, scene_bvh, ray)
        wall = r_wall <= range_scene
        if wall === true
            if side === 't'
                # photon escapes through reference plane
                continue_iterating = false
            elseif (side === 'b')
                # photon is absorbed by a black surface underneath the scene
                continue_iterating = false
            else # thus the condition: ' !(side ==='t') && !(side === 'b') ' is true
                # photon cyclicly re-enters from the opposite wall
                ray = cycle_ray(walls, ray)
                count_while_loop -= 1 # don't count these events toward the number of cycles
            end
        else # wall is false, hence ray must have hit the (objects in the) scene
            hit = true
	    geometry = geometries[oidx]
            mtl = geometry.mtl[eidx]
            pidx = assets[aidx].pidx
            shader = palettes[pidx][mtl]
            scatter, ray, r = shader(I_scene.x, I_scene.y, I_scene.z, -tray.dx, -tray.dy, -tray.dz, nrml.x, nrml.y, nrml.z)
            if scatter === false
                continue_iterating = false
            end
	end
    end
    return ray
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
function independent_dijkstra(mapf::MAPF, edge_weights::AbstractVector=mapf.edge_weights)
    @assert minimum(edge_weights) >= -eps()

    dijkstra_states = Dict(
        d => my_dijkstra_shortest_paths(
            mapf.rev_graph,
            d;
            edge_indices=mapf.rev_edge_indices,
            edge_weights=edge_weights,
        ) for d in unique(mapf.destinations)
    )

    A = nb_agents(mapf)
    solution = Vector{Path}(undef, A)
    for a in 1:A
        s, d, t0 = mapf.sources[a], mapf.destinations[a], mapf.starting_times[a]
        parents = dijkstra_states[d].parents
        t, v = t0, s
        path = [(t, v)]
        while v != d
            v = parents[v]
            t += 1
            push!(path, (t, v))
        end
        solution[a] = path
    end

    return solution
end

function independent_dijkstra(mapf::MAPF, edge_weights::AbstractMatrix)
    @assert minimum(edge_weights) >= -eps()
    A = nb_agents(mapf)
    solution = Vector{Path}(undef, A)
    for a in 1:A
        s, d, t0 = mapf.sources[a], mapf.destinations[a], mapf.starting_times[a]
        dijkstra_state = my_dijkstra_shortest_paths(
            mapf.rev_graph,
            d;
            edge_indices=mapf.rev_edge_indices,
            edge_weights=view(edge_weights, :, a),
        )
        parents = dijkstra_state.parents
        t, v = t0, s
        path = [(t, v)]
        while v != d
            v = parents[v]
            t += 1
            push!(path, (t, v))
        end
        solution[a] = path
    end

    return solution
end

function independent_topological_sort(
    mapf::MAPF, edge_weights::AbstractVector=mapf.edge_weights; T::Integer
)
    rev_g = mapf.rev_graph

    V = nv(rev_g)
    A = nb_agents(mapf)

    solution = Vector{Path}(undef, A)
    dists = OffsetMatrix{Float64}(undef, 0:T, V)
    parents = Matrix{Int}(undef, T, V)

    for d in unique(mapf.destinations)
        for v in vertices(rev_g)
            dists[0, v] = v == d ? 0.0 : Inf
        end
        for t in 1:T, v in 1:V
            dists[t, v] = Inf
            parents[t, v] = 0
            for u in inneighbors(rev_g, v)
                dist_du = dists[t - 1, u]
                weight_uv = edge_weights[mapf.rev_edge_indices[u, v]]
                if dist_du + weight_uv < dists[t, v]
                    dists[t, v] = dist_du + weight_uv
                    parents[t, v] = u
                end
            end
        end
        for a in 1:A
            mapf.destinations[a] == d || continue
            s = mapf.sources[a]
            t0 = mapf.starting_times[a]
            tf = minimum(t for t = 1:T if dists[t, s] < Inf)
            t, v = tf, s
            path = [(t0 + tf - t, v)]
            while v != d
                v = parents[t, v]
                t -= 1
                push!(path, (t0 + tf - t, v))
            end
            solution[a] = path
        end
    end
    return solution
end

function independent_astar(mapf::MAPF, edge_weights::AbstractVector=mapf.edge_weights)
    A = nb_agents(mapf)
    solution = Vector{Path}(undef, A)
    for a in 1:A
        s, d, t0 = mapf.sources[a], mapf.destinations[a], mapf.starting_times[a]
        dist = mapf.distances_to_destinations[d]
        heuristic(v) = dist[v]
        solution[a] = temporal_astar(
            mapf.graph,
            s,
            d,
            t0;
            edge_indices=mapf.edge_indices,
            edge_weights=edge_weights,
            heuristic=heuristic,
        )
    end
    return solution
end

function independent_astar(
    mapf::MAPF, constraints, edge_weights::AbstractVector=mapf.edge_weights
)
    A = nb_agents(mapf)
    solution = Vector{Path}(undef, A)
    for a in 1:A
        s, d, t0 = mapf.sources[a], mapf.destinations[a], mapf.starting_times[a]
        dist = mapf.distances_to_destinations[d]
        heuristic(v) = dist[v]
        forbidden_vertices = constraints[a]
        solution[a] = temporal_astar(
            mapf.graph,
            s,
            d,
            t0;
            edge_indices=mapf.edge_indices,
            edge_weights=edge_weights,
            forbidden_vertices=forbidden_vertices,
            heuristic=heuristic,
        )
    end
    return solution
end


text:
"""
    discretediag(chains::Chains{<:Real}; sections, kwargs...)

Discrete diagnostic where `method` can be
`[:weiss, :hangartner, :DARBOOT, MCBOOT, :billinsgley, :billingsleyBOOT]`.
"""
function MCMCDiagnosticTools.discretediag(
    chains::Chains{<:Real};
    sections = _default_sections(chains),
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Compute statistics.
    between_chain_vals, within_chain_vals = MCMCDiagnosticTools.discretediag(
        _chains.value.data; kwargs...
    )

    # Create dataframes
    parameters = (parameters = names(_chains),)
    between_chain_df = ChainDataFrame(
        "Chisq diagnostic - Between chains", merge(parameters, between_chain_vals),
    )
    within_chain_dfs = map(1:size(_chains, 3)) do i
        vals = map(val -> val[:, i], within_chain_vals)
        return ChainDataFrame("Chisq diagnostic - Chain $i", merge(parameters, vals))
    end
    dfs = vcat(between_chain_df, within_chain_dfs)

    return dfs
end


text:
    function PLacExampleUtility(ins)

        # Definition of the inputs for the ODEs
        if length(ins) == 1
            IPTG2,IPTG3 = ins[1];
        else
            IPTG2,IPTG3 = ins;
        end
        inputsMC = [0,IPTG2,IPTG3];

        # Solve ODEs
        solMC = PLacExample_SolveAll(tsMC, pD1MC, spMC, inputsMC, ivss1MC, sampsMC, pre1MC);

        # Extracte wanted vectors (observables) with time reduction
        Obs1_MC = solMC[:,2,:];
    Obs2_MC = 3 .*solMC[:,4,:];

        EntObs1_MC = zeros(size(solMC)[1]);
    EntObs2_MC = zeros(size(solMC)[1]);

        dists = [Beta, Exponential, LogNormal, Normal, Gamma, Laplace, Pareto, Rayleigh, Cauchy, Uniform];
        dists2 = [Exponential, LogNormal, Normal, Gamma, Laplace, Pareto, Rayleigh, Cauchy, Uniform];

        names = ["timePoint"];

        for j in 1:size(solMC)[1]
            fitts1 = Dict();
        bestfit1 = Dict();
        bestfitInd1 = zeros(1,1);
        fitts2 = Dict();
        bestfit2 = Dict();
        bestfitInd2 = zeros(1,1);

            try
                fitts1[names[1]] = fit.(dists, Ref(Obs1_MC[j,:]));
                bestfitInd1[1] = findmax(loglikelihood.(fitts1[names[1]], Ref(Obs1_MC[j,:])))[2];
                bestfit1[names[1]] = fitts1[names[1]][findmax(loglikelihood.(fitts1[names[1]], Ref(Obs1_MC[j,:])))[2]];
            catch
                fitts1[names[1]] = fit.(dists2, Ref(Obs1_MC[j,:]));
                bestfitInd1[1] = findmax(loglikelihood.(fitts1[names[1]], Ref(Obs1_MC[j,:])))[2]+1;
                bestfit1[names[1]] = fitts1[names[1]][findmax(loglikelihood.(fitts1[names[1]], Ref(Obs1_MC[j,:])))[2]];
            end
            try
                fitts2[names[1]] = fit.(dists, Ref(Obs2_MC[j,:]));
                bestfitInd2[1] = findmax(loglikelihood.(fitts2[names[1]], Ref(Obs2_MC[j,:])))[2];
                bestfit2[names[1]] = fitts2[names[1]][findmax(loglikelihood.(fitts2[names[1]], Ref(Obs2_MC[j,:])))[2]];
            catch
                fitts2[names[1]] = fit.(dists2, Ref(Obs2_MC[j,:]));
                bestfitInd2[1] = findmax(loglikelihood.(fitts2[names[1]], Ref(Obs2_MC[j,:])))[2]+1;
                bestfit2[names[1]] = fitts2[names[1]][findmax(loglikelihood.(fitts2[names[1]], Ref(Obs2_MC[j,:])))[2]];
            end

            EntObs1_MC[j] = entropy(bestfit1[names[1]]);
    EntObs2_MC[j] = entropy(bestfit2[names[1]]);

        end

        HES = zeros(1,2);
        HES[1,1] = sum(EntObs1_MC);
    HES[1,2] = sum(EntObs2_MC);

        util = sum(HES)/2

        return(util)
    end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
import MPI

export distribute, MPI_rank, MPI_size

if !MPI.Initialized()
    MPI.Init()
end

MPI_rank = MPI.Comm_rank(MPI.COMM_WORLD)
MPI_size = MPI.Comm_size(MPI.COMM_WORLD)

function distribute(size, comm_size, rank)
    """
    Compute the first element and size for a given rank

    Parameters
    ----------
    size : Int
        Length of array
    comm_size : Int
        Number of MPI  processes
    rank: Int
        MPI rank

    Returns
    -------
    size : Int
        Size
    offset : Int
        Offsets for the MPI process
    """
    if comm_size > size
        error("comm_size is larger than size!")
    end

    base = div(size, comm_size)
    leftover = size % comm_size
    if rank+1 <= leftover
        size = base+1
        start = rank * (base+1) + 1
    else
        size = base
        start = (base+1) * leftover + (rank-leftover) * base + 1
    end
    return start, size
end


text:
# This file includes and example file for 2D integral

using FractalTools
using GeometryBasics
using Makie

# Generate data
f(x, y) = [
    x^2 + y^2 + 1,
    x^2 - y^2 + 1,
]
ngon = Triangle(
    Point(BigFloat(0.), BigFloat(0.)),
    Point(BigFloat(1.), BigFloat(0.)),
    Point(BigFloat(0.5), BigFloat(1.)))
npts = 100
pts = getdata(f, ngon, npts)
interp = integrate(pts, HInterp2D(0.001 * ones(2,2)))


text:
using BlockArrays, BandedMatrices, BlockBandedMatrices, FillArrays, ArrayLayouts, LinearAlgebra, Test
import BlockBandedMatrices: MemoryLayout, ColumnMajor, BroadcastStyle, BlockBandedStyle, blockrowsupport, blockcolsupport
import Base.Broadcast: materialize!

@testset "BlockBandedMatrix" begin
    @test BroadcastStyle(BlockBandedMatrix) == BlockBandedStyle()

    @testset "BlockBandedMatrix constructors" begin
        l , u = 1,1
        N = M = 4
        cols = rows = 1:N

        @test Matrix(BlockBandedMatrix(Zeros(sum(rows),sum(cols)), rows, cols, (l,u))) ==
            Array(BlockBandedMatrix(Zeros(sum(rows),sum(cols)), rows, cols, (l,u))) ==
            zeros(Float64, 10, 10)

        @test Matrix(BlockBandedMatrix{Int}(Zeros(sum(rows),sum(cols)), rows,cols, (l,u))) ==
            zeros(Int, 10, 10)

        @test Matrix(BlockBandedMatrix(Eye(sum(rows)), rows,cols, (l,u))) ==
            Matrix{Float64}(I, 10, 10)

        @test Matrix(BlockBandedMatrix{Int}(Eye(sum(rows)), rows,cols, (l,u))) ==
            Matrix{Int}(I, 10, 10)

        @test Matrix(BlockBandedMatrix(I, rows,cols, (l,u))) ==
            Matrix{Float64}(I, 10, 10)

        @test Matrix(BlockBandedMatrix{Int}(I, rows,cols, (l,u))) ==
            Matrix{Int}(I, 10, 10)
    end

    @testset "BlockBandedMatrix block indexing" begin
        l , u = 1,1
        N = M = 4
        cols = rows = 1:N
        A = BlockBandedMatrix{Int}(undef, rows,cols, (l,u))
        A.data .= 1:length(A.data)

        @test A[1,1] == 1
        @test A[1,3] == 10

        @test blockbandwidth(A,1)  == 1
        @test blockbandwidths(A) == (l,u)

        # check views of blocks are indexing correctly

        @test A[Block(1), Block(1)] isa Matrix
        @test A[Block(1), Block(1)] == A[Block(1,1)] == BlockArrays.getblock(A, 1, 1) == Matrix(view(A, Block(1,1)))
        @test A[1,1] == view(A,Block(1),Block(1))[1,1] == view(A,Block(1,1))[1,1] == A[Block(1,1)][1,1]  == A[Block(1),Block(1)][1,1] == 1
        @test A[2,1] == view(A,Block(2),Block(1))[1,1] == view(A,Block(2,1))[1,1] == 2
        @test A[3,1] == view(A,Block(2),Block(1))[2,1] == 3
        @test A[4,1] == 0
        @test A[1,2] == view(A,Block(1,2))[1,1] == 4
        @test A[1,3] == view(A,Block(1,2))[1,2] == view(A,Block(1,2))[2] == 10

        @test view(A, Block(3),Block(1)) ≈ [0,0,0]
        @test_throws BandError view(A, Block(3),Block(1))[1,1] = 4
        @test_throws BlockBoundsError view(A, Block(5,1))

        # test blocks
        V = view(A, Block(1,1))
        @test_throws BoundsError V[2,1]

        V = view(A, Block(3,4))
        @test V[3,1] == 45
        V[3,1] = -7
        @test V[3,1] == -7
        @test Matrix(V) isa Matrix{Int}
        @test Matrix{Float64}(V) isa Matrix{Float64}
        @test Matrix{Float64}(Matrix(V)) == Matrix{Float64}(V)
        @test A[4:6,7:10] ≈ Matrix(V)

        A[1,1] = -5
        @test A[1,1] == -5
        A[1,3] = -6
        @test A[1,3] == -6

        A[Block(3,4)] = Matrix(Ones{Int}(3,4))
        @test A[Block(3,4)] == Matrix(Ones{Int}(3,4))

        ## Bug in setindex!
        ret = BlockBandedMatrix(Zeros{Float64}((4,6)), [2,2], [2,2,2], (0,2))

        V = view(ret, Block(1), Block(2))
        V[1,1] = 2
        @test ret[1,2] == 0
    end

    @testset "blockcol/rowsupport" begin
        l , u = 1,2
        N = M = 4
        cols = rows = 1:N
        A = BlockBandedMatrix{Int}(undef, rows,cols, (l,u))
        A.data .= 1:length(A.data)

        @test blockrowsupport(A, 1) == Block.(1:3)
        @test blockrowsupport(A, 2) == Block.(1:4)
        @test blockrowsupport(A, 3) == Block.(2:4)
        @test blockrowsupport(A, Block.(3:4)) == Block.(2:4)

        @test blockcolsupport(A, 1) == Block.(1:2)
        @test blockcolsupport(A, 2) == Block.(1:3)
        @test blockcolsupport(A, 3) == Block.(1:4)
        @test blockcolsupport(A, 4) == Block.(2:4)
        @test blockcolsupport(A, Block.(1:2)) == Block.(1:3)
    end

    @testset "block-banded matrix interface for blockranges" begin
        l , u = 1,1
        N = M = 4
        cols = rows = 1:N
        A = BlockBandedMatrix{Int}(undef, rows,cols, (l,u))
        A.data .= 1:length(A.data)

        V = view(A, Block.(2:3), Block.(3:4))
        @test isblockbanded(V)

        B = BlockBandedMatrix(V)
        @test B isa BlockBandedMatrix
        @test blockbandwidths(V) == blockbandwidths(B) == (2,0)
        @test B == V == A[Block.(2:3), Block.(3:4)]

        @test A[Block.(2:3), Block.(3:4)] isa BlockBandedMatrix

        x = randn(size(B,2))
        y = similar(x, size(B,1))
        @test all((similar(y) .= MulAdd(B, x)) .=== (similar(y) .= MulAdd(V,x)))
    end

    @testset "BlockBandedMatrix indexing" begin
        l , u = 1,1
        N = M = 10
        cols = rows = 1:N
        A = BlockBandedMatrix{Float64}(undef, rows,cols, (l,u))
        A.data .= 1:length(A.data)

        A[1,1] = 5
        @test A[1,1] == 5

        @test_throws BandError A[1,4] = 5
        A[1,4] = 0
        @test A[1,4] == 0

        @test A[1:10,1:10] ≈ Matrix(A)[1:10,1:10]

        l , u = 2,1
        N = M = 5
        cols = rows = 1:N

        A = BlockBandedMatrix{Int}(undef, rows,cols, (l,u))
        A.data .= 1:length(A.data)

        @test A[1,2] == 7
        A[1,2] = -5
        @test A[1,2] == -5
    end

    @testset "BlockBandedMatrix BLAS arithmetic" begin
        l , u = 1,1
        N = M = 10
        cols = rows = fill(100,N)
        A = BlockBandedMatrix{Float64}(undef, rows,cols, (l,u))
            A.data .= 1:length(A.data)

        V = view(A, Block(N,N))
        @test MemoryLayout(typeof(V)) == ColumnMajor()

        Y = zeros(cols[N], cols[N])
        @time BLAS.axpy!(2.0, V, Y)
        @test Y ≈ 2A[Block(N,N)]

        Y = BandedMatrix(Zeros(cols[N], cols[N]), (0, 0))
        @test_throws BandError BLAS.axpy!(2.0, V, Y)

        AN = A[Block(N,N)]
        @time BLAS.axpy!(2.0, V, V)
        @test A[Block(N,N)] ≈ 3AN

        A = BlockBandedMatrix(Ones{Float64}((4,6)), [2,2], [2,2,2], (0,2))
        B = BlockBandedMatrix(Ones{Float64}((6,6)), [2,2,2], [2,2,2], (0,1))
        @test sum(A) == 20
        @test sum(B) == 20
        C = BlockBandedMatrix{Float64}(undef, [2,2], [2,2,2], (0,3))
        @test all(mul!(C,A,B) .=== materialize!(MulAdd(1.0,A,B,0.0,similar(C))) .===
                    A*B)
        AB = A*B
        @test AB isa BlockBandedMatrix
        @test Matrix(AB) ≈ Matrix(A)*Matrix(B)
    end

    @testset "BlockBandedMatrix fill and copy" begin
        l , u = 1,1
        N = M = 10
        cols = rows = 1:N
        A = BlockBandedMatrix{Float64}(undef, rows,cols, (l,u))
            A.data .= randn(length(A.data))

        fill!(view(A, Block(2,1)), 2.0)
        @test A[Block(2,1)] ≈ fill(2.0, 2)

        @test_throws BandError fill!(view(A, Block(3,1)), 2.0)

        fill!(view(A, Block(3,1)), 0.0)
        @test A[Block(3,1)] ≈ zeros(3)

        @test_throws BandError fill!(A, 2.0)
        fill!(A, 0.0)
        @test Matrix(A) == zeros(size(A))
    end

    @testset "BlockBandedMatrix type inferrence bug (#9)" begin
        s = BlockBandedMatrices.BlockBandedSizes([1, 2], [1, 2], 0, 0)
        f(s) = s.block_starts.data
        @inferred(f(s))
    end
end


text:
module SPGBox

include("./SPGBoxResult.jl")
include("./VAux.jl")
include("./pr_gradnorm.jl")
include("./spgbox_main.jl")
export spgbox!, spgbox

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# models that can be fitted currently
const FITTABLE = filter(isstationary, setdiff(subtypes(Variogram), (NuggetEffect,NestedVariogram)))

"""
    VariogramFitAlgo

An algorithm for fitting theoretical variograms.
"""
abstract type VariogramFitAlgo end

"""
    WeightedLeastSquares()
    WeightedLeastSquares(w)

Fit theoretical variogram using weighted least squares with weighting
function `w` (e.g. h -> 1/h). If no weighting function is provided,
bin counts of empirical variogram are normalized and used as weights.
"""
struct WeightedLeastSquares <: VariogramFitAlgo
  weightfun::Union{Function,Nothing}
end

WeightedLeastSquares() = WeightedLeastSquares(nothing)

"""
    fit(V, γ, [algo])
    fit(V, γ, [weigthfun])

Fit theoretical variogram type `V` to empirical variogram `γ`
using algorithm `algo`. Default algorithm is `WeightedLeastSquares`.

Alternatively pass the weighting function `weightfun` directly
to the fitting procedure.
"""
function fit(V::Type{<:Variogram}, γ::EmpiricalVariogram,
             algo::VariogramFitAlgo=WeightedLeastSquares())
  # dispatch appropriate implementation
  vario, err = fit_impl(V, γ, algo)

  vario
end

"""
    fit(Variogram, γ, [algo])
    fit(Variogram, γ, [weightfun])

Fit all subtypes of `Variogram` to empirical variogram `γ` and
return the one with minimum error as defined by the algorithm `algo`.

Alternatively pass the weighting function `weightfun` directly
to the fitting procedure.
"""
function fit(::Type{Variogram}, γ::EmpiricalVariogram,
             algo::VariogramFitAlgo=WeightedLeastSquares())
  # fit each variogram type
  res = [fit_impl(V, γ, algo) for V in FITTABLE]
  γs, es = first.(res), last.(res)

  # return best candidate
  γs[argmin(es)]
end

function fit_impl(V::Type{<:Variogram}, γ::EmpiricalVariogram,
                  algo::WeightedLeastSquares)
  # values and distance
  x, y, n = values(γ)
  d = distance(γ)

  # discard invalid bins
  x = x[n .> 0]
  y = y[n .> 0]
  n = n[n .> 0]

  # evaluate weights
  f = algo.weightfun
  w = f ≠ nothing ? map(f, x) : n / sum(n)

  # objective function
  J(p) = begin
    g = V(range=p[1], sill=p[2], nugget=p[3], distance=d)
    sum(w[i]*(g(x[i]) - y[i])^2 for i in eachindex(x))
  end

  # auxiliary variables
  xmax, ymax = maximum(x), maximum(y)

  # initial guess
  pₒ = [xmax/3, .95*ymax, 1e-6]

  # box constraints
  l  = [0., 0., 0.]
  u  = [xmax, ymax, ymax]

  # solve optimization problem
  sol = Optim.optimize(J, l, u, pₒ)
  err = Optim.minimum(sol)
  p   = Optim.minimizer(sol)

  # optimal variogram
  vario = V(range=p[1], sill=p[2], nugget=p[3], distance=d)

  vario, err
end

# convenient methods with weighting function as third argument
fit(V::Type{<:Variogram}, γ::EmpiricalVariogram, weightfun::Function) =
  fit(V, γ, WeightedLeastSquares(weightfun))

fit(V::Type{Variogram}, γ::EmpiricalVariogram, weightfun::Function) =
  fit(V, γ, WeightedLeastSquares(weightfun))


text:
push!(LOAD_PATH,"../src/")

using FrameConnection
using Documenter

makedocs(
    modules=[FrameConnection],
    sitename="FrameConnection.jl",
    format=Documenter.HTML()
)

deploydocs(;
    repo="github.com/LidkeLab/FrameConnection.jl.git",
    devbranch = "main"
)


text:
#=
struct SPSSExport
    df::DataFrame
    code::String
end

function spssexportbyitemgroup(odm::ODMRoot, soid::AbstractString, moid::AbstractString)

    cld    = findclinicaldata(odm, soid, moid)
    bmd    = buildmetadata(odm, soid, moid)
    cldtab = clinicaldatatable(cld)
    evlist = eventlist(bmd)
    iglist = itemgrouplist(bmd; optional = true)
    items  = itemlist(bmd; optional = true)
    ############################################################################
code   = """GET
FILE='filename.sav'.
DATASET NAME main WINDOW=FRONT.
"""
    uiglist = unique(cldtab[!, :ItemGroupOID])

    leftjoin!(cldtab, items[!, ["OID", "SASFieldName"]]; on = "ItemOID" => "OID")

    for ig in eachrow(iglist)
        if ig[:OID] in uiglist
            igc = itemgroupcontent(bmd, ig.OID; optional = true)
            code *= """
DATASET ACTIVATE main.
DATASET COPY sample WINDOW=HIDDEN.
DATASET ACTIVATE sample.
SELECT IF (ItemGroupOID = '$(ig[:OID])').
EXECUTE.
RENAME VARIABLE (ItemOID = IOID).

SORT CASES BY SubjectKey StudyEventOID IOID.
CASESTOVARS
  /ID=SubjectKey StudyEventOID
  /INDEX=IOID
  /GROUPBY=VARIABLE.

"""
            for i in eachrow(igc)
                code *= """RENAME VARIABLE (Value.$(i.OID) = $(i.SASFieldName)).\n"""
            end
            code *= """VARIABLE LABELS"""
            for i in eachrow(igc)
                code *= """\n $(i.SASFieldName) '$(i.Comment)'"""
            end
            code *= """.\n"""
            for i in eachrow(igc)
                itemdef = findelement(bmd, :ItemDef, i.OID)
                if countelements(itemdef, :CodeListRef) > 0
                    codelist = findelements(itemdef, :CodeListRef)
                    code *= """VALUE LABELS $(i.SASFieldName)"""
                    for cl in codelist
                        cldef = findelement(bmd, :CodeList, attribute(cl, "CodeListOID"))
                        for clel in cldef.el
                            if name(clel) == :CodeListItem
                                code *= """\n '$(attribute(clel, "CodedValue"))' '$(findelement(findelement(clel, :Decode), :TranslatedText).el[1].content)'"""
                            end
                        end
                    end
                    code *= """.\n"""
                end
            end
            code *= """SAVE OUTFILE=
    '$(ig[:OID]).sav'
  /COMPRESSED.
DATASET CLOSE sample. \n"""
        end
    end

    SPSSExport(cldtab, code)
end
=#


text:
# Store the coefficient of pressure in and compare it against the
# module PressureMod
# export PressureData, calcPressureCoeff, writeSurfacePressureCoeff, readSurfacePressureCoeff

function calcPressureCoeff(params, q)

  # Assumption of calcFeeStream is speed of sound, c = 1 and rho_inf = 1

  p_inf = 1/params.gamma # Free stream pressure. Refer Pulliam's notes
  m_inf = 0.5*params.Ma*params.Ma

  press = calcPressure(params, q)

  pressCoeff = (press - p_inf)/m_inf

  return pressCoeff
end

@doc """
### writeSurfacePressureCoeff

Writes the pressure coefficient for multiple boundary edges to file. Every rank
writes for its portion of the geometric edge. the files are stored in a
separate directory called pressCoeff. The function deletes any existing
`pressCoeff` directory and creates a new one every time its called.

**Arguments**

*  `mesh` : Abstract mesh type
*  `sbp`  : Summation-By-Parts Operator
*  `eqn`  : Euler equation object
*  `opts` : Options dictionary
*  `g_edges` : geometric edges over which the coefficient of pressure is to be
               computed

"""

function writeSurfacePressureCoeff(mesh::AbstractDGMesh{Tmsh},
         sbp::AbstractOperator, eqn::EulerData{Tsol}, opts, g_edges, nodal_pressCoeff) where {Tmsh, Tsol}

  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
  if mesh.isDG
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end

  my_rank = MPI.Comm_rank(eqn.comm)  # Get the rank of the process

  if my_rank == 0 # Delete existing directory & Create new directory
    if isdir("./pressCoeff")
      rm("./pressCoeff"; recursive=true)
    end
    mkdir("./pressCoeff")
  end

  MPI.Barrier(eqn.comm)

  for itr = 1:length(g_edges)
    g_edge_number = g_edges[itr] # Extract geometric edge number
    filename = string("pressCoeff_edge", g_edge_number, "_rank", my_rank,".dat")
    f = open(joinpath("./pressCoeff",filename), "w")
    # get the boundary array associated with the geometric edge
    itr2 = 0
    for itr2 = 1:mesh.numBC
      if findfirst(mesh.bndry_geo_nums[itr2],g_edge_number) > 0
        break
      end
    end

    start_index = mesh.bndry_offsets[itr2]
    end_index = mesh.bndry_offsets[itr2+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)
    q2 = zeros(Tsol, mesh.numDofPerNode)

    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      global_facenum = idx_range[i]
      for j = 1:mesh.sbpface.numnodes
        q = sview(eqn.q_bndry, :, j, global_facenum)
        convertToConservative(eqn.params, q, q2)
        pressCoeff = calcPressureCoeff(eqn.params, q2)
        println(f, pressCoeff)
        nodal_pressCoeff[itr][j,i] = pressCoeff
      end  # End for j = 1:mesh.sbpface.numnodes
    end    # End for i = 1:nfaces

    close(f)  # Close filestream

  end   # End for itr = 1:length(edges)

  return nothing
end # End evalSurfacePressureCoeff

function readSurfacePressureCoeff(mesh::AbstractDGMesh{Tmsh},
         sbp::AbstractOperator, eqn::EulerData{Tsol}, opts, g_edges, nodal_pressCoeff) where {Tmsh, Tsol}

  my_rank = MPI.Comm_rank(eqn.comm)  # Get the rank of the process

  for itr = 1:length(g_edges)
    g_edge_number = g_edges[itr] # Extract geometric edge number

    # Read in the pressure coefficient from file
    filename = string("pressCoeff_edge", g_edge_number, "_rank", my_rank,".dat")
    f = open(joinpath("./pressCoeff",filename), "r")
    Cp = readdlm(f)
    close(f)

    # get the boundary array associated with the geometric edge
    itr2 = 0
    for itr2 = 1:mesh.numBC
      if findfirst(mesh.bndry_geo_nums[itr2],g_edge_number) > 0
        break
      end
    end

    start_index = mesh.bndry_offsets[itr2]
    end_index = mesh.bndry_offsets[itr2+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)
    q2 = zeros(Tsol, mesh.numDofPerNode)

    itr3 = 1
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      global_facenum = idx_range[i]
      for j = 1:mesh.sbpface.numnodes
        nodal_pressCoeff[itr][j,i] = Cp[itr3]
        itr3 += 1
      end
    end # fir i = 1:nfaces
  end # End for itr = 1:length(g_edges)

  return nothing
end

# end # End module Pressure


text:
function fibonacci(x::Int)
  if x <= 1
      return x
  else
      return fibonacci(x - 1) + fibonacci(x - 2)
  end
  return 0
end

max_stars_repo_name: iabhimanyu/Algorithms

max_issues_repo_name: iabhimanyu/Algorithms

max_forks_repo_name: iabhimanyu/Algorithms


text:
"""
Implementation of Differential Evolution: DE/**selection**/n/**recombination**

The constructor takes following keyword arguments:

- `populationSize`: The size of the population
- `F`: the differentiation (mutation) scale factor (default: 0.9). It's usually defined in range ``F \\in (0, 1+]``
- `n`: the number of differences used in the perturbation (default: 1)
- `selection`: the selection strategy function (default: [`random`](@ref))
- `recombination`: the recombination functions (default: [`uniformbin(0.5)`](@ref))
- `K`: the recombination scale factor (default: 0.5*(F+1))
"""
@kwdef struct DE{F1,F2} <: AbstractOptimizer
    populationSize::Integer = 50
    F::Real = 0.9
    n::Integer = 1
    K::Real = 0.5 * (F + 1)
    selection::F1 = random
    recombination::F2 = uniformbin(0.5)
end
population_size(method::DE) = method.populationSize
default_options(method::DE) = (abstol = 1e-10,)
summary(m::DE) = "DE/$(m.selection)/$(m.n)/$(m.recombination)"

mutable struct DEState{T,IT} <: AbstractOptimizerState
    N::Int
    fitness::Vector{T}
    fittest::IT
end
value(s::DEState) = minimum(s.fitness)
minimizer(s::DEState) = s.fittest

"""Initialization of ES algorithm state"""
function initial_state(method::DE, options, objfun, population)
    T = typeof(value(objfun))
    individual = first(population)
    N = length(individual)
    fitness = fill(maxintfloat(T), method.populationSize)

    # setup initial state
    return DEState(N, fitness, copy(individual))
end

function update_state!(objfun, constraints, state, population::AbstractVector{IT}, method::DE, itr) where {IT}

    # setup
    Np = method.populationSize
    n = method.n
    F = method.F

    offspring = Array{IT}(undef, Np)

    # select base vectors
    bases = method.selection(state.fitness, Np)

    # select target vectors
    for (i,b) in enumerate(bases)
        # mutation
        base = population[b]
        offspring[i] = copy(base)
        # println("$i => base:", offspring[i])

        targets = randexcl(1:Np, [i], 2 * n)
        offspring[i] = differentiation(offspring[i], @view population[targets]; F=F)
        # println("$i => mutated:", offspring[i], ", targets:", targets)

        # recombination
        offspring[i], _ = method.recombination(offspring[i], base)
        # println("$i => recombined:", offspring[i])
    end

    # Create new generation
    fitidx = 0
    minfit = Inf
    for i in 1:Np
        o = apply!(constraints, offspring[i])
        v = value(objfun, o) + penalty(constraints, o)
        if (v <= state.fitness[i])
            population[i] = o
            state.fitness[i] = v
            if v < minfit
                minfit = v
                fitidx = i
            end
        end
    end

    # set best individual
    if fitidx > 0
        state.fittest = population[fitidx]
    end

    return false
end


text:
#v1.0 command file layout:
#JSON dict with these entries:
#"analog waveform" => Dict{String,Any}(Pair{String,Any}("positioner1",Any[3,"positioner1_001",2,"positioner1_001",3,"positioner1_001"]))
#	(This dict maps constant analog channel names to a list of subsequences (defined in the "wave list" dict entry).  For each wave entry in the list, the number of repetitions
#	of that entry is specified in the previous list entry, so the ordering is [wave1count, wave1id, wave2count, wave2id...]
#	see hardware.jl for currently allowed names
#  "digital pulse"   => Dict{String,Any}(Pair{String,Any}("laser1",Any[3,"laser1_001",2,"laser1_001",3,"laser1_001"]),Pair{String,Any}("camera2",Any[3,"camera1_001",2,"ca…
#	(this dict follows the same conventions as the analog waveform entry)
#	see hardware.jl for currently allowed names
#  "metadata"        => Dict{String,Any}(Pair{String,Any}("frames",20),Pair{String,Any}("sample rate",10000),Pair{String,Any}("bi-direction",false),Pair{String,Any}("expo…
#	(another dict with this  metadata:
	#  "frames per stack"       (integer)
	#  "samples per second"  (integer, samples per second)
	#  "bi-direction" (true or false)
	#  "exposure time in seconds"     (floating point)
	#  "sample num"   (integer)
	#  "stacks"       (integer)
        #  "rig"          (string)
	#TODO: eliminate redundancies between these parameters and those recorded in the .imagine header (once we move to always using command files in the backend)

#  "version"         => "v1.0"
#  "wave list"       => Dict{String,Any}(Pair{String,Any}("stimulus1_001",Any[100,0,6000,1,2650,0]),Pair{String,Any}("laser_bidirection",Any[601,0,5000,1,1000,0,5000,
	#This is a dict of building block sequences that are composed to create the sequences in "analog waveform" and "digital pulse"
	#Each of these sequences is a list of integers.  Every other element of the list specifies a voltage or (if digital) a true/false value.  The intervening elements
	#store the number of repetitions of the value.  This is useful to compress signals that remain constant much of the time.
	#The order is [samp1count, samp1value, samp2count, samp2value...]
	#Currently we don't require that entries in the wave dict have consistent key strings, or even that they get used at all.
	#TODO: when writing a set of ImagineSignals, check which wave entries get used and only write those

#instead of adding a CommandList type that shares a single wave lookup dict between multiple ImagineSignals...
#combine into one dict before writing to file. do this by adding the sequence_lookup entries of each command to a growing dict.  Throw error if
#pointers to dicts, arrays, or array values are not equal
function combine_lookups(coms::Vector{ImagineSignal})
    output = Dict()
    for c in coms
        sl = sequence_lookup(c)
        for k in keys(sl)
            if !haskey(output, k)
                output[k] = sl[k]
            elseif pointer(output[k]) != pointer(sl[k])
                if !all(output[k].==sl[k])
                    error("Multiple commands contain different versions of the same sequence $k")
                end
            end
        end
    end
    return output
end

compress_seqnames(c::ImagineSignal) = compress(sequence_names(c))

cam_num_from_name(nm::String) = parse(Int, last(nm))

function initialize_outdict(rig::String, seq_lookup::Dict, which_cams,
                            nstacks, frames_per_stack, exp_time, samp_rate::TTI,
                            nsamps, exp_trig_mode, isbidi) where{TTI<:HasInverseTimeUnits}
    out_dict = Dict{String,Any}()
    ana_dict = Dict{String,Any}()
    dig_dict = Dict{String,Any}()
    out_dict[VERSION_KEY] = VERSION_STRING
    out_dict[ANALOG_KEY] = ana_dict
    out_dict[DIGITAL_KEY] = dig_dict
    out_dict[COMPONENT_KEY] = seq_lookup
    sampr = convert(Int, ustrip(upreferred(samp_rate)))
    @assert isa(sampr, Integer) #TODO: Check for samprate beyond DAQ's capability
    out_dict[METADATA_KEY] = Dict("samples per second" => sampr,
                                            "sample num" => nsamps,
                                            "rig" => rig,
                                            "generated from" => "ImagineInterface")
    for c in which_cams
        camdict = Dict{String, Any}()
        out_dict[METADATA_KEY][c] = camdict
        i = cam_num_from_name(c)
        #in the case that only one value was supplied, assume equal for both cameras
        #Note that the three lines below wouldn't work more than two cameras
        camdict["frames per stack"] = frames_per_stack[min(i, length(frames_per_stack))]
        camdict["stacks"] = nstacks[min(i, length(nstacks))]
        camdict["exposure time in seconds"] = ustrip(uconvert(Unitful.s, exp_time[min(i, length(exp_time))]))
        camdict["exposure trigger mode"] = exp_trig_mode[min(i, length(nstacks))]
        camdict["bidirectional"] = isbidi[min(i, length(nstacks))]
    end
    out_dict["version"] = "v1.1"
    return out_dict
end

function _write_commands!(out_dict, coms)
    dig_dict = out_dict[DIGITAL_KEY]
    ana_dict = out_dict[ANALOG_KEY]
    for c in coms
        if isdigital(c)
            dig_dict[name(c)] = Dict{String,Any}("daq channel"=>daq_channel(c))
            if isoutput(c) #only write the sequence field for outputs
                dig_dict[name(c)]["sequence"] = compress(sequence_names(c))
            elseif !isempty(sequence_names(c))
                @warn "Writing the $(name(c)) channel as an input even though it has one or more sample sequences"
            end
        else
            ana_dict[name(c)] = Dict{String,Any}("daq channel"=>daq_channel(c))
            if isoutput(c) #only write the sequence field for outputs
                ana_dict[name(c)]["sequence"] = compress(sequence_names(c))
            end
        end
    end
    return out_dict
end

function get_missing_monitors(coms_used)
    output = similar(coms_used, 0)
    for c in coms_used
        if hasmonitor(c)
            if findname(coms_used, monitor_name(c)) != nothing
                push!(output, getmonitor(c))
            end
        end
    end
    return output
end

function check_cam_param_counts(cams, params, description::String)
    if length(cams) > 1  && length(params) == 1
        @warn "Two camera output signals were supplied with only one $description argument.  Applying this argument to both cameras"
    elseif length(cams) == 1 && length(params) == 2
        error("Two $description parameters were supplied but only one camera is in use")
    elseif length(cams) != length(params)
        error("Invalid camera and parameter combination")
    end
end

function check_cam_meta(coms::Vector{ImagineSignal}, nstacks, nframes_per_stack, exp_time, exp_trig_mode)
    cams = getcameras(coms)
    check_cam_param_counts(cams, nstacks, "nstacks")
    check_cam_param_counts(cams, nframes_per_stack, "nframes_per_stack")
    check_cam_param_counts(cams, exp_time, "exp_time")
    allowed_modes = ["External Start"; "External Control"; "Fast External Control"]
    for m in exp_trig_mode
        if !in(m, allowed_modes)
            error("Exposure trigger mode $m is unrecognized.  Only these modes are allowed: $(allowed_modes)")
        end
    end
    check_cam_param_counts(cams, exp_trig_mode, "exp_trig_mode")
    frame_counts = zeros(length(cams))
    for i = 1:length(cams)
        frame_counts[i] = count_pulses(cams[i])
        expected_count = nframes_per_stack[min(i,length(nframes_per_stack))] * nstacks[min(i, length(nstacks))]
        if frame_counts[i] != nframes_per_stack[min(i,length(nframes_per_stack))] * nstacks[min(i, length(nstacks))]
            error("Expected $expected_count frames but only counted $(frame_counts[i]) exposure trigger pulses for $(name(cams[i]))")
        end
    end
end

function write_commands(fname::String, coms::Vector{ImagineSignal}, nstacks::NS,
                        nframes_per_stack::NF, exp_time::EXP;
                        exp_trig_mode = ["External Start";], isbidi=false, vertical_lines=nothing,
                        skip_validation=false) where{NS<:Union{Int, Vector{Int}}, NF<:Union{Int, Vector{Int}}, EXP<:Union{HasTimeUnits, Vector{HasTimeUnits}}}
    @assert splitext(fname)[2] == ".json"
    if isa(exp_trig_mode, AbstractString)
        exp_trig_mode = [exp_trig_mode;]
    end
    isused = map(x-> !isoutput(x) || !isempty(x), coms)
    coms_used = coms[isused]
    if !skip_validation
        print("Validating signals before writing...\n")
        check_cam_meta(coms_used, nstacks, nframes_per_stack, exp_time, exp_trig_mode)
        validate_all(coms_used; check_is_sufficient = true, vertical_lines=vertical_lines)
        print("...finished validating signals\n")
    end
    which_cams = map(name, getcameras(coms))
    rig = rig_name(first(coms_used))
    seq_lookup = combine_lookups(coms_used)
    if rig == "ocpi-2" && findname(coms_used, "all lasers")==nothing
        ref_sig = first(getoutputs(coms_used))
        rt = rigtemplate("ocpi-2"; sample_rate = samprate(ref_sig))
        all_las = getname(rt, "all lasers")
        all_las.sequence_lookup = seq_lookup
        append!(all_las, randstring(12), trues(length(ref_sig)))
        check_laser(all_las)
        push!(coms_used, all_las)
    end
    mons = get_missing_monitors(coms_used)
    out_dict = initialize_outdict(rig, seq_lookup, which_cams, nstacks, nframes_per_stack, exp_time, samprate(coms[1]), length(coms[1]), exp_trig_mode, isbidi)
    if !isempty(mons)
        print("Adding the following required monitors (inputs) to the command file: $(map(name, mons)) \n")
    end
    _write_commands!(out_dict, coms_used)
    print("Writing these commands as requested: $(map(name, coms_used)) \n")
    _write_commands!(out_dict, mons)
    f = open(fname, "w")
    JSON.print(f, out_dict)
    close(f)
end

function JSON.Writer.show_element(io::JSON.Writer.JSONContext, s, x::RepeatedValue)
    JSON.show_element(io, s, x.n)
    JSON.show_element(io, s, x.value)
end

max_forks_repo_name: HolyLab/ImagineInterface


text:
import Distributions
using MatrixCompletion
import Random
using TimerOutputs

const to = TimerOutput()

function log_simulation_result(dist::ExponentialFamily, completed_matrix, truth_matrix, type_tracker; io = Base.stdout)
  predicted = predict(dist, forward_map(dist,
                                        completed_matrix[type_tracker[convert(Symbol, dist)]]))
  truth     = truth_matrix[type_tracker[convert(Symbol, dist)]]
  summary   = provide(Diagnostics{Any}(),
                      reference  = truth,
                      input_data = predicted)
  @printf(io, "\nSummary about %s\n%s\n", string(convert(Symbol, dist)), repeat("-", 80))
  show(io, MIME("text/plain"), summary)
  print(io, "\n")
end

@testset "$(format("ADMM Algorithm: Small Input[Gaussian + Bernoulli]"))" begin
  let
    Random.seed!(65536)
    for i in 1:200
      @printf("small case: rank = %d\n", i)
      io = open("./test_result/gaussian_bernoulli/small/rank"*string(2*i)*".txt", "w")
      truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = i), 400, 200),
                               (FixedRankMatrix(Distributions.Bernoulli(0.5),  rank = i),  400, 200)])
      sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
      input_matrix     = sample_model.draw(truth_matrix)
      @timeit to  "Gaussian + Bernoulli" completed_matrix, type_tracker = complete(A            = input_matrix,
                                                                                maxiter      = 200,
                                                                                ρ            = 0.3,
                                                                                use_autodiff = false,
                                                                                gd_iter      = 3,
                                                                                debug_mode   = false,
                                                                                project_rank = 2 * i,
                                                                                io           = io)
      log_simulation_result(Bernoulli(), completed_matrix, truth_matrix, type_tracker, io = io)
      log_simulation_result(Gaussian(),  completed_matrix, truth_matrix, type_tracker, io = io)
      show(io, MIME("text/plain"), to)
      close(io)
    end
  end
end

@testset "$(format("ADMM Algorithm: Medium Input[Gaussian + Bernoulli]"))" begin
  let
    Random.seed!(65536)
    for i in 1:10:500
      @printf("medium case: rank = %d\n", i)
      io = open("./test_result/gaussian_bernoulli/small/rank"*string(2*i)*".txt", "w")
      truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = i), 2000, 1000),
                               (FixedRankMatrix(Distributions.Bernoulli(0.5),  rank = i), 2000, 1000)])
      sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
      input_matrix     = sample_model.draw(truth_matrix)
      @timeit to  "Gaussian + Bernoulli" completed_matrix, type_tracker = complete(A            = input_matrix,
                                                                                maxiter      = 200,
                                                                                ρ            = 0.3,
                                                                                use_autodiff = false,
                                                                                gd_iter      = 3,
                                                                                debug_mode   = false,
                                                                                project_rank = 2 * i,
                                                                                io           = io)
      log_simulation_result(Bernoulli(), completed_matrix, truth_matrix, type_tracker, io = io)
      log_simulation_result(Gaussian(),  completed_matrix, truth_matrix, type_tracker, io = io)
      show(io, MIME("text/plain"), to)
      close(io)
    end
  end
end


text:
@testset "109.convert-sorted-list-to-binary-search-tree.jl" begin
    @test sorted_list_to_bst(ListNode{Int64}([-10, -3, 0, 5, 9])) == TreeNode{Int64}(0, TreeNode{Int64}(-3, TreeNode{Int64}(-10, nothing, nothing), nothing), TreeNode{Int64}(9, TreeNode{Int64}(5, nothing, nothing), nothing))
end


text:
# This file is a part of SimilaritySearch.jl

using Dates

"""
    abstract type NeighborhoodReduction end

Overrides `Base.reduce(::NeighborhoodReduction, res::KnnResult, index::SearchGraph)` to postprocess `res` using some criteria.
Called from `find_neighborhood`, and returns a new KnnResult struct (perhaps a copy of res) since `push_neighborhood` captures
the reference of its output.
"""
abstract type NeighborhoodReduction end
abstract type LocalSearchAlgorithm end

"""
    abstract type Callback end

Abstract type to trigger callbacks after some number of insertions.
SearchGraph stores the callbacks in `callbacks` (a dictionary that associates symbols and callback objects);
A SearchGraph object controls when callbacks are fired using `callback_logbase` and `callback_starting`

"""
abstract type Callback end

### Basic operations on the index

"""
    @with_kw mutable struct Neighborhood

Determines the size of the neighborhood, \$k\$ is adjusted as a callback, and it is intended to affect previously inserted vertices.
The neighborhood is designed to consider two components \$k=in+out\$, i.e. _in_coming and _out_going edges for each vertex.
- The \$out\$ size is computed as \$minsize + \\log(logbase, n)\$ where \$n\$ is the current number of indexed elements; this is computed searching
for \$out\$  elements in the current index.
- The \$in\$ size is unbounded.
- reduce is intended to postprocess neighbors (after search process, i.e., once out edges are computed); do not change \$k\$ but always must return a copy of the reduced result set.

Note: Set \$logbase=Inf\$ to obtain a fixed number of \$in\$ nodes; and set \$minsize=0\$ to obtain a pure logarithmic growing neighborhood.

"""
@with_kw mutable struct Neighborhood
    ksearch::Int32 = 2
    logbase::Float32 = 2
    minsize::Int32 = 2
    reduce::NeighborhoodReduction = SatNeighborhood()
end

Base.copy(N::Neighborhood; ksearch=N.ksearch, logbase=N.logbase, minsize=N.minsize, reduce=copy(N.reduce)) =
    Neighborhood(; ksearch, logbase, minsize, reduce)

struct NeighborhoodSize <: Callback end

"""
    struct SearchGraph <: AbstractSearchContext

SearchGraph index. It stores a set of points that can be compared through a distance function `dist`.
The performance is determined by the search algorithm `search_algo` and the neighborhood policy.
It supports callbacks to adjust parameters as insertions are made.

- `hints`: Initial points for exploration (empty hints imply using random points)

Note: Parallel insertions should be made through `append!` or `index!` function with `parallel_block > 1`

"""
@with_kw struct SearchGraph{DistType<:SemiMetric, DataType<:AbstractDatabase, SType<:LocalSearchAlgorithm}<:AbstractSearchContext
    dist::DistType = SqL2Distance()
    db::DataType = VectorDatabase()
    links::Vector{Vector{Int32}} = Vector{Int32}[]
    locks::Vector{Threads.SpinLock} = Threads.SpinLock[]
    hints::Vector{Int32} = Int32[]
    search_algo::SType = BeamSearch()
    neighborhood::Neighborhood = Neighborhood()
    verbose::Bool = true
end

@with_kw struct SearchGraphCallbacks
    hints::Union{Nothing,Callback} = DisjointHints()
    neighborhood::Union{Nothing,Callback} = NeighborhoodSize()
    hyperparameters::Union{Nothing,Callback} = OptimizeParameters(kind=ParetoRecall())
    logbase::Float32 = 1.5
    starting::Int32 = 8
end

@inline Base.length(g::SearchGraph) = length(g.locks)
include("visitedvertices.jl")

Base.copy(g::SearchGraph;
        dist=g.dist,
        db=g.db,
        links=g.links,
        locks=g.locks,
        hints=g.hints,
        search_algo=copy(g.search_algo),
        neighborhood=copy(g.neighborhood),
        verbose=true
) = SearchGraph(; dist, db, links, locks, hints, search_algo, neighborhood, verbose)

## search algorithms

"""
    SearchGraphPools(results=GlobalKnnResult, vstates=GlobalVisitedVertices, beams=GlobalBeamKnnResult)

A set of pools to alleviate memory allocations in `SearchGraph` construction and searching. Relevant on multithreading scenarious where distance functions, `evaluate`
can call other metric indexes that can use these shared resources (globally defined).

Each pool is a vector of `Threads.nthreads()` preallocated objects of the required type.
"""
struct SearchGraphPools{VisitedVerticesType}
    beams::Vector{KnnResultShift}
    satnears::Vector{KnnResult}
    vstates::VisitedVerticesType
end

@inline function getvstate(len, pools::SearchGraphPools)
    @inbounds v = pools.vstates[Threads.threadid()]
    _init_vv(v, len)
end

@inline function getbeam(bsize::Integer, pools::SearchGraphPools)
    @inbounds reuse!(pools.beams[Threads.threadid()], bsize)
end

@inline function getsatknnresult(pools::SearchGraphPools)
    reuse!(pools.satnears[Threads.threadid()], 1)
end

getpools(::SearchGraph; beams=GlobalBeamKnnResult, satnears=GlobalSatKnnResult, vstates=GlobalVisitedVertices) = SearchGraphPools(beams, satnears, vstates)

include("beamsearch.jl")
## parameter optimization and neighborhood definitions
include("optbs.jl")
include("neighborhood.jl")
include("hints.jl")

"""
    search(index::SearchGraph, q, res; hints=index.hints, pools=getpools(index))

Solves the specified query `res` for the query object `q`.
"""
function search(index::SearchGraph, q, res::KnnResult; hints=index.hints, pools=getpools(index))
    if length(index) > 0
        search(index.search_algo, index, q, res, hints, pools)
    else
        (res=res, cost=0)
    end
end

"""
    execute_callbacks(index::SearchGraph, n=length(index), m=n+1)

Process all registered callbacks
"""
function execute_callbacks(callbacks::SearchGraphCallbacks, index::SearchGraph, n=length(index), m=n+1; force=false)
    if force || (n >= callbacks.starting && ceil(Int, log(callbacks.logbase, n)) != ceil(Int, log(callbacks.logbase, m)))
        callbacks.hints !== nothing && execute_callback(callbacks.hints, index)
        callbacks.neighborhood !== nothing && execute_callback(callbacks.neighborhood, index)
        callbacks.hyperparameters !== nothing && execute_callback(callbacks.hyperparameters, index)
    end
end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
module AutoryBroadcastMacros
    export @T
    export @no_op

    using ReTest
    using RowSheets
    using GenSymx
    using DataStructures: OrderedDict, OrderedSet

    include("utils.jl")
    include("rowmacros.jl")
    include("functionwrappers.jl")

end