text:
"""
    module Events

Provides the abstract [`Event`](#) type and concrete event types.

Events in [`TrainingPhase`](#) and [`ValidationPhase`](#):

- [`EpochBegin`](#) and [`EpochEnd`](#), called at the beginning and end of each
  epoch.
- [`StepBegin`](#) and [`StepEnd`](#), called at the beginning and end of each
  batch.
- [`LossBegin`](#), called after the forward pass but before the loss calculation.

[`TrainingPhase`](#) only:

- [`BackwardBegin`](#), called after forward pass and loss calculation but before gradient
  calculation.
- [`BackwardEnd`](#), called after gradient calculation but before parameter update.

"""
module Events

"""
    abstract type Event

Abstract type for events that callbacks can hook into
"""
abstract type Event end

"""
    EpochBegin()

[`Event`](#) called at the beginning of an epoch.
"""
struct EpochBegin <: Event end

"""
    EpochEnd()

[`Event`](#) called at the end of an epoch.
"""
struct EpochEnd <: Event end

"""
    StepBegin()

[`Event`](#) called at the beginning of a batch.
"""
struct StepBegin <: Event end

"""
    StepEnd()

[`Event`](#) called at the end of a batch.
"""
struct StepEnd <: Event end

"""
    LossBegin()

[`Event`](#) called between calculating `y_pred` and calculating loss
"""
struct LossBegin <: Event end

"""
    BackwardBegin()

[`Event`](#) called between calculating loss and calculating gradients
"""
struct BackwardBegin <: Event end

"""
    BackwardEnd()

[`Event`](#) called between calculating gradients and updating parameters.
"""
struct BackwardEnd <: Event end

export
    # abstract
    Event,
    # concrete
    EpochBegin, EpochEnd,
    StepBegin, StepEnd,
    LossBegin,
    BackwardBegin, BackwardEnd

end # module

using .Events


text:
abstract type HMCAdaptor end

struct NoAdaptor <: HMCAdaptor end

struct MassMatrixAdaptor <: HMCAdaptor end

@with_kw struct StepSizeAdaptor <: HMCAdaptor
    target_acceptance::Float64 = 0.8
end

@with_kw struct NaiveHMCAdaptor <: HMCAdaptor
    target_acceptance::Float64 = 0.8
end

@with_kw struct StanHMCAdaptor <: HMCAdaptor
    target_acceptance::Float64 = 0.8
    n_adapts::Int64 = 500
end

function AHMCAdaptor(
    adaptor::NoAdaptor,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator
)
    return AdvancedHMC.NoAdaptation()
end

function AHMCAdaptor(
    adaptor::MassMatrixAdaptor,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator
)
    return AdvancedHMC.MassMatrixAdaptor(metric)
end

function AHMCAdaptor(
    adaptor::StepSizeAdaptor,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator
)
    return AdvancedHMC.StepSizeAdaptor(adaptor.target_acceptance, integrator)
end

function AHMCAdaptor(
    adaptor::NaiveHMCAdaptor,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator
)
    mma = AdvancedHMC.MassMatrixAdaptor(metric)
    ssa = AdvancedHMC.StepSizeAdaptor(adaptor.target_acceptance, integrator)
    return AdvancedHMC.NaiveHMCAdaptor(mma, ssa)
end

function AHMCAdaptor(
    adaptor::StanHMCAdaptor,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator,
)
    mma = AdvancedHMC.MassMatrixAdaptor(metric)
    ssa = AdvancedHMC.StepSizeAdaptor(adaptor.target_acceptance, integrator)
    stan_adaptor = AdvancedHMC.StanHMCAdaptor(mma, ssa)

    # Initialize state using Stan defaults
    # See: https://mc-stan.org/docs/2_18/reference-manual/hmc-algorithm-parameters.html
    AdvancedHMC.initialize!(stan_adaptor.state, 75, 50, 25, Int(adaptor.n_adapts))
    return stan_adaptor
end

ext: jl


text:
module CitationNetworks

using LightGraphs
using LightGraphs.SimpleGraphs
using LightGraphs.Traversals
using LightGraphs.Traversals: TraversalState, TraversalAlgorithm, BreadthFirst, traverse_graph!, topological_sort
import LightGraphs.Traversals: newvisitfn!, revisitfn!
using LightGraphs.ShortestPaths
using LightGraphs.ShortestPaths: AbstractGraphAlgorithm, AbstractGraphResult
using MetaGraphs

# change LG to 2.0 when available

include("weights_spc.jl")
include("mainpath_global.jl")
include("mainpath_local.jl")
include("mainpath_gbfp.jl")
include("genetic_knowper.jl")
include("convenience.jl")

export weights_spc
export mainpath
export ForwardLocal, BackwardLocal, StandardGlobal, GBFP, FBMP
export genetic_knowper
export set_mp_prop!

end # module


text:
#!/opt/julia0.6/julia

include("/home/ponet/Documents/PhD/Coding/Julia/DFTools.jl/src/DFTypes.jl")
include("/home/ponet/Documents/PhD/Coding/Julia/DFTools.jl/src/DFFileprc.jl")
include("/home/ponet/Documents/PhD/Coding/Julia/DFTools.jl/src/DFSupCalc.jl")
include("/home/ponet/Documents/PhD/Coding/Julia/DFTools.jl/src/DFTbCalc.jl")
include("/home/ponet/Documents/PhD/Coding/Julia/DFTools.jl/src/DFHamiCalc.jl")
include("/home/ponet/Documents/PhD/Coding/Julia/DFTools.jl/src/DFUtils.jl")
# include("src/DFGraphics.jl")
include("/home/ponet/Documents/PhD/Coding/Julia/DFTools.jl/src/DFVisualize.jl")

using GLVisualize, GLWindow, GeometryTypes, Colors, GLAbstraction,GLVisualize.GLVisualizeShader

x = create_TB_model("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/","/home/ponet/Documents/PhD/GeTe/SOC/GeTe_bands.out",[[PhysAtom(0.0,0.0,-0.0239129,-0.155854) for i=1:4]...,[PhysAtom(0.0,0.0,5.5540692,0.318205) for i=1:4]...],Float64);

visualize_wfc(x.wfcs[1],0.35)

# using Meshing,GeometryTypes,ColorTypes,GLVisualize

# sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
#     sqrt(sum(dot(v,v))) - 1 # sphere
# end;

# # meshing marching_cubes test
# m = marching_cubes(sdf,0.1);
# it_colors = [[RGBA(1.0,0.0,0.0,0.6) for i=1:3];[RGBA(0.0,0.0,1.0,0.6) for i=1:3]];
# m_test = GLNormalVertexcolorMesh(vertices=m.vertices,faces=m.faces,color=[it_colors[(i-1)%6+1] for i = 1:length(m.vertices)]);

# screen = glscreen()
# _view(visualize(m_test),camera=:perspective)
# renderloop(screen)


text:
@doc """
    rOfDelay(τ::Float64,dt::Float64,ord::Int64)
    rOfDelay(τ::Float64,dt::Float64,ord::ORDER)

Calculate the delay resolution for given delay `τ`, time resolution `dt` and order `ord` defined by an integer or `ORDER` type.

## Arguments
- `τ::Float64`: time delay to discretise
- `dt::Float64`: time resolution (length)
- `ord::Int64`/`ord::ORDER`: order of the discretisation
""" rOfDelay
function rOfDelay(τ::Real, dt::Real, ord::Integer)
    Int64((τ + 100eps(τ)) ÷ dt + ord ÷ 2)
  # trunc(Int64,τ/dt+ord/2.+1000*eps())
end
function rOfDelay(τ::Real, ord::Union{FullDiscretization,SemiDiscretization})
    rOfDelay(τ, ord.Δt, methodorder(ord))
end
function nStepOfLength(τ::Real,dt::Real)
    Int64((τ + 100eps(τ)) ÷ dt)
end
@doc """
    subMXrange(i::Int64,d::Int64)

Return the index range of the submatrix in the large discretisation matrix (used for constructing the deterministic matrix)

## Arguments
- `i::Int64`: the index of the state space vector the submatrix is multiplied with
- `d::Int64`: dimension of the state space
""" subMxRange
function subMxRange(i::Int64, d::Int64)
    (1:d, i * d + 1:(i + 1) * d)
end
@doc """
    subMXArray(i::Int64,d::Int64)

Return the index range of the submatrix in the large discretisation matrix as a 1 dimensional array (used for constructing the stochastic matrix)

## Arguments
- `i::Int64`: the index of the state space vector the submatrix is multiplied with
- `d::Int64`: dimension of the state space
""" subMXArray
function subMXArray(i::Int64, d::Int64)
    sMXrng = subMxRange(i, d);
    [(i, j) for i in sMXrng[1],j in sMXrng[2]]
end

@doc """
    subVecRange(d::Int64)

Return the index range of the subvector in the large discretised additive vector (it is always added to the present state)

## Arguments
- `d::Int64`: dimension of the state space
""" subVecRange
function subVecRange(d::Int64)
    (1:nDim, 1:1)
end
@doc """
    subVecArray(d::Int64)

Return the index range of the subvector in the large discretised additive vector as a 1 dimensional array (it is always added to the present state)

## Arguments
- `d::Int64`: dimension of the state space
""" subVecArray
function subVecArray(d::Int64)
    sMXrng = subVecRange(d);
    [(i, j) for i in sMXrng[1],j in sMXrng[2]]
end

@doc """
    addSubmatrixToResult!(result::AbstractArray,subMX::detSubMX)

Add the deterministic submatrix `subMX` to the deterministic mapping matrix
Return the index range of the subvector in the large discretised additive vector as a 1 dimensional array (it is always added to the present state)

## Arguments
- `result::AbstractArray`: mapping matrix to add the subMX to
- `subMX::detSubMX`: submatrix to add to the mapping matrix
""" addSubmatrixToResult!
function addSubmatrixToResult!(result::AbstractArray, subMX::SubMX)
    for i in eachindex(subMX.ranges, subMX.MXs)
        result[subMX.ranges[i]...] .+= subMX.MXs[i]
    end
end
function addSubvectorToResults!(results::AbstractArray, subV::SubV)
    results[1:size(subV.V, 1)] .+= subV.V
end
# function addSubvectorToResults!(results::AbstractArray, subV::SubV)
#     d = size(subV.Vs[1], 1)
#     for j in eachindex(subV.Vs, results)
#         for i in eachindex(subV.Vs[j])
#             results[j][1:d] .+= subV.Vs[j][:]
#         end
#     end
# end

# Towards higher order ############################################################
# Building the Lagrangian polynomial for higher order methods
function lagr_atom(i::Real, k::Real, l::Real, τ::Real, r::Real, h::Real, t::Real)
    ((t - τ) - (i - r) * h - l * h) / ((k - l) * h)
end

function lagr_el(q::Real, i::Real, k::Real, τ::Real, r::Real, h::Real, t::Real)
    range = [0:k - 1; k + 1:q]
    prod([lagr_atom(i, k, l, τ, r, h, t) for l in range])
end

function lagr_atom0(k::Real, l::Real, τerr::Real, h::Real, t::Real)
    (t - l * h - τerr) / ((k - l) * h)
  #((t-τ)-(l-r)*h)/((k-l)*h)
end

function lagr_el0(q::Real, k::Real, τerr::Real, h::Real, t::Real)
    range = [0:(k - 1); (k + 1):q]
    prod([lagr_atom0(k, l, τerr, h, t) for l in range])
end

# Handling multiple time-steps
# Reduction of the mapping matrices and -vectors on a time interval
function prodl(mappingMXs::Vector{TM}) where TM<:AbstractMatrix
    mx0 = deepcopy(mappingMXs[1])
    for i in 2:length(mappingMXs)
        mx0 .= mappingMXs[i] * mx0
    end
    return mx0
end
function prodl(mappingMXs::Vector{TM},idxs::AbstractVector{<:Integer}) where TM<:AbstractMatrix
    mx0 = deepcopy(mappingMXs[1][idxs,idxs])
    for i in 2:length(mappingMXs)
        mx0 .= mappingMXs[i][idxs,idxs] * mx0
    end
    return mx0
end
# function prodl(mappingMXs::Vector{TM}) where TM<:AbstractMatrix
#     prod(reverse(mappingMXs))
# end
function reduce_additive(mappingMXs::Vector{TM}, mappingVs::Vector{TV}) where {TM<:AbstractMatrix,TV<:AbstractVector}
    mappingV = deepcopy(mappingVs[1]);
    for i in 2:length(mappingVs)
        mappingV .= mappingMXs[i] * mappingV .+ mappingVs[i]
    end
    return mappingV
end
function reduce_additive(mappingMXs::Vector{TM}, mappingVs::Vector{TV}, idxs::AbstractVector{<:Integer}) where {TM<:AbstractMatrix,TV<:AbstractVector}
    mappingV = mappingVs[1][idxs];
    for i in 2:length(mappingVs)
        mappingV .= mappingMXs[i][idxs,idxs] * mappingV .+ mappingVs[i][idxs]
    end
    return mappingV
end

function reduce_additive(mappingMX0::TM, mappingV0::TV, q::Int) where {TM<:AbstractMatrix,TV<:AbstractVector}
    mappingV = deepcopy(mappingV0);
    for i in 2:q
        mappingV .= mappingMX0 * mappingV .+ mappingV0
    end
    return mappingV
end


text:
## Topology functions
"""
```
findConnectedNode(network::DislocationNetwork, node, connection)
```
Use a `node`'s `connection` to find the other node in the link. Returns a tuple with the link corresponding to the `connection` of `node`, the row of links the connected node appears, and the connected node.

## Examples
We want to find the node connected to `node` 5 via `connection` 3.

```julia
julia> link, rowColLink, connectedNode = findConnectedNode(network, 5, 10)
julia> link
10
julia> oppRowLink
1
julia> connectedNode
25
```
Means that the link corresponding to `connection` 3 of `node` 5 is `network.links[:, 10]`; the node connected to `node` 5 is found on `network.links[oppRowLink, 10]` and is node 25, so the node we're looking for is `network.links[oppRowLink, 10] == 25` and the node whose connection we are looking for is `network.links[3-oppRowLink, 10] == 5`. In this case `oppRowLink == 1`, so through its 3rd connection, node 5 is the second node on link 10, where it is connected to node 25.
"""
function findConnectedNode(network::DislocationNetwork, node, connection)
    links = network.links
    connectivity = network.connectivity

    link = connectivity[2 * connection, node] # Find link that corresponds to node's connection.
    colLink = connectivity[2 * connection .+ 1, node] # Find whether node corresponds to the first or second entry of the link corresponding to connection.
    rowColLink = 3 .- colLink # The other node corresponds to the other entry in the link.
    connectedNode = links[rowColLink, link] # Connected node.

    return link, rowColLink, connectedNode
end
include("DislocationTopologyRemove.jl")
include("DislocationTopologyAdd.jl")
include("DislocationTopologySpecial.jl")
include("DislocationTopologyCollision.jl")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
using Test

@testset "Maximum spanning tree" begin
  @testset "Basic" begin
    graph = complete_graph(5)
    rewards = Dict(Edge(1, 2) => 121.0, Edge(1, 3) => 10.0, Edge(1, 4) => 10.0, Edge(1, 5) => 10.0, Edge(2, 3) => 121.0, Edge(2, 4) => 10.0, Edge(2, 5) => 10.0, Edge(3, 4) => 121.0, Edge(3, 5) => 10.0, Edge(4, 5) => 121.0)

    i = SpanningTreeInstance(graph, rewards)
    s = st_prim(i)
    @test s.instance == i
    @test length(s.tree) == 4
    @test Edge(1, 2) in s.tree
    @test Edge(2, 3) in s.tree
    @test Edge(3, 4) in s.tree
    @test Edge(4, 5) in s.tree
  end

  @testset "Basic" begin
    graph = complete_graph(5)
    rewards = Dict(Edge(2, 5) => 0.0, Edge(3, 5) => 0.0, Edge(4, 5) => 1.0, Edge(1, 2) => 0.0, Edge(2, 3) => 0.0, Edge(1, 4) => 0.0, Edge(2, 4) => 0.0, Edge(1, 5) => 0.0, Edge(1, 3) => 0.0, Edge(3, 4) => 0.0)

    i = SpanningTreeInstance(graph, rewards)
    s = st_prim(i)
    @test s.instance == i
    @test length(s.tree) == 4
    @test length(unique(s.tree)) == 4 # Only unique edges.
    @test Edge(4, 5) in s.tree # Only edge with nonzero cost.
  end
end

@testset "Budgeted maximum spanning tree" begin
  @testset "Basic" begin
    graph = complete_graph(3)
    costs = Dict(Edge(1, 2) => 1.0, Edge(1, 3) => 0.5, Edge(2, 3) => 3.0)
    weights = Dict(Edge(1, 2) => 0, Edge(1, 3) => 2, Edge(2, 3) => 0)

    ε = 0.0001
    budget = 1

    # Lagrangian relaxation.
    i = BudgetedSpanningTreeInstance(graph, costs, weights, budget)
    lagrangian = st_prim_budgeted_lagrangian_search(i, ε)
    @test lagrangian.λ ≈ 0.25 atol=ε
    @test lagrangian.value ≈ 3.75 atol=ε
    @test length(lagrangian.tree) == 2
    @test Edge(1, 3) in lagrangian.tree
    @test Edge(2, 3) in lagrangian.tree
    @test _budgeted_spanning_tree_compute_weight(i, lagrangian.tree) >= budget

    # Helpers.
    a = [Edge(1, 2), Edge(1, 3)]
    b = [Edge(1, 3), Edge(1, 4)]
    @test _solution_symmetric_difference_size(a, b) == 2
    res_a, res_b = _solution_symmetric_difference(a, b)
    @test res_a == [Edge(1, 2)] # Elements that are in a but not in b
    @test res_b == [Edge(1, 4)] # Elements that are in b but not in a

    # Complete approximation algorithm.
    sol = st_prim_budgeted_lagrangian_refinement(i)
    @test sol.instance == i
    s = sol.tree
    @test length(s) == 2
    @test Edge(1, 3) in s # Only important edge in this instance: the only one to have a non-zero weight.
    @test _budgeted_spanning_tree_compute_weight(i, s) >= budget
  end

  @testset "Conformity" begin
    # More advanced tests to ensure the algorithm works as expected.

    # 1.
    graph = complete_graph(2)
    r = Dict(Edge(1, 2) => 0.0)
    w = Dict(Edge(1, 2) => 5)
    i = BudgetedSpanningTreeInstance(graph, r, w, 0)
    s = st_prim_budgeted_lagrangian_refinement(i)
    @assert s.tree == [Edge(1, 2)]

    i = BudgetedSpanningTreeInstance(graph, r, w, 20)
    s = st_prim_budgeted_lagrangian_refinement(i)
    @assert s.tree == Edge{Int}[]
  end
end


text:
# test tools for matrix smoothing

T = Float32
@testset "SNR" begin
A = rand(100,100) .+ 10. .* sin.(0.:2π/100:2π) .- 0.5
A .= T.(A)
# test high SNR
@test all(snr(A) .> 10)
# test type stability
@test eltype(snr(A)) == eltype(A)

C = CorrData()
C.corr = A
@test all(snr(C) .> 10)
end

@testset "smoothing" begin
A = rand(Float32,1000,10)
presmooth = mean(abs.(diff(A,dims=1)),dims=1)
el = eltype(A)
# test return type
@test smooth!(A) == nothing
# test type stability
@test eltype(A) == el
# test smoothness
@test all(presmooth .> mean(abs.(diff(A,dims=1)),dims=1))

@test eltype(smooth(A)) == eltype(A)
@test size(smooth(A)) == size(A)
# test smooth with window lendth of 3
B = smooth(A)
@test isapprox(mean(A[3:9,1]),B[6,1])
end

@testset "abs_max" begin
A = (rand(T,1000,10) .- T(0.5)) .* 1000
@test abs_max!(A) == nothing
@test maximum(A) <= 1.
@test minimum(A) >= -1.
A[:,1] .= 0
@test_throws DomainError abs_max!(A)

# test allocating
A = (rand(T,1000,10) .- T(0.5)) .* 1000
@test size(abs_max(A)) == size(A)
@test eltype(abs_max(A)) == eltype(A)
B = abs_max(A)
@test maximum(B) <= 1.
@test minimum(B) >= -1.

# test CorrData
C = CorrData()
A = (rand(T,1000,10) .- T(0.5)) .* 1000
C.corr = A
@test abs_max!(C) == nothing
@test maximum(C.corr) <= 1.
@test minimum(C.corr) >= -1

# test CorrData
C = CorrData()
A = (rand(T,1000,10) .- T(0.5)) .* 1000
C.corr = A
@test isa(abs_max(C),CorrData)
c = abs_max(C)
@test maximum(c.corr) <= 1.
@test minimum(c.corr) >= -1
end

@testset "standardize" begin
A = (rand(T,1000,10) .- T(0.5)) .* 1000
@test standardize!(A) == nothing
@test all(std(A,dims=1) .≈ T(1.))

# test allocating
A = (rand(T,1000,10) .- T(0.5)) .* 1000
@test size(standardize(A)) == size(A)
@test eltype(standardize(A)) == eltype(A)
B = standardize(A)
@test all(std(B,dims=1) .≈ T(1.))

# test CorrData
C = CorrData()
A = (rand(T,1000,10) .- T(0.5)) .* 1000
C.corr = A
@test standardize!(C) == nothing
@test all(std(C.corr,dims=1) .≈ T(1.))

# test CorrData
C = CorrData()
A = (rand(T,1000,10) .- T(0.5)) .* 1000
C.corr = A
@test isa(standardize(C),CorrData)
c = standardize(C)
@test all(std(c.corr,dims=1) .≈ T(1.))
end

@testset "mad" begin
A = rand(100,10)
A .= T.(A)
# test high mad
@test all(mad(A) .< 1)
# test type stability
@test eltype(mad(A)) == eltype(A)
@test size(mad(A),2) .== size(A,2)

C = CorrData()
C.corr = A
@test all(mad(C) .< 1)
end

@testset "std_threshold" begin
A = rand(T,1000,10)
A[100,1] = T(1000)
@test std_threshold(A,5.) == collect(2:10)
@test std_threshold(A,5) == collect(2:10)
end


text:
module LSRK4
"""
 Low storage Runge-Kutta
 M.H. Carpenter and C.A. Kennedy.
 Fourth-order 2N-storage Runge-Kutta schemes.
 TechnicalReport NASA TM-109112,
 4th order 5-4 solution 3
 y = lsrk4(g,y,t,dt)
"""

# Butcher Tableaus
a = [0
     -567301805773/1357537059087
     -2404267990393/2016746695238
     -3550918686646/2091501179385
     -1275806237668/842570457699]
b = [1432997174477/9575080441755
    5161836677717/13612068292357
    1720146321549/2090206949498
    3134564353537/4481467310338
    2277821191437/14882151754819]
c = [0
    1432997174477/9575080441755
    2526269341429/6820363962896
    2006345519317/3224310063776
    2802321613138/2924317926251]

"""
 Input:
     g: g(y,t) function to integrate in time
     y: Solution vector at current time t
    dy: Temporary storage of the same size as y
     t: Current time
    dt: Time step

 Output:
     y: Solution at time t + dt
"""
function lsrk4(g::Function,
               y::AbstractArray,
               dy::AbstractArray,
               t::Float64, dt::Float64)

        dy = 0 * dy

        for k in 1:size(a, 1)

            # Set rates
            @inbounds dy = a[k] .* dy + g(y, t + c[k] * dt)

            # Update fields
            @inbounds y = y + dt .* b[k] .* dy

        end
        return y
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
using LightGraphs.SimpleGraphs: SimpleEdgeIter

"""
    ColoredEdgeIter

The function [`edges`](@ref) returns a `SimpleEdgeIter` for `AbstractSimpleGraphs`.
The iterates are in lexicographical order, smallest first. The iterator is valid for
one pass over the edges, and is invalidated by changes to the graph.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = path_graph(3);

julia> es = edges(g)
SimpleEdgeIter 2

julia> e_it = iterate(es)
(Edge 1 => 2, (1, 2))

julia> iterate(es, e_it[2])
(Edge 2 => 3, (2, 3))
```
"""
struct ColoredEdgeIter{G} <: AbstractEdgeIter
    g::G
end

eltype(::Type{SimpleEdgeIter{CGT}}) where {CGT<:AbstractColoredGraph} = ColoredEdge{T}

function iterate(eit::ColoredEdgeIter{G}, state=((one(eltype(eit.g)), 1),1)) where {T,G<:AbstractColoredGraph{T}}
    state_uncolored, col_id = state

    iter_res = iterate(edges(uncolored(eit.g)), state_uncolored)
    isnothing(iter_res) && return nothing

    # iterate colors
    edge, state_uncolored_next = iter_res

    col_edge = colored(eit.g, edge, col_id)
    if col_id < length(colors(eit.g, edge))
        return col_edge, (state_uncolored, col_id +1)
    end

    # otherwise iterate edges and reset to color # 1
    return col_edge, (state_uncolored_next, 1)
end

length(eit::ColoredEdgeIter) = ne(eit.g)

function _isequal(e1::ColoredEdgeIter, e2)
    k = 0
    for e in e2
        has_edge(e1.g, e) || return false
        k += 1
    end
    return k == ne(e1.g)
end

in(e, es::ColoredEdgeIter) = has_edge(es.g, e)

show(io::IO, eit::ColoredEdgeIter) = write(io, "ColoredEdgeIter $(ne(eit.g))")


text:
# This code separately visualize the 3D radiation pattern
# of electric and magnetic resonances.
# Because the amplitude of scattering intensity is normalized,
# the setting of wavelength is not very important...

using jlmie
using Plots

# settings
# structure
nmat = 3.5  # refractive index of a sphere
radius = 100*1e-9

# environment
nenv = 1.00

# wavelenth (not very important)
lbd0ED = 618*1e-9  # ED
lbd0MD = 831*1e-9  # MD
lbd0 = [lbd0ED lbd0MD]

# angular range
theta = range(0,π,length=37) # per 5 degrees with length 37
phi = range(0,2π,length=73) # per 5 degrees with length 73

# other settings
n = 1
EorH = 2  # 1: Eelectric, 2: Magnemtic
Isff = zeros(length(theta),length(phi))
for i = 1:length(theta)
    for j = 1:length(phi)
        Isff[i,j] = jlmie_Isff_EM_n(nmat,radius,lbd0[EorH],nenv,theta[i],phi[j],n)[EorH]
    end
end
Isff = Isff/maximum(Isff)
mat_x = Isff.* (sin.(theta) * cos.(phi)')
mat_y = Isff.* (sin.(theta) * sin.(phi)')
mat_z = Isff.* (cos.(theta) * cos.(phi*0)')

# Just a sphere
# mat_x = 1 .* (sin.(theta) * cos.(phi)')
# mat_y = 1 .* (sin.(theta) * sin.(phi)')
# mat_z = 1 .* (cos.(theta) * cos.(phi*0)')

# show results
pyplot()
# show results
plt = plot(mat_x,mat_y,mat_z,
        fill_z=abs.(Isff),
        st=:surface,
        camera=(45,30), # azimuth, elevate
        xlabel="x",
        ylabel="y",
        zlabel="z",
        # colorbar_title="Scattering intensity (arb. units)",
        xlims=(-1,1),
        ylims=(-1,1),
        zlims=(-1,1),
)
display(plt)

# save
if false
    runningfilename = splitext(splitpath(@__FILE__)[end])[1]
    outfilename = ".\\result\\" * runningfilename
    savefig(plt,outfilename)
    println("Graph has been saved in result directory")
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
module Curlcurls

export Curlcurl

using ..Lattices
using LinearAlgebra
using SparseArrays

struct Curlcurl{N}
    cc0::SparseMatrixCSC{ComplexF64,Int}
    cc1::SparseMatrixCSC{ComplexF64,Int}
    cc2::SparseMatrixCSC{ComplexF64,Int}
end

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

Base.conj(cc::Curlcurl{N}) where N = Curlcurl{N}(cc.cc0,-cc.cc1,cc.cc2)

function Base.show(io::IO,cc::Curlcurl{N}) where N
    print(io,"Curlcurl in $(N)D (call with ")
    printstyled(io,"ky",color=:cyan)
    print(io,", ")
    printstyled(io,"kz",color=:cyan)
    print(io,")")
end

end


text:
using PyPlot
using LsqFit
using DelimitedFiles

# Load Config File in /FemtoLabScripts/Data
up_power = readdlm(joinpath(@__DIR__,"../Data/uc_power_20201118.config"), header=true)

# Load Values for max_power, power_value and sample_power
m = up_power[1][:,1]
v = up_power[1][:,2]
s = up_power[1][:,3]

# Define 2D-Model s = p1 + p2*m + p3*m^2 + p4*v + p5*v^2
x = [m v]
@. model(x,p) = p[1] + p[2]*x[:,1] + p[3]*x[:,1]^2 + p[4]*x[:,2] + p[5]*x[:,2]^2

# Set Starting Values
p0 = ones(5)

# Fit the data
fit = curve_fit(model, x, s, p0)

# Generate Data Points for Calculated Values
L = 11
mcalc = range(12, 30, length=L)
vcalc = range(0, 30, length=L)
mvcalc = Array{Float64,2}(undef, 0, 2)
for mm in mcalc, vv in vcalc
    global mvcalc = vcat(mvcalc, [mm vv])
end
scalc = model(mvcalc, fit.param)
scalc[scalc .< 0] .= NaN

# Plot the data
println("Plotting data...")
pygui(true)
figure()
scatter3D(m,v,s, label="real Sample power", color="C1")
mesh(mcalc, vcalc, reshape(scalc, (L,L)))
# scatter3D(model(x,fit.param),m,v, label="calculated Sample Power")
xlabel("Max Power")
ylabel("Power Value")
zlabel("Sample Power")
legend()

println("...done!")

# Print calibration equation
param = round.(fit.param,digits=5)
println("------------------------------------------------------------------")
println("Equation to calculate Sample Power with Max Power and Power Value:")
println("s = $(param[1]) + ($(param[2]))*m + ($(param[3]))*m**2 + ($(param[4]))*v + ($(param[5]))*v**2;")
println("------------------------------------------------------------------")

max_stars_repo_name: PhysChemUniDue/FemtoLabScripts

max_issues_repo_name: PhysChemUniDue/FemtoLabScripts

max_forks_repo_name: PhysChemUniDue/FemtoLabScripts


text:
module LiquidDSP

const rrrf = ("rrrf", Float32,Float32,Float32)
const crcf = ("crcf", Complex64,Float32,Complex64)
const cccf = ("cccf", Complex64,Complex64,Complex64)

include("../deps/deps.jl")
include("util.jl")
include("liquid_types.jl")

include("Filter/filter.jl")

end # module

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
# ---
# title: 850. Rectangle Area II
# id: problem850
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Segment Tree, Line Sweep
# link: <https://leetcode.com/problems/rectangle-area-ii/description/>
# hidden: true
# ---
#
# We are given a list of (axis-aligned) `rectangles`.  Each `rectangle[i] = [x1,
# y1, x2, y2] `, where (x1, y1) are the coordinates of the bottom-left corner,
# and (x2, y2) are the coordinates of the top-right corner of the `i`th
# rectangle.
#
# Find the total area covered by all `rectangles` in the plane.  Since the
# answer may be too large, **return it modulo 10^9 + 7**.
#
# ![](https://s3-lc-
# upload.s3.amazonaws.com/uploads/2018/06/06/rectangle_area_ii_pic.png)
#
# **Example 1:**
#
#
#
#     Input: [[0,0,2,2],[1,0,2,3],[1,0,3,1]]
#     Output: 6
#     Explanation: As illustrated in the picture.
#
#
# **Example 2:**
#
#
#
#     Input: [[0,0,1000000000,1000000000]]
#     Output: 49
#     Explanation: The answer is 10^18 modulo (10^9 + 7), which is (10^9)^2 = (-7)^2 = 49.
#
#
# **Note:**
#
#   * `1 <= rectangles.length <= 200`
#   * `rectanges[i].length = 4`
#   * `0 <= rectangles[i][j] <= 10^9`
#   * The total area covered by all rectangles will never exceed `2^63 - 1` and thus will fit in a 64-bit signed integer.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#Generate animation from filename
function generate_animation_heatmap( fname::String , gifname::String )

    #Get max iteration
    max_iter , tmax = get_t_max(fname)

    #Read Variables
    fid         = h5open( fname , "r")

    R_orbit  = read_attribute(fid, "Rorbit")[1]
    Ω       = read_attribute(fid, "Omega")[1]
    R1      = read_attribute(fid, "R1")[1]
    R2      = read_attribute(fid, "R2")[1]

    #Get initial data
    t , x , y , ψ , _ =    get_fields( fname , 1 )

    #Iteration
    iter = 1:max_iter

    #Print Infomation
    println("---------------------------------------------------")
    println("Starting Frame Sequence...")

    #Nodes to control data
    field   = Node( ψ )
    ttt     = Node( t )
    mtitle  = Node(string("t = " , round( t , digits = 6) , "M" ))

    #Center of regions and index of array
    x_orbit = R_orbit * cos( -Ω * t)
    y_orbit = R_orbit * sin( -Ω * t)

    c1 =  Node( Circle( Point2f0( x_orbit ,  y_orbit) , R1 ) )
    c2 =  Node( Circle( Point2f0( -x_orbit ,  -y_orbit) , R2 ) )

    #Create Figure
    fig = Figure(   resolution = (900, 900) , font = "Helvetica")
    ax = fig[1,1] =  Axis(fig, title = mtitle)

    sc = heatmap!(ax , x , y , field , colorrange=(-5,5) , colormap = :vik ) # color map :vik :roma :seismic :twilight :diverging_bkr_55_10_c35_n256

    poly!( c1 , color = RGBA{Float32}(0.0f0, 0.0f0 , 0.0f0 , 1.0f0)   )
    poly!( c2 , color = RGBA{Float32}(0.0f0, 0.0f0 , 0.0f0 , 1.0f0)   )

    #ProgressBar for animation monitoring
    anim_prog_bar = Progress( max_iter )

    record(fig, gifname , iter; framerate = 24) do iter

        #Update progress bar
        ProgressMeter.update!(anim_prog_bar ,  iter)

        #Update data
        t , _ , _ , ψ , _ =    get_fields( fname , iter)

        x_orbit = R_orbit * cos( Ω * t)
        y_orbit = R_orbit * sin( Ω * t)

        c1[] = Circle( Point2f0( x_orbit ,  y_orbit) , R1 )
        c2[] = Circle( Point2f0( -x_orbit ,  -y_orbit) , R2 )

        field[]  = ψ
        ttt[]     = t
        mtitle[] = string("t = " , round(t , digits = 6) , "M" )

     end

     println("Saving Animation...")

end


text:
module DbgMacro

export @dbg, @dumpct, @dumprt, @qn

"""
    @dbg <expr1> <expr2> <expr3> ...

Display all the expressions in the same way as `@show` does, each on a separate line,
preceded by the location in the format `module:file:line`.

The output goes to `stderr`.

Useful for debugging.

# Examples

```jldoctest
julia> m = [1 2; 3 4];

julia> @dbg 1+2 "Hello" m
Main:none:1  1 + 2 = 3
Main:none:1  "Hello" = "Hello"
Main:none:1  m = [1 2; 3 4]
```
"""
macro dbg(exprs...)
    mod = __module__
    file = relpath(string(__source__.file))
    line = __source__.line
    location = "$mod:$file:$line  "
    esc(quote
        $((:(Base.println(Base.stderr, $location,
            $(ex isa String ? repr(ex) : string(ex)),
            " = ", Base.repr($ex))) for ex in exprs)...)
    end)
end

export @dumpct, @dumprt

"""
    @dumpct <expression>

Dump the expression at compile-time.

# Examples

```jldoctest
julia> function foo(x)
           @dumpct :x + x + \$x
           x
       end
Expr
  head: Symbol call
  args: Array{Any}((4,))
    1: Symbol +
    2: QuoteNode
      value: Symbol x
    3: Symbol x
    4: Expr
      head: Symbol \$
      args: Array{Any}((1,))
        1: Symbol x
foo (generic function with 1 method)

julia> foo(42)
42
```
"""
macro dumpct(ex)
    dump(ex; maxdepth=32)
end

"""
    @dumprt <expression>

Dump the expression at run-time.

# Examples

```jldoctest
julia> function foo(x)
           @dumprt :x + x + \$x
           x
       end
foo (generic function with 1 method)

julia> foo(42)
Expr
  head: Symbol call
  args: Array{Any}((4,))
    1: Symbol +
    2: QuoteNode
      value: Symbol x
    3: Symbol x
    4: Expr
      head: Symbol \$
      args: Array{Any}((1,))
        1: Symbol x
42
```
"""
macro dumprt(ex)
    :(dump($(QuoteNode(ex)); maxdepth=32))
end

"""
    @qn <expression>

Return the expression, without \$ interpolation.

# Examples

```jldoctest
julia> @qn foo(1, x, :y, \$z)
:(foo(1, x, :y, \$(Expr(:\$, :z))))
```
"""
macro qn(ex)
    QuoteNode(ex)
end

end


text:
# this_dir=string(dictGlobalSettings["Julia Code Folder"],"\\src");cd(this_dir);push!(LOAD_PATH,this_dir)
this_dir = "C:\\Users\\bernhard.konig\\Documents\\ASync\\irobot\\TariffWatch New Repo\\algorithms\\Julia\\Code\\dev\\src";cd(this_dir);push!(LOAD_PATH, this_dir)
# include(string(this_dir,"\\DTM.jl"))
using DTM
using DataFrames

#=
reffile="A:\\JuliaTestData\\bagging_small.csv"
reffile="C:\\Users\\bernhard.konig\\Documents\\ASync\\irobot\\data_tmp\\bagging_small.csv"
eltypev=DTM.define_eltypevector(DataFrames.readtable(reffile,nrows=100),DTM.global_const_shift_cols,35)
reference_db=readtable(reffile,eltypes=eltypev,nrows=-1,normalizenames=true) =#

tmpf = "A:\\CH\\03Testing October 2017\\motordata_10k.csv"
tmpf = "C:\\Users\\bernhard.konig\\Documents\\ASync\\irobot\\data_tmp\\motordata_200k.csv"
df_tmp = readtable(tmpf);

expl = ["AFF_MM" "FLG_BERSANI" "f_anz_cli_Rca" "f_citt" "anno" "f_cod_fraz" "f_flg_rinnovo_new" "f_num_gar" "F_rin_cons" "f_num_pol_att" "anno_pat" "fonte" "FLG_CNT_INT" "f_can_pol" "f_profes" "f_flg_cla_tpgu" "f_sex" "f_eta" "f_anz_cli_Zuco" "f_etavei" "f_alim" "F_stat_civ_gran" "f_pagamento_new" "f_att_na_old" "f_att_pulito" "f_cod_marc" "f_ant_rit" "f_reg" "COND_ABI" "UT_AUTO" "KM_AN" "SG_PR" "numpol_att_motor" "Ass" "Rit" "inc_fur" "tut" "Nat" "Soc" "ktot" "Crist" "Inf" "bluk" "eta_conseg_patente" "anni_di_patente"][:]

# prep df
@warn("add option forced_categorical_variables to specify columns which are interpreted as categorical variables even if they are numerical (e.g. [1 2 4 9 1][:] ...)")
@time (df_prepped, this_sett) = prepare_dataframe_for_dtm(df_tmp, numcol="n_claims_artificial", denomcol="exposure", independent_vars=expl);

# remove all rows which have some NA entries
# an alternative is of course to replace the NAs with some specific value
origrows = size(df_prepped, 1)
for i = 1:size(df_prepped, 2)
    df_prepped = df_prepped[.!isna.(df_prepped[i]),:]
end
@show droppedrows = origrows - size(df_prepped, 1)
@show size(df_prepped)

@time updateSettingsMod!(this_sett, model_type="boosted_tree", minWeight=-.15, randomw=.1, learningRate=.1, subsampling_features_prop=1, iterations=50, subsampling_prop=1.00)

outfileStringOnly = "boosti"
datafolder = "C:\\temp\\"
outfilename = string("c:\\temp\\", outfileStringOnly)
dataFilename = "C:\\temp\\where_to_save_jld_file"
key, trn, numeratortrn, denominatortrn, weighttrn, trn_numfeatures, trn_charfeatures_PDA, keyval, val, val_numfeatures, val_charfeatures_PDA, numeratorval, denominatorval, weightval, mappings, sett, num_levels, char_levels, all_levels, all_levels_as_string_vector, names_and_levels, candMatWOMaxValues = prep_data_from_df_modernized(df_prepped, this_sett, datafolder, outfilename, dataFilename, outfileStringOnly)

result_of_runmodel, resulting_model = run_model_actual(key, trn, numeratortrn, denominatortrn, weighttrn, trn_numfeatures, trn_charfeatures_PDA, keyval, val, val_numfeatures, val_charfeatures_PDA, numeratorval, denominatorval, weightval, mappings, datafolder, outfilename, outfileStringOnly, DTM.global_const_shift_cols, sett, num_levels, char_levels, all_levels, all_levels_as_string_vector, names_and_levels, candMatWOMaxValues, dataFilename)

# a


text:
using Printf
using BenchmarkTools
function burgers()
    N = 10001
    U = Vector{Float64}(undef,N)
    u = Vector{Float64}(undef,N)
    x = Vector{Float64}(undef,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    x .= a:dx:b
    dt = 0.5*dx
    tmax = (1-mod(1,dt))/dt
    for i=1:N
        u[i] = sin(x[i])
    end
    U[1] = u[1]
    t = 0*dt
    vel = open("Julia.dat", "w")
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"u\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," \n")
    for i=1:N
        write(vel, @sprintf("%1.9e",x[i])," ",@sprintf("%1.9e",u[i])," \n")
    end
    for j=1:tmax
        for i=2:N-1
            U[i] = u[i]-dt/dx*(0.5*u[i]^2-0.5*(u[i-1]^2))
        end
        u .= U
        if mod(j,100)==0
            t = j*dt
            write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
            write(vel, "variables = \"x\", \"u\""," \n")
            write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," \n")
            for i=1:N
                write(vel, @sprintf("%1.9e",x[i])," ",@sprintf("%1.9e",u[i])," \n")
            end
        end
    end
    t = tmax*dt
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"u\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," \n")
    for i=1:N
        write(vel, @sprintf("%1.9e",x[i])," ",@sprintf("%1.9e",u[i])," \n")
    end
    close(vel)
end

function burgers_Vec()
    N = 10001
    U = Vector{Float64}(undef,N)
    u = Vector{Float64}(undef,N)
    x = Vector{Float64}(undef,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    x .= a:dx:b
    dt = 0.5*dx
    tmax = (1-mod(1,dt))/dt
    u .= sin.(x)
    U[1] = u[1]
    t = 0*dt
    vel = open("Julia.dat", "w")
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"u\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," \n")
    for i=1:N
        write(vel, @sprintf("%1.9e",x[i])," ",@sprintf("%1.9e",u[i])," \n")
    end
    for j=1:tmax
        U[2:N-1] .= u[2:N-1].-dt./dx.*(0.5.*u[2:N-1].^2 .- 0.5.*(u[1:N-2].^2))
        u .= U
        if mod(j,100)==0
            write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
            write(vel, "variables = \"x\", \"u\""," \n")
            write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," \n")
            for i=1:N
                write(vel, @sprintf("%1.9e",x[i])," ",@sprintf("%1.9e",u[i])," \n")
            end
        end
    end
    t = tmax*dt
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"u\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," \n")
    for i=1:N
        write(vel, @sprintf("%1.9e",x[i])," ",@sprintf("%1.9e",u[i])," \n")
    end
    close(vel)
end

function burgers_wof()
    N = 10001
    U = Vector{Float64}(undef,N)
    u = Vector{Float64}(undef,N)
    x = Vector{Float64}(undef,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    x .= a:dx:b
    dt = 0.5*dx
    tmax = (1-mod(1,dt))/dt
    for i=1:N
        u[i] = sin(x[i])
    end
    U[1] = u[1]
    for j=1:tmax
        for i=2:N-1
            U[i] = u[i]-dt/dx*(0.5*u[i]^2-0.5*u[i-1]^2)
        end
        u .= U
    end
end

function burgers_Vec_wof()
    N = 10001
    U = Vector{Float64}(undef,N)
    u = Vector{Float64}(undef,N)
    x = Vector{Float64}(undef,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    x .= a:dx:b
    dt = 0.5*dx
    tmax = (1-mod(1,dt))/dt
    u .= sin.(x)
    U[1] = u[1]
    for j=1:tmax
        U[2:N-1] .= u[2:N-1].-dt./dx.*(0.5.*u[2:N-1].^2 .- 0.5.*(u[1:N-2].^2))
        u .= U
    end
end

results1 = @benchmark burgers()
results2 = @benchmark burgers_Vec()
results3 = @benchmark burgers_wof()
results4 = @benchmark burgers_Vec_wof()


text:
using Photon
using Knet:relu, softmax

# Define a model with fully connected layers.
model = Sequential(
      Dense(256, relu),
      Dense( 10),
      softmax,
)

# Lets use the CPU for this, even if we have a GPU.
setContext(device=:cpu)

# We can call the model to predict a minibatch
x = randn(Float32, 10, 16)
y = model(x)

# Or use the predict function for a single sample.
predict(model, randn(10))

println("### Done ###")


text:
using Optim
using ScenTreesMakie

#############################
# Regularization Ansatz

function myfunc(x)

    trr1.state = x # update states Vector
    #trr1.state[1] = trr2.state[1]
    dist = nestedWasserstein(trr1,trr2,2.0).^2.0

    return dist

end

###################
function myConstraint(x)

    # lower bound such that process is monotonically increasing
    lb = trr1.state[trr1.parent[2:end]]
    tmp = sum( max.([0;lb] .- x,0))

    #upper bound such that fitted tree does not exceed original tree (significantly)
    ub = maximum(trr2.state);
    tmp2 = sum(max.(x .- ub,0) )

    return( tmp + tmp2 )
end

function myRegularizedFunction(x)
    return( myfunc(x) + 1000*myConstraint(x) )

end

trr1 = Tree(404);
trr2 = Tree(404);
n = length(trr1.state);
trr1.state[2:end] = max.(trr1.state[trr1.parent[2:end]] .- trr1.state[2:end],0) + trr1.state[2:end];
x0 = trr1.state;

result = Optim.optimize(myRegularizedFunction, x0, ParticleSwarm(n_particles = 5), Optim.Options(iterations=2000, show_trace=false)  )

#result = Optim.optimize(myRegularizedFunction, x0, NelderMead(), Optim.Options(g_tol = 1e-12, show_trace=true)  )
Optim.minimizer(result)

plot(trr1)
plot(trr2)


text:
function build_julia_copyright_header_buffer(intermediate_dictionary::Dict{String,Any})

    # What is the current year?
    current_year = string(Dates.year(now()))

    # Get comment data from
    buffer = ""
    buffer *= "# ----------------------------------------------------------------------------------- #\n"
    buffer *= "# Copyright (c) $(current_year) Varnerlab\n"
    buffer *= "# Robert Frederick Smith School of Chemical and Biomolecular Engineering\n"
    buffer *= "# Cornell University, Ithaca NY 14850\n"
    buffer *= "#\n"
    buffer *= "# Permission is hereby granted, free of charge, to any person obtaining a copy\n"
    buffer *= "# of this software and associated documentation files (the \"Software\"), to deal\n"
    buffer *= "# in the Software without restriction, including without limitation the rights\n"
    buffer *= "# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n"
    buffer *= "# copies of the Software, and to permit persons to whom the Software is\n"
    buffer *= "# furnished to do so, subject to the following conditions:\n"
    buffer *= "#\n"
    buffer *= "# The above copyright notice and this permission notice shall be included in\n"
    buffer *= "# all copies or substantial portions of the Software.\n"
    buffer *= "#\n"
    buffer *= "# THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n"
    buffer *= "# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n"
    buffer *= "# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n"
    buffer *= "# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n"
    buffer *= "# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n"
    buffer *= "# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n"
    buffer *= "# THE SOFTWARE.\n"
    buffer *= "# ----------------------------------------------------------------------------------- #\n"

    # return -
    return buffer
end


text:
################################################################################
#
#   FACE CENTERED CUBIC (FCC) LATTICE
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcellFCC(
            unitcell_type  :: Type{U},
            implementation :: Val{1}
        ) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[0, 0.5, 0.5],
            Float64[0.5, 0, 0.5],
            Float64[0.5, 0.5, 0]
        ],
        # sites
        S[
            newSite(S, Float64[0,0,0], getDefaultLabelN(LS,1))
        ],
        # bonds
        B[
            newBond(B, 1,1, getDefaultLabel(LB), (+1,0,0)),
            newBond(B, 1,1, getDefaultLabel(LB), (-1,0,0)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,+1,-1)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,-1,+1)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,+1,0)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,-1,0)),
            newBond(B, 1,1, getDefaultLabel(LB), (-1,0,+1)),
            newBond(B, 1,1, getDefaultLabel(LB), (+1,0,-1)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,0,+1)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,0,-1)),
            newBond(B, 1,1, getDefaultLabel(LB), (-1,+1,0)),
            newBond(B, 1,1, getDefaultLabel(LB), (+1,-1,0))
        ]
    )
end

# Implementation
# - implementation 2
# - labels <: Any
function getUnitcellFCC(
            unitcell_type  :: Type{U},
            implementation :: Val{2}
        ) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[0, 0.5, 0.5],
            Float64[0.5, 0, 0.5],
            Float64[0.5, 0.5, 0]
        ],
        # sites
        S[
            newSite(S, Float64[0,0,0], getDefaultLabelN(LS,1))
        ],
        # bonds
        B[
            newBond(B, 1,1, getDefaultLabelN(LB,1), (+1,0,0)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (-1,0,0)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (0,+1,-1)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (0,-1,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (0,+1,0)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (0,-1,0)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (-1,0,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (+1,0,-1)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (0,0,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (0,0,-1)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (-1,+1,0)),
            newBond(B, 1,1, getDefaultLabelN(LB,1), (+1,-1,0)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (+1,+1,-1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (+1,-1,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (-1,+1,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (-1,-1,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (-1,+1,-1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (+1,-1,-1))
        ]
    )
end

# Implementation
# - implementation 4
# - labels <: Any
function getUnitcellFCC(
            unitcell_type  :: Type{U},
            implementation :: Val{4}
        ) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[0, 0.5, 0.5],
            Float64[0.5, 0, 0.5],
            Float64[0.5, 0.5, 0]
        ],
        # sites
        S[
            newSite(S, Float64[0,0,0], getDefaultLabelN(LS,1))
        ],
        # bonds
        B[
            newBond(B, 1,1, getDefaultLabelX(LB), (+1,0,0)),
            newBond(B, 1,1, getDefaultLabelX(LB), (-1,0,0)),
            newBond(B, 1,1, getDefaultLabelX(LB), (0,+1,-1)),
            newBond(B, 1,1, getDefaultLabelX(LB), (0,-1,+1)),
            newBond(B, 1,1, getDefaultLabelY(LB), (0,+1,0)),
            newBond(B, 1,1, getDefaultLabelY(LB), (0,-1,0)),
            newBond(B, 1,1, getDefaultLabelY(LB), (-1,0,+1)),
            newBond(B, 1,1, getDefaultLabelY(LB), (+1,0,-1)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (0,0,+1)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (0,0,-1)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (-1,+1,0)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (+1,-1,0))
        ]
    )
end

# Implementation
# - implementation 5
# - labels <: Any
function getUnitcellFCC(
            unitcell_type  :: Type{U},
            implementation :: Val{5}
        ) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[0, 0.5, 0.5],
            Float64[0.5, 0, 0.5],
            Float64[0.5, 0.5, 0]
        ],
        # sites
        S[
            newSite(S, Float64[0,0,0], getDefaultLabelN(LS,1))
        ],
        # bonds
        B[
            newBond(B, 1,1, getDefaultLabelX(LB), (+1,0,0)),
            newBond(B, 1,1, getDefaultLabelX(LB), (-1,0,0)),
            newBond(B, 1,1, getDefaultLabelX(LB), (0,+1,-1)),
            newBond(B, 1,1, getDefaultLabelX(LB), (0,-1,+1)),
            newBond(B, 1,1, getDefaultLabelY(LB), (0,+1,0)),
            newBond(B, 1,1, getDefaultLabelY(LB), (0,-1,0)),
            newBond(B, 1,1, getDefaultLabelY(LB), (-1,0,+1)),
            newBond(B, 1,1, getDefaultLabelY(LB), (+1,0,-1)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (0,0,+1)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (0,0,-1)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (-1,+1,0)),
            newBond(B, 1,1, getDefaultLabelZ(LB), (+1,-1,0)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (+1,+1,-1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (+1,-1,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (-1,+1,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (-1,-1,+1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (-1,+1,-1)),
            newBond(B, 1,1, getDefaultLabelN(LB,2), (+1,-1,-1))
        ]
    )
end


text:
using DelimitedFiles, DataStructures

function get_features(instance_path::String)
  data = load_matpower(instance_path)
  costs_generators = Dict{String, Tuple{Float64,Float64}}() # bus g => (cg, kg)
  loads = Dict{String, Complex}() # bus => Sl= Pl + im*Ql
  shunts = Dict{String, Complex}() # bus n => g - im *b
  bounds_voltage = Dict{String, Tuple{Float64, Float64}}() # bus => (vmin,vmax)
  bounds_realpower = Dict{String, Tuple{Float64, Float64}}() # bus => (Pmin,Pmax)
  bounds_imagpower = Dict{String, Tuple{Float64, Float64}}() # bus => (Qmin,Qmax)
  max_current = Dict{Tuple{String,String}, Float64}() # line=(bus,bus) => imax

  bus_id_line=SortedDict{Int, Int}()
  bus_id_name=SortedDict{Int, String}()

  checkfor(data, 2, "mpc.baseMVA")
  baseMVA = data[2,3]

  ## Building bus load and shunt information
  i_debut, i_fin = find_numarray(1, data)
  bustype = data[i_debut:i_fin, 2]                # weather bus is 1:"PQ" (generators are not accounted for) or 2:"PV"
  checkfor(data, i_debut-1, "mpc.bus")
  for i=i_debut:i_fin
    id = i-i_debut+1
    busname = "BUS_$id"
    bounds_voltage[busname] = (data[i,13], data[i,12])
    loads[busname] = data[i,3] + im*data[i,4]
    if data[i,5] == data[i,6] == 0
      #no shunts
    else
      shunts[busname] = data[i,5] - im*data[i,6]
    end
    bus_id_line[data[i,1]] = id
    bus_id_name[data[i,1]] = busname
  end

  ## Adding bus generator information
  gen2bus = SortedDict{Int, Int}()
  line2busgen = SortedDict{Int, Tuple{Int, Int}}()

  i_debut, i_fin = find_numarray(i_fin+1, data)
  checkfor(data, i_debut-1, "mpc.gen")

  genind, prevgen = 0, 0
  for i=i_debut:i_fin
    gen2bus[i-i_debut+1] = bus_id_line[data[i,1]]
    busname = bus_id_name[data[i,1]]

    S_min = S_max = 0
    if data[i, 1] == prevgen
      println("Several generators at bus $busname, $(data[i,1])")
      genind += 1
      (Pmin, Pmax) = bounds_realpower[busname]
      (Qmin, Qmax) = bounds_imagpower[busname]
      S_min = Pmin + im*Qmin
      S_max = Pmax + im*Qmax
      # S_min = bus[busname]["Gen"].power_min
      # S_max = bus[busname]["Gen"].power_max
    else
      prevgen = data[i, 1]
      genind = 1
    end
    line2busgen[i-i_debut+1] = (data[i, 1], genind)

    if data[i, 8] > 0 #generator is on
      S_min += data[i,10] + im*data[i,5] #Smin = Pmin + i Qmin
      S_max += data[i,9] + im*data[i,4] #Smax = Pmax + i Qmax
    end
    bounds_realpower[busname] = (real(S_min), real(S_max))
    bounds_imagpower[busname] = (imag(S_min), imag(S_max))

  end

  ## building link information
  i_debut, i_fin = find_numarray(i_fin+1, data)
  checkfor(data, i_debut-1, "mpc.branch")
  for i=i_debut:i_fin
    if data[i, 11] == 0
      #@warn("$link $(linkname.orig)⟶$(linkname.dest) breaker out of service !")
    else
      rs, xs, bc = data[i,3:5]
      Smax = data[i,6]
      τ, θ = data[i, 9:10]
      if Smax != 0
        max_current[(bus_id_name[data[i,1]], bus_id_name[data[i,2]])] = Smax/baseMVA
      end
    end
  end

  i_debut, i_fin = find_numarray(i_fin+1, data)
  if data[i_debut-1, 1] == "mpc.areas"
    @warn("Not loading mpc.areas data.")
    i_debut = i_fin+1
    i_debut, i_fin = find_numarray(i_fin+1, data)
  end

  bus_line_id = SortedDict([val=>key for (key,val) in bus_id_line])

  ## Adding generator cost information
  checkfor(data, i_debut-1, "mpc.gencost")
  genind, cur_genind = 0, data[i_debut, 1]
  for i=i_debut:i_fin
    buslineid, _ = line2busgen[i-i_debut+1]
    busid = bus_id_line[buslineid]

    busname = "BUS_$busid"

    cost_degree = data[i, 4]
    cost_coeffs = data[i, 5:(5+cost_degree-1)]
    if cost_degree == 2
      costs_generators[busname] = (cost_coeffs[1], cost_coeffs[2])
    elseif cost_degree == 3
      costs_generators[busname] = (cost_coeffs[2], cost_coeffs[3])
    end
  end
  return costs_generators, loads, shunts, bounds_voltage, bounds_realpower, bounds_imagpower, max_current
end

## Utils functions
function find_numarray(i_start, data)
  i_debut = i_start
  while !isa(data[i_debut, 1], Int)
    i_debut+=1
  end
  i_fin=i_debut
  while !isa(data[i_fin,1], SubString)
    i_fin += 1
  end
  i_debut, i_fin-1
end

checkfor(data, line_ind, name) = (data[line_ind, 1] == name) || error("Expected ", name, " at line ", line_ind, ", got ", data[line_ind,1], " instead.")

function load_matpower(filename)
  instance_name = split(filename, '.')[1]

  touch(instance_name*".temp")
  f = open(filename)
  out = open(instance_name*".temp", "w")

  # removing all ';' at end of lines
  while !eof(f)
    line = readline(f)
    if length(line) > 0 && line[1] != '%' && line[1] != 'f'
      s = split(line, ";")
      println(out, s[1])
    end
  end
  close(f)
  close(out)

  data = DelimitedFiles.readdlm(instance_name*".temp")
  rm(instance_name*".temp")
  return data
end

costs_generators, loads, shunts, bounds_voltage, bounds_realpower, bounds_imagpower, max_current =
get_features("D:\\repo\\data\\data_Matpower\\matpower\\case1888rte.m")
println(bounds_realpower["BUS_802"])
println(bounds_imagpower["BUS_802"])
# println(loads)


text:
module EvaluationCF

using Persa

using Random: shuffle

include("experimental/resampling.jl")
include("experimental/holdout.jl")
include("experimental/kfolds.jl")

include("metrics/accuracy.jl")

end # module


text:
using TimerOutputs
timeroutput = TimerOutput()
using CuArrays, CUDAnative
using JuMP
using ForwardDiff
using acopf

# case="data/case9241pegase"
# case="data/case1354pegase"
case="data/case9"
# case="data/case30"
# case="data/case118"
function jump_model(opfdata, max_iter)
  println("Ipopt through JuMP")
  @timeit timeroutput "model" begin
    opfmodel, Pg, Qg, Va, Vm = acopf.jumpmodel.model(opfdata; max_iter = max_iter)
  end
  @timeit timeroutput "solve" begin
    opfmodel,status = acopf.jumpmodel.solve(opfmodel,opfdata)
  end

  if status==MOI.LOCALLY_SOLVED
    acopf.jumpmodel.outputAll(opfmodel,opfdata, Pg, Qg, Va, Vm)
  end
  return Pg, Qg, Vm, Va
end
function main()
  max_iter=100
  arraytype = Array
  @timeit timeroutput "load data" begin
  opfdata = acopf.jumpmodel.opf_loaddata(case)
  end
  Pg, Qg, Vm, Va = jump_model(opfdata, max_iter)
  Pg0, Qg0, Vm0, Va0 = acopf.jumpmodel.initialPt_IPOPT(opfdata)
  println("Ipopt interface")
  @timeit timeroutput "warmup" begin
  acopf.IpoptTest.test(Pg0, Qg0, Vm0, Va0, timeroutput, opfdata, arraytype; max_iter = max_iter)
  end
  # reset_timer!(timeroutput)
  # reset_timer!()
  @timeit timeroutput "benchmark" begin
  acopf.IpoptTest.test(Pg0, Qg0, Vm0, Va0, timeroutput, opfdata, arraytype; max_iter = max_iter)
  end
  # t1sPg, t1sPg = acopf.benchmark(opfdata, Pg, Qg, Vm, Va, size(Pg,1), size(Pg,1), 100, timeroutput)
  acopf.benchmark(opfdata, Pg0, Qg0, Vm0, Va0, 3, 3, 100, timeroutput, arraytype)
  # println("Objective: ", ForwardDiff.value.(t1sPg))
  # println("Objective gradient: ", ForwardDiff.partials.(t1sPg).values)
  # println("Objective Hessian: ", [i.values for i in ForwardDiff.partials.(ForwardDiff.partials.(t2sPg).values)])
  # println("Constraint Hessian: ", ForwardDiff.partials.(t2srbalconst))
  show(timeroutput)
  print_timer()
  return
end

main()


text:
import ACE: standardevaluator, graphevaluator
using SparseArrays: SparseMatrixCSC, sparse
using LinearAlgebra: mul!

struct RPIBasis{L, T, BOP, NZ, TIN} <: IPBasis
   pibasis::PIBasis{BOP, NZ, TIN}
   A2Bmaps::NTuple{NZ, SparseMatrixCSC{T, Int}}
   Bz0inds::NTuple{NZ, UnitRange{Int}}
end


text:
using Printf
using LinearAlgebra
using LinearAlgebra: givensAlgorithm

struct SymmetricTridiagonalFactorization{T} <: Factorization{T}
    uplo::Char
    factors::Matrix{T}
    τ::Vector{T}
    diagonals::SymTridiagonal
end

Base.size(S::SymmetricTridiagonalFactorization, i::Integer) = size(S.factors, i)

struct EigenQ{T} <: AbstractMatrix{T}
    uplo::Char
    factors::Matrix{T}
    τ::Vector{T}
end

function Base.getproperty(S::SymmetricTridiagonalFactorization, s::Symbol)
    if s == :Q
        return EigenQ(S.uplo, S.factors, S.τ)
    else
        return getfield(S, s)
    end
end

Base.size(Q::EigenQ) = (size(Q.factors, 1), size(Q.factors, 1))
function Base.size(Q::EigenQ, i::Integer)
    if i < 1
        throw(ArgumentError(""))
    elseif i < 3
        return size(Q.factors, 1)
    else
        return 1
    end
end

function LinearAlgebra.lmul!(Q::EigenQ, B::StridedVecOrMat)
    m, n = size(B)
    if size(Q, 1) != m
        throw(DimensionMismatch(""))
    end
    if Q.uplo == 'L'
        for k in length(Q.τ):-1:1
            for j in 1:size(B, 2)
                b = view(B, :, j)
                vb = B[k + 1, j]
                for i in (k + 2):m
                    vb += Q.factors[i, k]'B[i, j]
                end
                τkvb = Q.τ[k]'vb
                B[k + 1, j] -= τkvb
                for i in (k + 2):m
                    B[i, j] -= Q.factors[i, k]*τkvb
                end
            end
        end
    elseif Q.uplo == 'U'
        for k in length(Q.τ):-1:1
            for j in 1:size(B, 2)
                b = view(B, :, j)
                vb = B[k + 1, j]
                for i in (k + 2):m
                    vb += Q.factors[k, i]*B[i, j]
                end
                τkvb = Q.τ[k]'vb
                B[k + 1, j] -= τkvb
                for i in (k + 2):m
                    B[i, j] -= Q.factors[k, i]'*τkvb
                end
            end
        end
    else
        throw(ArgumentError("Q.uplo is neither 'L' or 'U'. This should never happen."))
    end
    return B
end

function LinearAlgebra.rmul!(A::StridedMatrix, Q::EigenQ)
    m, n = size(A)
    if size(Q, 1) != n
        throw(DimensionMismatch(""))
    end
    if Q.uplo == 'L'
        for k in 1:length(Q.τ)
            for i in 1:size(A, 1)
                a = view(A, i, :)
                va = A[i, k + 1]
                for j in (k + 2):n
                    va += A[i, j]*Q.factors[j, k]
                end
                τkva = va*Q.τ[k]'
                A[i, k + 1] -= τkva
                for j in (k + 2):n
                    A[i, j] -= τkva*Q.factors[j, k]'
                end
            end
        end
    elseif Q.uplo == 'U' # FixMe! Should consider reordering loops
        for k in 1:length(Q.τ)
            for i in 1:size(A, 1)
                a = view(A, i, :)
                va = A[i, k + 1]
                for j in (k + 2):n
                    va += A[i, j]*Q.factors[k, j]'
                end
                τkva = va*Q.τ[k]'
                A[i, k + 1] -= τkva
                for j in (k + 2):n
                    A[i, j] -= τkva*Q.factors[k, j]
                end
            end
        end
    else
        throw(ArgumentError("Q.uplo is neither 'L' or 'U'. This should never happen."))
    end
    return A
end

Base.Array(Q::EigenQ) = lmul!(Q, Matrix{eltype(Q)}(I, size(Q, 1), size(Q, 1)))

function _updateVectors!(c, s, j, vectors)
    @inbounds for i = 1:size(vectors, 1)
        v1 = vectors[i, j + 1]
        v2 = vectors[i, j]
        vectors[i, j + 1] = c*v1 + s*v2
        vectors[i, j]     = c*v2 - s*v1
    end
end

function eigvals2x2(d1::Number, d2::Number, e::Number)
    r1 = (d1 + d2)/2
    r2 = hypot(d1 - d2, 2*e)/2
    return r1 + r2, r1 - r2
end
function eigQL2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matrix)
    dj = d[j]
    dj1 = d[j + 1]
    ej = e[j]
    r1 = (dj + dj1)/2
    r2 = hypot(dj - dj1, 2*ej)/2
    λ = r1 + r2
    d[j] = λ
    d[j + 1] = r1 - r2
    e[j] = 0
    c = ej/(dj - λ)
    if isfinite(c) # FixMe! is this the right fix for overflow?
        h = hypot(one(c), c)
        c /= h
        s = inv(h)
    else
        c = one(c)
        s = zero(c)
    end

    _updateVectors!(c, s, j, vectors)

    return c, s
end

function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), debug::Bool=false) where T<:Real
    d = S.dv
    e = S.ev
    n = length(d)
    blockstart = 1
    blockend = n
    iter = 0
    @inbounds begin
        for i = 1:n - 1
            e[i] = abs2(e[i])
        end
        while true
            # Check for zero off diagonal elements
            for be in blockstart + 1:n
                if abs(e[be - 1]) <= tol*sqrt(abs(d[be - 1]))*sqrt(abs(d[be]))
                    blockend = be - 1
                    break
                end
                blockend = n
            end

            # Deflate?
            if blockstart == blockend
                # Yes
                blockstart += 1
            elseif blockstart + 1 == blockend
                debug && println("Yes, but after sreolving 2x2 block")
                d[blockstart], d[blockend] = eigvals2x2(d[blockstart], d[blockend], sqrt(e[blockstart]))
                e[blockstart] = zero(T)
                blockstart += 1
            else
                # Calculate shift
                sqrte = sqrt(e[blockstart])
                μ = (d[blockstart + 1] - d[blockstart])/(2*sqrte)
                r = hypot(μ, one(T))
                μ = d[blockstart] - sqrte/(μ + copysign(r, μ))

                # QL bulk chase
                singleShiftQLPWK!(S, μ, blockstart, blockend)

                debug && @printf("QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %e, rotations: %d\n", blockstart, blockend, e[blockstart], e[blockend-1], μ, iter += blockend - blockstart)
            end
            if blockstart >= n
                break
            end
        end
    end
    sort!(d)
end

function eigQL!(S::SymTridiagonal{T};
                vectors::Matrix = zeros(T, 0, size(S, 1)),
                tol = eps(T),
                debug::Bool=false) where T<:Real
    d = S.dv
    e = S.ev
    n = length(d)

    if size(vectors, 2) != n
        throw(DimensionMismatch("eigenvector matrix must have $(n) columns but had $(size(vectors, 2))"))
    end
    if size(vectors, 1) > n
        throw(DimensionMismatch("eigenvector matrix must have at most $(n) rows but had $(size(vectors, 1))"))
    end

    blockstart = 1
    blockend = n
    @inbounds begin
        while true
            # Check for zero off diagonal elements
            for be = blockstart+1:n
                if abs(e[be-1]) <= tol*sqrt(abs(d[be-1]))*sqrt(abs(d[be]))
                    blockend = be - 1
                    break
                end
                blockend = n
            end
            # Deflate?
            if blockstart == blockend
                # Yes!
                blockstart += 1
            elseif blockstart + 1 == blockend
                debug && println("Yes, but after solving 2x2 block")
                eigQL2x2!(d, e, blockstart, vectors)
                blockstart += 1
            else
                # Calculate shift
                μ = (d[blockstart + 1] - d[blockstart])/(2*e[blockstart])
                r = hypot(μ, one(T))
                μ = d[blockstart] - (e[blockstart]/(μ + copysign(r, μ)))

                # QL bulk chase
                singleShiftQL!(S, μ, blockstart, blockend, vectors)
                debug && @printf("QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %f\n", blockstart, blockend, e[blockstart], e[blockend-1], μ)
            end
            if blockstart >= n
                break
            end
        end
    end
    p = sortperm(d)
    return d[p], vectors[:,p]
end

function eigQR!(S::SymTridiagonal{T};
                vectors::Matrix = zeros(T, 0, size(S, 1)),
                tol = eps(T),
                debug::Bool=false) where T<:Real
    d = S.dv
    e = S.ev
    n = length(d)
    blockstart = 1
    blockend = n
    @inbounds begin
        while true
            # Check for zero off diagonal elements
            for be in (blockstart + 1):n
                if abs(e[be - 1]) <= tol*sqrt(abs(d[be - 1]))*sqrt(abs(d[be]))
                    blockend = be - 1
                    break
                end
                blockend = n
            end
            # Deflate?
            if blockstart == blockend
                # Yes!
                blockstart += 1
            elseif blockstart + 1 == blockend
                debug && println("Yes, but after solving 2x2 block")
                eigQL2x2!(d, e, blockstart, vectors)
                blockstart += 1
            else
                # Calculate shift
                μ = (d[blockend - 1] - d[blockend])/(2*e[blockend - 1])
                r = hypot(μ, one(T))
                μ = d[blockend] - (e[blockend - 1]/(μ + copysign(r, μ)))

                # QR bulk chase
                singleShiftQR!(S, μ, blockstart, blockend, vectors)

                debug && @printf("QR, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, d[blockend]: %f, μ: %f\n", blockstart, blockend, e[blockstart], e[blockend-1], d[blockend], μ)
            end
            if blockstart >= n
                break
            end
        end
    end
    p = sortperm(d)
    return d[p], vectors[:,p]
end

# Own implementation based on Parlett's book
function singleShiftQLPWK!(S::SymTridiagonal,
                           shift::Number,
                           istart::Integer = 1,
                           iend::Integer = length(S.dv))
    d = S.dv
    e = S.ev
    n = length(d)
    γi = d[iend] - shift
    π = abs2(γi)
    ci = one(eltype(S))
    s = zero(eltype(S))
    @inbounds for i = iend-1:-1:istart
        ei = e[i]
        ζ  = π + ei
        if i < iend-1
            e[i+1] = s*ζ
        end
        ci1    = ci
        ci     = π/ζ
        s      = ei/ζ
        di     = d[i]
        γi1    = γi
        γi     = ci*(di - shift) - s*γi1
        d[i+1] = γi1 + di - γi
        π      = ci == 0 ? ei*ci1 : γi*γi/ci
    end
    e[istart] = s*π
    d[istart] = shift + γi
    S
end

# Own implementation based on Parlett's book
function singleShiftQL!(S::SymTridiagonal,
                        shift::Number,
                        istart::Integer = 1,
                        iend::Integer = length(S.dv),
                        vectors = zeros(eltype(S), 0, size(S, 1)))
    d, e   = S.dv, S.ev
    n      = length(d)
    γi     = d[iend] - shift
    π      = γi
    ci, si = one(eltype(S)), zero(eltype(S))
    @inbounds for i = (iend - 1):-1:istart
        ei        = e[i]
        ci1       = ci
        si1       = si
        ci, si, ζ = givensAlgorithm(π, ei)
        if i < iend - 1
            e[i + 1] = si1*ζ
        end
        di       = d[i]
        γi1      = γi
        γi       = ci*ci*(di - shift) - si*si*γi1
        d[i + 1] = γi1 + di - γi
        π        = ci == 0 ? -ei*ci1 : γi/ci

        # update eigen vectors
        _updateVectors!(ci, si, i, vectors)
    end
    e[istart] = si*π
    d[istart] = shift + γi
    S
end

# Own implementation based on Parlett's book
function singleShiftQR!(S::SymTridiagonal,
                        shift::Number,
                        istart::Integer = 1,
                        iend::Integer = length(S.dv),
                        vectors = zeros(eltype(S), 0, size(S, 1)))
    d, e   = S.dv, S.ev
    n      = length(d)
    γi     = d[istart] - shift
    π      = γi
    ci, si = one(eltype(S)), zero(eltype(S))
    for i = (istart + 1):iend
        ei        = e[i - 1]
        ci1       = ci
        si1       = si
        ci, si, ζ = givensAlgorithm(π, ei)
        if i > istart + 1
            e[i - 2] = si1*ζ
        end
        di       = d[i]
        γi1      = γi
        γi       = ci*ci*(di - shift) - si*si*γi1
        d[i - 1] = γi1 + di - γi
        π        = ci == 0 ? -ei*ci1 : γi/ci

        # update eigen vectors
        _updateVectors!(ci, -si, i - 1, vectors)
    end
    e[iend - 1] = si*π
    d[iend]     = shift + γi
    S
end

function zeroshiftQR!(A::Bidiagonal{T}) where T
    d = A.dv
    e = A.ev
    n = length(d)
    oldc = one(T)
    olds = oldc
    c = oldc
    for i = 1:n-1
        c, s, r = givensAlgorithm(d[i]*c,e[i])
        if i > 1
            e[i-1] = olds*r
        end
        oldc, olds, d[i] = givensAlgorithm(oldc*r,d[i+1]*s)
    end
    h = d[n]*c
    e[n-1] = h*olds
    d[n] = h*oldc
    return A
end

symtri!(A::Hermitian) = A.uplo == 'L' ? symtriLower!(A.data) : symtriUpper!(A.data)
symtri!(A::Symmetric{T}) where {T<:Real} = A.uplo == 'L' ? symtriLower!(A.data) : symtriUpper!(A.data)

# Assume that lower triangle stores the relevant part
function symtriLower!(AS::StridedMatrix{T},
                      τ = zeros(T, size(AS, 1) - 1),
                      u = Vector{T}(undef, size(AS, 1))) where T
    n = size(AS, 1)

    # We ignore any non-real components of the diagonal
    @inbounds for i in 1:n
        AS[i,i] = real(AS[i,i])
    end

    @inbounds for k = 1:(n - 2 + !(T<:Real))

        τk = LinearAlgebra.reflector!(view(AS, (k + 1):n, k))
        τ[k] = τk

        for i in (k + 1):n
            u[i] = AS[i, k + 1]
        end
        for j in (k + 2):n
            ASjk = AS[j, k]
            for i in j:n
                u[i] += AS[i, j]*ASjk
            end
        end
        for j in (k + 1):(n - 1)
            tmp = zero(T)
            for i in j+1:n
                tmp += AS[i,j]'AS[i,k]
            end
            u[j] += tmp
        end

        vcAv = u[k + 1]
        for i in (k + 2):n
            vcAv += AS[i, k]'u[i]
        end
        ξτ2 = real(vcAv)*abs2(τk)/2

        u[k + 1] = u[k + 1]*τk - ξτ2
        for i in (k + 2):n
            u[i] = u[i]*τk - ξτ2*AS[i, k]
        end

        AS[k + 1, k + 1] -= 2real(u[k + 1])
        for i in (k + 2):n
            AS[i, k + 1] -= u[i] + AS[i, k]*u[k + 1]'
        end
        for j in (k + 2):n
            ASjk = AS[j, k]
            uj = u[j]
            AS[j,j] -= 2real(uj*ASjk')
            for i = (j + 1):n
                AS[i, j] -= u[i]*ASjk' + AS[i, k]*uj'
            end
        end
    end
    SymmetricTridiagonalFactorization('L', AS, τ, SymTridiagonal(real(diag(AS)), real(diag(AS, -1))))
end

# Assume that upper triangle stores the relevant part
function symtriUpper!(AS::StridedMatrix{T},
                      τ = zeros(T, size(AS, 1) - 1),
                      u = Vector{T}(undef, size(AS, 1))) where T
    n = LinearAlgebra.checksquare(AS)

    # We ignore any non-real components of the diagonal
    @inbounds for i in 1:n
        AS[i,i] = real(AS[i,i])
    end

    @inbounds for k = 1:(n - 2 + !(T<:Real))
        # This is a bit convoluted method to get the conjugated vector but conjugation is required for correctness of arrays of quaternions. Eventually, it should be sufficient to write vec(x') but it currently (July 10, 2018) hits a bug in LinearAlgebra
        τk = LinearAlgebra.reflector!(vec(transpose(view(AS, k, (k + 1):n)')))
        τ[k] = τk'

        for j in (k + 1):n
            tmp = AS[k + 1, j]
            for i in (k + 2):j
                tmp += AS[k, i]*AS[i, j]
            end
            for i in (j + 1):n
                tmp += AS[k, i]*AS[j, i]'
            end
            u[j] = τk'*tmp
        end

        vcAv = u[k + 1]
        for i in (k + 2):n
            vcAv += u[i]*AS[k, i]'
        end
        ξ = real(vcAv*τk)

        for j in (k + 1):n
            ujt       = u[j]
            hjt       = j > (k + 1) ? AS[k, j] : one(ujt)
            ξhjt      = ξ*hjt
            for i in (k + 1):(j - 1)
                hit = i > (k + 1) ? AS[k, i] : one(ujt)
                AS[i, j] -= hit'*ujt + u[i]'*hjt - hit'*ξhjt
            end
            AS[j, j] -= 2*real(hjt'*ujt) - abs2(hjt)*ξ
        end
    end
    SymmetricTridiagonalFactorization('U', AS, τ, SymTridiagonal(real(diag(AS)), real(diag(AS, 1))))
end

_eigvals!(A::SymmetricTridiagonalFactorization;
          tol = eps(real(eltype(A))),
          debug = false) = eigvalsPWK!(A.diagonals, tol = tol, debug = debug)

_eigvals!(A::SymTridiagonal;
          tol = eps(real(eltype(A))),
          debug = false) = eigvalsPWK!(A, tol = tol, debug = debug)

_eigvals!(A::Hermitian;
          tol = eps(real(eltype(A))),
          debug = false) = eigvals!(symtri!(A), tol = tol, debug = debug)

LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization;
                       tol = eps(real(eltype(A))),
                       debug = false) = _eigvals!(A, tol = tol, debug = debug)

LinearAlgebra.eigvals!(A::SymTridiagonal;
                       tol = eps(real(eltype(A))),
                       debug = false) = _eigvals!(A, tol = tol, debug = debug)

LinearAlgebra.eigvals!(A::Hermitian;
                       tol = eps(real(eltype(A))),
                       debug = false) = _eigvals!(A, tol = tol, debug = debug)

_eigen!(A::SymmetricTridiagonalFactorization;
        tol = eps(real(eltype(A))),
        debug = false) =
    LinearAlgebra.Eigen(eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol, debug = debug)...)

_eigen!(A::SymTridiagonal;
        tol = eps(real(eltype(A))),
        debug = false) =
    LinearAlgebra.Eigen(eigQL!(A, vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), tol = tol, debug = debug)...)

_eigen!(A::Hermitian;
        tol = eps(real(eltype(A))),
        debug = false) = _eigen!(symtri!(A), tol = tol, debug = debug)

LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization;
                     tol = eps(real(eltype(A))),
                     debug = false) = _eigen!(A, tol = tol, debug = debug)

LinearAlgebra.eigen!(A::SymTridiagonal;
                     tol = eps(real(eltype(A))),
                     debug = false) = _eigen!(A, tol = tol, debug = debug)

LinearAlgebra.eigen!(A::Hermitian;
                     tol = eps(real(eltype(A))),
                     debug = false) = _eigen!(A, tol = tol, debug = debug)

function eigen2!(A::SymmetricTridiagonalFactorization;
                 tol = eps(real(float(one(eltype(A))))),
                 debug = false)
    V = zeros(eltype(A), 2, size(A, 1))
    V[1] = 1
    V[end] = 1
    eigQL!(A.diagonals, vectors = rmul!(V, A.Q), tol = tol, debug = debug)
end

function eigen2!(A::SymTridiagonal;
                 tol = eps(real(float(one(eltype(A))))),
                 debug = false)
    V = zeros(eltype(A), 2, size(A, 1))
    V[1] = 1
    V[end] = 1
    eigQL!(A, vectors = V, tol = tol, debug = debug)
end

eigen2!(A::Hermitian;
        tol = eps(float(real(one(eltype(A))))),
        debug = false) = eigen2!(symtri!(A), tol = tol, debug = debug)

eigen2(A::SymTridiagonal;
       tol = eps(float(real(one(eltype(A))))),
       debug = false) = eigen2!(copy(A), tol = tol, debug = debug)

eigen2(A::Hermitian     , tol = eps(float(real(one(eltype(A))))), debug = false) = eigen2!(copy(A), tol = tol, debug = debug)

# First method of each type here is identical to the method defined in
# LinearAlgebra but is needed for disambiguation
function LinearAlgebra.eigvals(A::Hermitian{<:Real})
    T = typeof(sqrt(zero(eltype(A))))
    return eigvals!(LinearAlgebra.copy_oftype(A, T))
end
function LinearAlgebra.eigvals(A::Hermitian{<:Complex})
    T = typeof(sqrt(zero(eltype(A))))
    return eigvals!(LinearAlgebra.copy_oftype(A, T))
end
function LinearAlgebra.eigvals(A::Hermitian)
    T = typeof(sqrt(zero(eltype(A))))
    return eigvals!(LinearAlgebra.copy_oftype(A, T))
end
function LinearAlgebra.eigen(A::Hermitian{<:Real})
    T = typeof(sqrt(zero(eltype(A))))
    return eigen!(LinearAlgebra.copy_oftype(A, T))
end
function LinearAlgebra.eigen(A::Hermitian{<:Complex})
    T = typeof(sqrt(zero(eltype(A))))
    return eigen!(LinearAlgebra.copy_oftype(A, T))
end
function LinearAlgebra.eigen(A::Hermitian)
    T = typeof(sqrt(zero(eltype(A))))
    return eigen!(LinearAlgebra.copy_oftype(A, T))
end

# Aux (should go somewhere else at some point)
function LinearAlgebra.givensAlgorithm(f::Real, g::Real)
    h = hypot(f, g)
    return f/h, g/h, h
end


text:
using DataFrames
using Test

@testset "clean names " begin
    df = DataFrame(ok = 2:3, ok2 = 2:3, ok3=2:3)
    rename!(df, :ok => Symbol("ok-2"))

    @test Symbol.(names(cleannames!(df))) == [:ok_2, :ok2, :ok3]

    @test renamedups!([:ok, :ok_1, :ok_1]) == [:ok, :ok_1, :ok_1_1]
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
module LMDB_DBI
    using LMDB
    using Test

    const dbname = "testdb"
    key = 10
    val = "key value is "

    # Create dir
    mkdir(dbname)
    try

        # Procedural style
        env = create()
        try
            open(env, dbname)
            txn = start(env)
            dbi = open(txn)
            put!(txn, dbi, key+1, val*string(key+1))
            put!(txn, dbi, key, val*string(key))
            put!(txn, dbi, key+2, key+2)
            put!(txn, dbi, key+3, [key, key+1, key+2])
            @test isopen(txn)
            commit(txn)
            @test !isopen(txn)
            close(env, dbi)
            @test !isopen(dbi)
        finally
            close(env)
        end
        @test !isopen(env)

        # Block style
        create() do env
            set!(env, LMDB.NOSYNC)
            open(env, dbname)
            start(env) do txn
                open(txn, flags = Cuint(LMDB.REVERSEKEY)) do dbi
                    k = key
                    value = get(txn, dbi, k, String)
                    println("Got value for key $(k): $(value)")
                    @test value == val*string(k)
                    delete!(txn, dbi, k)
                    k += 1
                    value = get(txn, dbi, k, String)
                    println("Got value for key $(k): $(value)")
                    @test value == val*string(k)
                    delete!(txn, dbi, k, value)
                    @test_throws LMDBError get(txn, dbi, k, String)
                    k += 1
                    value = get(txn, dbi, k, Int)
                    println("Got value for key $(k): $(value)")
                    @test value == k
                    k += 1
                    value = get(txn, dbi, k, Vector{Int})
                    println("Got value for key $(k): $(value)")
                    @test value == [key, key+1, key+2]
                end
                abort(txn)
            end
        end
    finally
        rm(dbname, recursive=true)
    end
end


text:
# Base.include(Main, "dataStruct.jl")
# Base.MainInclude.include("../../RawFileofflineProc/dataStruct.jl")
function readin!(cfgdata::cfgStruct,CfgArray::Array{UInt16,1})
    cfgdata.meterperSeg =CfgArray[1] /5000
    cfgdata.chn2SegBegIdx =CfgArray[2]
    cfgdata.reflectorFactor =CfgArray[3] /10000
    cfgdata.attenutionFactor =CfgArray[4] /10000
    cfgdata.scanrate = CfgArray[5]
    cfgdata.calSamples = CfgArray[6]
    cfgdata.SegNumber = CfgArray[7]
end

function show( dataObj::RNG )
    println( "RNG Data :-----------------------------------" )
    println( " AlrmProcStartID:$(dataObj.AlrmProcStartID)")
    println( " AlrmProcEndtID1:$(dataObj.AlrmProcEndtID1)")
    println( " startID1:$(dataObj.startID1)")
    println( " endID1:$(dataObj.endID1)")
    println( " chnid:$(dataObj.chnid)")
    println( " ProtectSpaceWidth:$(dataObj.ProtectSpaceWidth)")
    println("----------------------")
end
###################
#############
function seg2Meters( CfgInfoObj::cfgStruct,segID::Int  )
    meter = CfgInfoObj.meterperSeg
    chn2Pos = CfgInfoObj.chn2SegBegIdx
    ## write the title s: 通道号下的光缆 时间、各段的米标 。
    # print( "时间" ,"," )
    # for segID = CumSumDataobj.AlrmProcStartID : CumSumDataobj.AlrmProcEndtID
    if segID  >= chn2Pos;
        return ( round( UInt,Float64(segID-chn2Pos)* meter) ,2 )
    else
        return (round( UInt,Float64( segID  )* meter ),1)
    end
end

function getRNG!( x::RNG,cfgObj::cfgStruct, startMeters::Int ,endmeters::Int ,ProtectMeters::Int,chnid::Int )
    ProtectSegs =  round(Int, Float64( ProtectMeters ) / cfgObj.meterperSeg)+1
    AlrmProcStartID =  round(Int, Float64( startMeters ) / cfgObj.meterperSeg)
    AlrmProcEndtID1 =   round(Int, Float64( endmeters ) / cfgObj.meterperSeg)

    if chnid == 1 ;
        AlrmProcStartID < 1 ? AlrmProcStartID = 1 : 0  ;
        AlrmProcEndtID1 >=  cfgObj.chn2SegBegIdx  ? AlrmProcEndtID1 = cfgObj.chn2SegBegIdx -1 : 0  ;

        startID1 = AlrmProcStartID - ProtectSegs
        endID1 = AlrmProcEndtID1 + ProtectSegs

        startID1 < 1 ? startID1 = 1 :  0 ;
        endID1 >=  cfgObj.chn2SegBegIdx  ? endID1 = cfgObj.chn2SegBegIdx : 0 ;
    else
        AlrmProcEndtID1 +=  cfgObj.chn2SegBegIdx
        AlrmProcStartID += cfgObj.chn2SegBegIdx

        AlrmProcEndtID1 >  cfgObj.SegNumber  ? AlrmProcEndtID1 = cfgObj.SegNumber :  0 ;
        AlrmProcStartID >  cfgObj.SegNumber  ? AlrmProcStartID = cfgObj.SegNumber :  0 ;

        startID1 = AlrmProcStartID - ProtectSegs
        endID1 = AlrmProcEndtID1 + ProtectSegs

        startID1 < cfgObj.chn2SegBegIdx ? startID1 = cfgObj.chn2SegBegIdx : 0 ;
        endID1 >  cfgObj.SegNumber  ? endID1 = cfgObj.SegNumber :  0 ;
        # println( AlrmProcEndtID1 ,"=======$AlrmProcStartID===$endID1 ===$startID1====")
    end

    x.ProtectSpaceWidth = ProtectSegs

    x.AlrmProcStartID = AlrmProcStartID
    x.AlrmProcEndtID1 = AlrmProcEndtID1
    x.startID1 = startID1
    x.endID1 = endID1
    x.chnid = chnid >= 2 ? chnid =2 : chnid =1;

end

function setValue!( x::cfgStruct  , y::CfgDataStruct)
    x.meterperSeg     = y.SpatialResolution
    x.chn2SegBegIdx   = y.LenCh1 +1
    x.reflectorFactor = y.Refractivity
    x.attenutionFactor=y.ATT
    x.scanrate        =y.ScanRate
    x.calSamples      = y.EnergyCnt
    x.SegNumber       = y.LenCh1 +y.LenCh2
end
function show( dataObj::CfgDataStruct   )
    println( "Config Data :-----------------------------------" )
    println( " SpatialResolution:$(dataObj.SpatialResolution)")
    println( " Refractivity:$(dataObj.Refractivity)")
    println( " OffsetM:$(dataObj.OffsetM)")
    println( " ATT:$(dataObj.ATT)")
    println( " ScanRate:$(dataObj.ScanRate)")
    println( " EnergyCnt:$(dataObj.EnergyCnt)")
    println( " LcLimitHigh:$(dataObj.LcLimitHigh)")
    println( " LcLimitLow:$(dataObj.LcLimitLow)")
    println( " StartCh1:$(dataObj.StartCh1)")
    println( " LenCh1:$(dataObj.LenCh1)")
    println( " StartCh2:$(dataObj.StartCh2)")
    println( " LenCh2:$(dataObj.LenCh2)")
    # println( "SpatialResolution:$(dataObj.SpatialResolution)")

    println("----------------------")
end
function show(cfgdata::cfgStruct)
    println(" meterPerSeg = $(cfgdata.meterperSeg) ")
    println(" chn2SegBegIdx = $(cfgdata.chn2SegBegIdx )")
    println(" reflectorFactor = $(cfgdata.reflectorFactor) ")
    println(" attenutionFactor = $(cfgdata.attenutionFactor) ")
    println(" calSamples = $(cfgdata.calSamples) ")
    println(" scanrate = $(cfgdata.scanrate) ")
    println(" SegNumber = $(cfgdata.SegNumber) ")
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
##

using FreeTypeAbstraction
using AprilTags
using Images, ImageView
using FileIO
using Optim

using ImageDraw, Colors, FixedPointNumbers

##

calibfiles = [
  ENV["HOME"]*"/data/cameracalib/CommonData/IMG_2909.jpg";
  ENV["HOME"]*"/data/cameracalib/CommonData/IMG_2910.jpg";
  ENV["HOME"]*"/data/cameracalib/CommonData/IMG_2911.jpg";
  ENV["HOME"]*"/data/cameracalib/CommonData/IMG_2913.jpg";
  ENV["HOME"]*"/data/cameracalib/CommonData/IMG_2914.jpg";
  ENV["HOME"]*"/data/cameracalib/CommonData/IMG_2916.jpg";
  ENV["HOME"]*"/data/cameracalib/CommonData/IMG_2920.jpg";
]

##

INITTEST = load(calibfiles[1])

f_width  = size(INITTEST,1) + 0.0
f_height = f_width  + 0.0
c_width  = size(INITTEST,2) / 2 + 0.0
c_height = size(INITTEST,1) / 2 + 0.0

s = 0.0

taglength = 0.03

##

arr = Vector{typeof(INITTEST)}()

for i in 1:length(calibfiles)
  push!(arr, load(calibfiles[i]))
end

## quick test

detector = AprilTagDetector()
tags = detector.(arr)

#

ARR = []

for tag in tags
  push!(ARR, [(;id=tg.id,p=tg.p, H=tg.H) for tg in tag])
end

##

resid = calcCalibResidualAprilTags!( arr, ARR, f_height=f_height, taglength=taglength )

##

obj = (f_w, f_h, c_w, c_h) -> calcCalibResidualAprilTags!(arr, ARR,
                                                          # tagList=40:40;
                                                          f_width=f_w,
                                                          f_height=f_h,
                                                          c_width=c_w,
                                                          c_height=c_h,
                                                          taglength=taglength )

obj_ = (fc_wh) -> obj(fc_wh...)

## current best guess

f_width_  = f_width
f_height_ = f_height
c_width_  = c_width
c_height_ = c_height

##

# check that it works
obj_([f_width_, f_height_, c_width_, c_height_])

##

result = optimize(obj_, [f_width_; f_height_; c_width_; c_height_ ], BFGS(), Optim.Options(iterations = 100, show_trace=true, x_tol=1e-8))

##

f_width_  = result.minimizer[1]
f_height_ = result.minimizer[2]
c_width_  = result.minimizer[3]
c_height_ = result.minimizer[4]

##

# iPhone 8 rear camera (coarse calibration)
# f_height_ = 3346.1894
# f_width_ = 3346.1894
# c_height_ = 2021.11068
# c_width_ = 1471.0241

# f_height_ = 3371.2553294118493
# f_width_ = 3353.696574041437
# c_height_ = 2007.7796750349364
# c_width_ = 1496.4523912712611

# f_height_ = 3370.4878918701756
# f_width_  = 3352.8348099534364
# c_height_ = 2005.641610450976
# c_width_  = 1494.8282013012076

# f_height_ = 3431.554669353193
# f_width_  = 3411.8526620805414
# c_height_ = 2016.7496065830883
# c_width_  = 1501.0475003688337

minim = obj(f_width_, f_height_, c_width_, c_height_)

## draw what is going on

SEL = 4

cimg_ = deepcopy(arr[SEL])
calcCornerProjectionsAprilTags!(cimg_, ARR[SEL],
                                # tagList=18:18,
                                taglength=0.0315,
                                f_width=f_width_,
                                f_height=f_height_,
                                c_width=c_width_,
                                c_height=c_height_,
                                dodraw=true )
#

##

ImageView.imshow(cimg_)

##

freeDetector!(detector)

##


text:
import Base: getindex
import Oceananigans.Fields: interior, interiorparent

using Base: @propagate_inbounds

using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑyᵃᶠᵃ, ℑzᵃᵃᶠ
using Oceananigans.Fields

struct LazyPrimitiveField{X, Y, Z, A, G, C, D} <: AbstractField{X, Y, Z, A, G}
          architecture :: A
                  grid :: G
    conservative_field :: C
               density :: D
end

LazyPrimitiveField(LX, LY, LZ, arch, grid, ρϕ, ρ) =
    LazyPrimitiveField{LX, LY, LZ, typeof(arch), typeof(grid), typeof(ρϕ), typeof(ρ)}(arch, grid, ρϕ, ρ)

@inline @propagate_inbounds getindex(f::LazyPrimitiveField{Cell, Cell, Cell}, I...) =
    @inbounds f.conservative_field[I...] / f.density[I...]

@inline @propagate_inbounds getindex(f::LazyPrimitiveField{Face, Cell, Cell}, I...) =
    @inbounds f.conservative_field[I...] / ℑxᶠᵃᵃ(I..., f.grid, f.density)

@inline @propagate_inbounds getindex(f::LazyPrimitiveField{Cell, Face, Cell}, I...) =
    @inbounds f.conservative_field[I...] / ℑyᵃᶠᵃ(I..., f.grid, f.density)

@inline @propagate_inbounds getindex(f::LazyPrimitiveField{Cell, Cell, Face}, I...) =
    @inbounds f.conservative_field[I...] / ℑzᵃᵃᶠ(I..., f.grid, f.density)

LazyVelocityFields(arch, grid, ρ, ρũ) =
    (u = LazyPrimitiveField(Face, Cell, Cell, arch, grid, ρũ.ρu, ρ),
     v = LazyPrimitiveField(Cell, Face, Cell, arch, grid, ρũ.ρv, ρ),
     w = LazyPrimitiveField(Cell, Cell, Face, arch, grid, ρũ.ρw, ρ))

function LazyTracerFields(arch, grid, ρ, ρc̃)
    c_names = [filter(c -> c != 'ρ', string(c)) for c in keys(ρc̃)]
    c_names = filter(s -> s != "", c_names) .|> Symbol |> Tuple  # Don't include the ρ tracer.

    c_fields = Tuple(
        LazyPrimitiveField(Cell, Cell, Cell, arch, grid, getproperty(ρc̃, Symbol(:ρ, c)), ρ)
        for c in c_names
    )

    return NamedTuple{c_names}(c_fields)
end

function interior(f::LazyPrimitiveField)
    data = zeros(size(f.conservative_field))
    for k in 1:f.grid.Nz, j in 1:f.grid.Ny, i in 1:f.grid.Nx
        data[i, j, k] = f[i, j, k]
    end
    return data
end

interiorparent(lf::LazyPrimitiveField) = interior(lf)


text:
using Documenter

deploydocs(repo = "github.com/uva-bi-sdad/DSPG-Manual.git",
           push_preview = true)

max_stars_repo_name: uva-bi-sdad/DSPG-Manual

max_issues_repo_name: uva-bi-sdad/DSPG-Manual

max_forks_repo_name: uva-bi-sdad/DSPG-Manual


text:
import Base: rand, length
import AbstractGPs: logpdf, AbstractMvNormal

export elbo, dtc
export SparseFiniteGP

"""
    SparseFiniteGP{T1<:FiniteGP, T2<:FiniteGP}

A finite-dimensional projection of an `AbstractGP` `f` at locations `x`, which uses a second
`FiniteGP` defined at a sparse set of inducing points [1] to do approximate inference.

This object has similar methods to an ordinary `FiniteGP`, but when you call `logpdf` on it,
it actually computes the `elbo`, and when you call `posterior` on it, `approx_posterior` is
called.

[1] - M. K. Titsias. "Variational learning of inducing variables in sparse Gaussian
processes". In: Proceedings of the Twelfth International Conference on Artificial
Intelligence and Statistics. 2009.

```jldoctest
julia> f = wrap(GP(Matern32Kernel()), GPC());

julia> fobs = f(rand(100));

julia> finducing = f(0:0.1:1);

julia> fxu = SparseFiniteGP(fobs, finducing);

julia> y = rand(fxu);

julia> logpdf(fxu, y) < logpdf(fobs, y)
true
```
"""
struct SparseFiniteGP{T1<:FiniteGP, T2<:FiniteGP} <: AbstractMvNormal
    fobs::T1
    finducing::T2
end

Base.length(f::SparseFiniteGP) = length(f.fobs)

mean(f::SparseFiniteGP) = mean(f.fobs)

const __covariance_error = "The covariance matrix of a sparse GP can often be dense and " *
    "can cause the computer to run out of memory. If you are sure you have enough " *
    "memory, you can use `cov(f.fobs)`."

cov(f::SparseFiniteGP) = error(__covariance_error)

marginals(f::SparseFiniteGP) = marginals(f.fobs)

rand(rng::AbstractRNG, f::SparseFiniteGP, N::Int) = rand(rng, f.fobs, N)
rand(f::SparseFiniteGP, N::Int) = rand(Random.GLOBAL_RNG, f, N)
rand(rng::AbstractRNG, f::SparseFiniteGP) = vec(rand(rng, f, 1))
rand(f::SparseFiniteGP) = vec(rand(f, 1))

elbo(f::SparseFiniteGP, y::AV{<:Real}) = elbo(VFE(f.finducing), f.fobs, y)

logpdf(f::SparseFiniteGP, y::AV{<:Real}) = elbo(VFE(f.finducing), f.fobs, y)

function logpdf(f::SparseFiniteGP, Y::AbstractMatrix{<:Real})
    return map(y -> logpdf(f, y), eachcol(Y))
end

function posterior(f::SparseFiniteGP, y::AbstractVector{<:Real})
    return posterior(AbstractGPs.VFE(f.finducing), f.fobs, y)
end


text:
using DataFrames
using CSV
using PyPlot

df = CSV.read("bullettin.csv");
names(df)

fig, ax = PyPlot.subplots(figsize=(18,15))
ax.plot(df.day, df.susceptible, lw=10, label="Susceptible")
ax.plot(df.day, df.infectus, lw=10, label="Infectus")
ax.plot(df.day, df.recovered, lw=10, label="Recovered")

ax.set_xlabel("Days", fontsize=35)
ax.set_ylabel("Count", fontsize=35)

ax.tick_params(labelsize=35)

ax.legend(fontsize=35)
PyPlot.show()


text:
module style
import Parameters: @with_kw

import Term: unspace_commas,
            NAMED_MODES,
            has_markup,
            OPEN_TAG_REGEX,
            replace_text,
            get_last_ANSI_code,
            CODES,
            ANSICode,
            tview

import ..color: AbstractColor,
            NamedColor,
            is_color,
            is_background,
            get_color,
            is_hex_color,
            hex2rgb

export apply_style

"""
Check if a string is a mode name
"""
is_mode(string) = string ∈ NAMED_MODES

# ---------------------------------------------------------------------------- #
#                                  MarkupStyle                                 #
# ---------------------------------------------------------------------------- #
"""
    MarkupStyle

Holds information about the style specification set out by a `MarkupTag`.
"""
@with_kw mutable struct MarkupStyle
    default::Bool                               = false
    bold::Bool                                  = false
    dim::Bool                                   = false
    italic::Bool                                = false
    underline::Bool                             = false
    blink::Bool                                 = false
    inverse::Bool                               = false
    hidden::Bool                                = false
    striked::Bool                               = false
    color::Union{Nothing,AbstractColor}         = nothing
    background::Union{Nothing,AbstractColor}    = nothing
end

"""
    MarkupStyle(tag::MarkupTag)

Builds a MarkupStyle definition from a MarkupTag.
"""
function MarkupStyle(markup)
    codes = split(unspace_commas(markup))

    style = MarkupStyle()
    for code in codes
        if is_mode(code)
            setproperty!(style, Symbol(code), true)
        elseif is_color(code)
            style.color = get_color(code)
        elseif is_background(code)
            style.background = get_color(code; bg = true)
        elseif code != "nothing"
            @debug "Code type not recognized: $code" tag tag.markup typeof(code)
        end
    end
    return style
end

# -------------------------------- apply style ------------------------------- #

"""
    get_style_codes(style::MarkupStyle)

Get `ANSICode`s corresponding to a `MarkupStyle`.
"""
function get_style_codes(style::MarkupStyle)
    # start applying styles
    style_init, style_finish = "", ""
    for attr in fieldnames(MarkupStyle)
        value = getfield(style, attr)
        # BACKGROUND
        if attr == :background
            if !isnothing(value)
                code = ANSICode(value; bg = true)
            else
                code = nothing
            end

            # COLOR
        elseif attr == :color
            if !isnothing(value)
                code = ANSICode(value; bg = false)
            else
                code = nothing
            end

            # MODES
        elseif attr != :tag && value == true
            code = CODES[attr]

        else
            if value != false && attr != :tag
                @debug "Attr/value not recognized or not set" attr value
            end
            continue
        end

        if !isnothing(code)
            style_init *= code.open
            style_finish *= code.close
            style_finish *= (occursin(code.close, style_finish) ? "" : code.close)
        end
    end

    return style_init, style_finish
end

"""
    apply_style(text)

Apply style to a piece of text.

Extract markup style information and insert the
appropriate ANSI codes to style a string.
"""
function apply_style(text)::String
    has_markup(text) || return text

    while has_markup(text)
        # get opening markup tag
        open_match = match(OPEN_TAG_REGEX, text)
        markup = open_match.match[2:end-1]

        # get style codes
        ansi_open, ansi_close = get_style_codes(MarkupStyle(markup))

        # replace markup with ANSI codes
        text = replace_text(text, max(open_match.offset-1, 0), open_match.offset + length(markup)+1, ansi_open)

        # get closing tag (including [/] or missing close)
        close_rx = r"(?<!\[)\[(?!\[)\/" * markup * r"\]"

        if !occursin(close_rx, text)
            text = text * "[/" * markup * "]"
        end

        # @assert occursin(close_rx, text) "Did not find closing regex tag for $markup starting at position $(open_match.offset) in '$text'"
        close_match = match(close_rx, text)

        # check if there was previous ansi style info
        prev_style = get_last_ANSI_code(tview(text, 1, open_match.offset-1))
        prev_style = occursin(prev_style, ansi_close) ? "" : prev_style
        # replace close tag
        text = replace_text(text, close_match.offset-1, close_match.offset + length(markup)+2, ansi_close * prev_style)
    end
    return text
end

end


text:
export OrthoNormal

doc"""
    OrthoNormal
"""
struct OrthoNormal
end

function (o::OrthoNormal)(::Type{T}, dim1::Int, dim2::Int) where T
    I = eye(dim2)
    lr = 0.1
    eps = 0.05 / (dim1+dim2)
    tries = 0
    while tries < 10
        Q = randn(dim1, dim2) / sqrt(dim2)
        for i = 1:100
            QTQmI = Q' * Q - I
            loss = sum(QTQmI .^ 2 / 2)
            Q2 = Q .^ 2
            a = abs.(Q2 .+ sum(Q2,1) .+ sum(Q2,2) - 1.0) + eps
            Q -= lr * Q * QTQmI ./ a
            if maximum(Q) > 1e6 || loss > 1e6 || isinf(loss)
                tries += 1
                lr /= 2.0
                break
            end
        end
        return Matrix{T}(Q)
    end
    #Q = xp.random.randn(input_size, output_size) / xp.sqrt(output_size)
    throw("Generation failed.")
end


text:
using AbstractIndices, Test, Statistics

using AbstractIndices.IndexCore

import Base: OneTo

#include("uniqueness_tests.jl")
include("matmul_tests.jl")
include("math_tests.jl")
include("indexing_tests.jl")
include("abstractindex_tests.jl")
include("broadcasting.jl")
#include("indicesarray_tests.jl")
include("reduce_tests.jl")
#include("pop_tests.jl") TODO