text:
#=
Model:  fiveDimSys.jl

This is a five-dimensional model taken from [Girard05]. This also appears as
Example 4.1. page 57 in the thesis [LeGuernic09].

See also Fig. 5.1. page 71 in the same thesis.

[Girard05]       -- Girard, Antoine. "Reachability of uncertain linear systems using zonotopes."
                    International Workshop on Hybrid Systems: Computation and Control.
                    Springer, Berlin, Heidelberg, 2005.
[LeGuernic09]    -- Le Guernic, Colas. Reachability analysis of hybrid systems
                    with linear continuous dynamics. Diss.
                    Université Joseph-Fourier-Grenoble I, 2009.
=#
using Reachability, Plots

function compute(input_options::Pair{Symbol,<:Any}...)
    # =====================
    # Problem specification
    # =====================
    D = [-1 -4 0 0 0; 4 -1 0 0 0; 0 0 -3 1 0; 0 0 -1 -3 0; 0 0 0 0 -2.]
    P = [0.6 -0.1 0.1 0.7 -0.2; -0.5 0.7 -0.1 -0.8 0; 0.9 -0.5 0.3 -0.6 0.1;
        0.5 -0.7 0.5 0.6 0.3; 0.8 0.7 0.6 -0.3 0.2]
    A = P * D * inv(P)

    # initial set
    X0 = BallInf([1.0, 0.0, 0.0, 0.0, 0.0], 0.1)

    # input set
    U = Ball2(zeros(5), 0.01)

    # instantiate continuous LTI system
    S = ContinuousSystem(A, X0, U) # @system x' = A*x + u, u ∈ U

    # ===============
    # Problem solving
    # ===============
    # define solver-specific options
    options = merge(Options(
        :mode => "reach",
        :blocks => [1],
        :plot_vars => [0, 1]
        ), Options(input_options...))
    result = solve(S, options)

    # ========
    # Plotting
    # ========
    if options[:mode] == "reach"
        println("Plotting...")
        tic()
        plot(result)
        @eval(savefig(@filename_to_png))
        toc()
    end
end # function

compute(:N => 100, :T => 1.0); # warm-up
compute(:δ => 0.01, :T => 5.0); # benchmark settings (long)

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
#!/usr/bin/env julia

module EnsembleKalmanInversionTests

    #Start Test Script
    using EnsembleKalmanInversion
    using Test
    using Random
    using Statistics

    # Run tests

    @testset "EnsembleKalmanInversion Tests" begin

        @testset "EnsembleKalmanInversion Tests" begin
            @testset "Linear Test" begin
            #This should certainly be a set of "convergence" tests
                G = randn(50,100)
                x = randn(100)
                σ = 0.1
                yt = G*x
                y = yt .+ σ * randn(50)
                η = sqrt((y-yt)'*(y-yt)/σ^2)
                prior() = randn(100)
                gmap(q) = G*q
                J = 100
                N = 100
                xe = eki(y, σ, η, J, N, prior, gmap; verbosity=1, ρ=0.5, ζ=2.0)
                @test sqrt((yt-G*xe)'*((yt-G*xe))/σ^2) <= 2.0*η
            end

            @testset "Linear Batched Test" begin
            #This should certainly be a set of "convergence" tests
                G1 = randn(10,100)
                G2 = randn(15,100)
                G3 = randn(13,100)
                G4 = randn(11,100)

                x = randn(100)
                σ = 0.1
                yt = [G1*x,G2*x,G3*x,G4*x]
                y = [G1*x.+σ*randn(10),G2*x.+σ*randn(15),G3*x.+σ*randn(13),G4*x.+σ*randn(11)]
                η = sqrt((y-yt)'*(y-yt)/σ^2)
                prior() = randn(100)
                function gmap(q,b)
                    if b <=1
                        return G1*q
                    elseif b <= 2
                        return G2*q
                    elseif b <= 3
                        return G3*q
                    else
                        return G4*q
                    end
                end
                J = 100
                N = 100
                xe = eki(y, σ, η, J, N, prior, gmap; verbosity=1, ρ=0.5, ζ=2.0, batched=true, batches=3)
                ye = [G1*xe,G2*xe,G3*xe,G4*xe]
                @test sqrt((vcat(yt...)-vcat(ye...))'*((vcat(y...)-vcat(ye...)))/σ^2) <= 2.0*η
            end
        end
    end
end

ext: jl


text:
using ExecutableSpecifications.Gherkin
using ExecutableSpecifications:
    ResultAccumulator, accumulateresult, issuccess, featureresults,
    FeatureResult, ScenarioResult, Given, SuccessfulStepExecution, Scenario, StepFailed

@testset "Result Accumulator" begin
    @testset "Accumulate results; One feature with a successful scenario; Total result is success" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), [])
        given = Given("some precondition")
        successfulscenario = ScenarioResult([SuccessfulStepExecution()], Scenario("Some scenario", [], [given]))
        featureresult = FeatureResult(feature, [successfulscenario])

        accumulateresult(accumulator, featureresult)

        @test issuccess(accumulator)
    end

    @testset "Accumulate results; One feature with a failing scenario; Total result is fail" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), [])
        given = Given("some precondition")
        scenario = ScenarioResult([StepFailed("")], Scenario("Some scenario", [], [given]))
        featureresult = FeatureResult(feature, [scenario])

        accumulateresult(accumulator, featureresult)

        @test issuccess(accumulator) == false
    end

    @testset "Accumulate results; One scenario with one successful and one failing step; Total result is fail" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), [])
        given = Given("some precondition")
        when = When("some action")
        scenario = ScenarioResult([SuccessfulStepExecution(), StepFailed("")], Scenario("Some scenario", [], [given, when]))
        featureresult = FeatureResult(feature, [scenario])

        accumulateresult(accumulator, featureresult)

        @test issuccess(accumulator) == false
    end

    @testset "Accumulate results; Two scenarios, one failing; Total result is fail" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), [])
        given = Given("some precondition")
        when = When("some action")
        successfulscenario = ScenarioResult([SuccessfulStepExecution()], Scenario("Some scenario", [], [given]))
        failingscenario = ScenarioResult([StepFailed("")], Scenario("Some other scenario", [],  [when]))
        featureresult = FeatureResult(feature, [successfulscenario, failingscenario])

        accumulateresult(accumulator, featureresult)

        @test issuccess(accumulator) == false
    end

    @testset "Accumulate results; One feature, one passing Scenario; One success and zero failures" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), [])
        given = Given("some precondition")
        successfulscenario = ScenarioResult([SuccessfulStepExecution()], Scenario("Some scenario", [], [given]))
        featureresult = FeatureResult(feature, [successfulscenario])

        accumulateresult(accumulator, featureresult)

        @test featureresults(accumulator)[1].n_success == 1
        @test featureresults(accumulator)[1].n_failure == 0
    end

    @testset "Accumulate results; Two scenarios, one failing; One success and one failure" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), [])
        given = Given("some precondition")
        when = When("some action")
        successfulscenario = ScenarioResult([SuccessfulStepExecution()], Scenario("Some scenario", [], [given]))
        failingscenario = ScenarioResult([StepFailed("")], Scenario("Some other scenario", [],  [when]))
        featureresult = FeatureResult(feature, [successfulscenario, failingscenario])

        accumulateresult(accumulator, featureresult)

        @test featureresults(accumulator)[1].n_success == 1
        @test featureresults(accumulator)[1].n_failure == 1
    end

    @testset "Accumulate results; Seven scenarios, two failing; 5 success and 2 failures" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), [])
        given = Given("some precondition")
        when = When("some action")
        successfulscenario = ScenarioResult([SuccessfulStepExecution()], Scenario("Some scenario", [], [given]))
        failingscenario = ScenarioResult([StepFailed("")], Scenario("Some other scenario", [],  [when]))
        featureresult = FeatureResult(feature,
            [successfulscenario, successfulscenario, successfulscenario, successfulscenario, successfulscenario,
             failingscenario, failingscenario])

        accumulateresult(accumulator, featureresult)

        @test featureresults(accumulator)[1].n_success == 5
        @test featureresults(accumulator)[1].n_failure == 2
    end

    @testset "Accumulate results; Two successful features; Both have no failures" begin
        accumulator = ResultAccumulator()

        feature1 = Gherkin.Feature(FeatureHeader("1", [], []), [])
        feature2 = Gherkin.Feature(FeatureHeader("1", [], []), [])
        given = Given("some precondition")
        successfulscenario = ScenarioResult([SuccessfulStepExecution()], Scenario("Some scenario", [], [given]))
        featureresult1 = FeatureResult(feature1, [successfulscenario])
        featureresult2 = FeatureResult(feature1, [successfulscenario])

        accumulateresult(accumulator, featureresult1)
        accumulateresult(accumulator, featureresult2)

        @test featureresults(accumulator)[1].n_success > 0
        @test featureresults(accumulator)[1].n_failure == 0
        @test featureresults(accumulator)[2].n_success > 0
        @test featureresults(accumulator)[2].n_failure == 0
    end
end


text:
using JuMP, Cbc

sudokuModel = Model(Cbc.Optimizer)

gridsize = 3
@variable(sudokuModel, x[1:gridsize^2,1:gridsize^2,1:gridsize^2], binary=true)
a = zeros(gridsize^2,gridsize^2,gridsize^2)

as = [
1 4 2;1 5 6;1 7 7;1 9 1;
2 1 6;2 2 8;2 5 7;2 8 9;
3 1 1;3 2 9;3 6 4;3 7 5;
4 1 8;4 2 2;4 4 1;4 8 4;
5 3 4;5 4 6;5 6 2;5 7 9;
6 2 5;6 6 3;6 8 2;6 9 8;
7 3 9;7 4 3;7 8 7;7 9 4;
8 2 4;8 5 5;8 8 3;8 9 6;
9 1 7;9 3 3;9 5 1;9 6 8
]

for n in 1:size(as,1)
    c = as[n,:]
    a[c[1],c[2],c[3]]=1
end

for d1 = 1:gridsize^2
    for d2 = 1:gridsize^2
        @constraints(sudokuModel, begin
            sum(x[d1,d2,:])==1
            sum(x[d1,:,d2])==1
            sum(x[:,d1,d2])==1
        end)
        for k = 1:gridsize^2
            @constraint(sudokuModel, x[d1,d2,k]>=a[d1,d2,k])
        end
    end
end
for s1 = 0:(gridsize-1)
    for s2 = 0:(gridsize-1)
        for k = 1:gridsize^2
            @constraint(sudokuModel, sum(x[3*s1+a,3*s2+b,k] for a=1:gridsize, b=1:gridsize)==1)
        end
    end
end

stats = JuMP.optimize!(sudokuModel)

sol = JuMP.value.(x)
ABC = sum(sol[5,2,:].*collect(1:gridsize^2))+sum(sol[3,6,:].*collect(1:gridsize^2))+sum(sol[9,8,:].*collect(1:gridsize^2))

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
using NBodyIPs, JuLIP, BenchmarkTools, StaticArrays
using Test

using NBodyIPs: BondLengthDesc, transform, transform_d, fcut, fcut_d,
                invariants, invariants_d, invariants_ed
using NBodyIPs.Polys: NBPoly

nbasis = [0, 10, 10, 10, 10]
nbasis = [0, 10, 50, 100, 300]

println("Setting up the test systems ...")
r0 = rnn(:Cu)
TRANSFORM = "r -> exp( - 3 * ((r/$r0) - 1))"
at = rattle!(bulk(:Cu, cubic=true) * 2, 0.05)
rcut3 = 3.1 * r0
D3 = BondLengthDesc(TRANSFORM, (:cos, 0.66*rcut3, rcut3) )
rcut4 = 2.1 * r0
D4 = BondLengthDesc(TRANSFORM, (:cos, 0.66*rcut4, rcut4) )
rcut5 = 1.5 * r0
D5 = BondLengthDesc(TRANSFORM, (:cos, 0.66*rcut5, rcut5) )
DD = [nothing, D3, D3, D4, D5]

random_nbody(N, ntup) = (
   (N <= 3) ? NBPoly( [tuple( [rand(1:4, ((N*(N-1))÷2)); 0]... ) for n = 1:ntup],
                     (0.1+rand(ntup))/factorial(N), DD[N] )
           : NBPoly( [tuple( rand(0:N, ((N*(N-1))÷2+1))... ) for n = 1:ntup],
                   (0.1+rand(ntup))/factorial(N), DD[N] ) )

BB = [ nothing; [ [random_nbody(N, 1) for _=1:nbasis[N]] for N = 2:5 ] ]
at = rattle!(bulk(:Cu, cubic=true) * 2, 0.05)

println("Performance Tests")
println("-----------------")
println("to get a rough estimate of the cost per atom take the `separate`")
println("timing and multiply by the factor provided")
for N = 2:4
   B = BB[N]
   println("[$N-body], length(B) == $(length(B)), factor = $(1/(length(at)*length(B)))")
   print(" E separate:" ); @btime ([energy(b, $at) for b in $B])
   print(" E combined:" ); @btime energy( $B, $at )
   print(" F separate:" ); @btime ([forces(b, $at) for b in $B])
   print(" F combined:" ); @btime forces( $B, $at )
   print(" S separate:" ); @btime ([virial(b, $at) for b in $B])
   print(" S combined:" ); @btime virial( $B, $at )
end

# B = BB[4]
#
# # @time [energy(b, at) for b in B]
# # @time [energy(b, at) for b in B]
# # @profile [energy(b, at) for b in B]
# # Base.Profile.print()
#
# @time [forces(b, at) for b in B];
# @time [forces(b, at) for b in B];
# @profile [forces(b, at) for b in B];
# Base.Profile.print()

# # THE CURRENT IMPLEMENTATION
# [2-body], length(B) == 10, factor = 0.003125
#  E separate:  26.563 ms (8693 allocations: 9.80 MiB)
#  E combined:  3.061 ms (804 allocations: 1002.56 KiB)
#  F separate:  44.766 ms (403043 allocations: 21.90 MiB)
#  F combined:  18.543 ms (395712 allocations: 13.08 MiB)
#  S separate:  30.726 ms (11363 allocations: 9.98 MiB)
#  S combined:  4.431 ms (5300 allocations: 1.22 MiB)
# [3-body], length(B) == 50, factor = 0.000625
#  E separate:  1.611 s (43403 allocations: 49.01 MiB)
#  E combined:  218.385 ms (804 allocations: 1002.94 KiB)
#  F separate:  2.663 s (2015203 allocations: 109.47 MiB)
#  F combined:  593.160 ms (1975312 allocations: 61.48 MiB)
#  S separate:  2.593 s (56803 allocations: 49.89 MiB)
#  S combined:  502.200 ms (23260 allocations: 2.19 MiB)
# [4-body], length(B) == 100, factor = 0.0003125
#  E separate:  4.788 s (86403 allocations: 34.84 MiB)
#  E combined:  312.957 ms (800 allocations: 356.28 KiB)
#  F separate:  7.037 s (1297203 allocations: 72.09 MiB)
#  F combined:  859.382 ms (1217008 allocations: 37.69 MiB)
#  S separate:  6.862 s (113203 allocations: 36.33 MiB)
#  S combined:  807.974 ms (45706 allocations: 2.50 MiB)
# [5-body], length(B) == 300, factor = 0.00010416666666666667
#  E separate:  2.188 s (259203 allocations: 78.31 MiB)
#  E combined:  28.564 ms (800 allocations: 268.41 KiB)
#  F separate:  2.774 s (1472403 allocations: 115.96 MiB)
#  F combined:  124.575 ms (1230208 allocations: 38.18 MiB)
#  S separate:  2.926 s (339604 allocations: 82.50 MiB)
#  S combined:  79.406 ms (135507 allocations: 6.44 MiB)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
struct Simulation{DT <: Number, IT <: Integrator{DT}}
    int::IT
    x::OffsetArray{DT,3,Array{DT,3}}

    function Simulation(x₀::AbstractArray{DT}, int::IT, nt::Int) where {DT, IT <: Integrator{DT}}
        x = OffsetArray(zeros(DT, size(x₀)..., nt+1), axes(x₀)..., 0:nt)
        x[:,:,0] .= x₀
        new{DT,IT}(int, x)
    end
end

function Simulation(equ::Equation{DT}, integrator::Type{<:Integrator}, Δt::DT, nt::Int; kwargs...) where {DT}
    int = integrator(equ, Δt; kwargs...)
    Simulation(equ.x₀, int, nt)
end

ndims(sim::Simulation) = ndims(sim.int.equ)
nsamples(sim::Simulation) = nsamples(sim.int.equ)
ntimesteps(sim::Simulation) = lastindex(sim.x,3)

eachsample(sim::Simulation) = axes(sim.x,2)
eachtimestep(sim::Simulation) = axes(sim.x,3)[1:end]

function (sim::Simulation{DT})() where {DT}
    for n in eachtimestep(sim)
        for i in eachsample(sim)
            @views sim.int(sim.x[:,i,n-1], sim.x[:,i,n])
        end
    end
    return sim.x
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
@cell "Layouting" [scatter, lines, surface, heatmap] begin
    p1 = scatter(rand(10), markersize = 1)
    p2 = lines(rand(10), rand(10))
    p3 = surface(0..1, 0..1, rand(100, 100))
    p4 = heatmap(rand(100, 100))
    x = 0:0.1:10
    p5 = lines(0:0.1:10, sin.(x))
    pscene = vbox(
        hbox(p1, p2),
        p3,
        hbox(p4, p5, sizes = [0.7, 0.3]),
        sizes = [0.2, 0.6, 0.2]
    )
end


text:
module Sim

    using DifferentialEquations

    function solveT(eqs, tspan)
        sys = ODESystem(eqs)
        tspan = (0.0,100.0)

        prob = ODEProblem(sys,u0,tspan,p,jac=true)
        sol = solve(prob)

        return sol
    end

    function solveNullcline(eqs)

        ode_f = ODEFunction(sys)

        substitute.(eqs[1].rhs,p)
    end

end

max_stars_repo_name: stevengogogo/ProgrammableAging

max_issues_repo_name: stevengogogo/ProgrammableAging

max_forks_repo_name: stevengogogo/ProgrammableAging


text:
module TestsInSrc

const TEMPLATE_PATH = joinpath(@__DIR__, "templates", "tests_in_src.jl")

end # module


text:
# ---
# title: 523. Continuous Subarray Sum
# id: problem523
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/continuous-subarray-sum/description/>
# hidden: true
# ---
#
# Given a list of **non-negative** numbers and a target **integer** k, write a
# function to check if the array has a continuous subarray of size at least 2
# that sums up to a multiple of **k** , that is, sums up to n*k where n is also
# an **integer**.
#
#
#
# **Example 1:**
#
#
#
#     Input: [23, 2, 4, 6, 7],  k=6
#     Output: True
#     Explanation: Because [2, 4] is a continuous subarray of size 2 and sums up to 6.
#
#
# **Example 2:**
#
#
#
#     Input: [23, 2, 6, 4, 7],  k=6
#     Output: True
#     Explanation: Because [23, 2, 6, 4, 7] is an continuous subarray of size 5 and sums up to 42.
#
#
#
#
# **Constraints:**
#
#   * The length of the array won't exceed 10,000.
#   * You may assume the sum of all the numbers is in the range of a signed 32-bit integer.
#
#
## @lc code=start
using LeetCode

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


text:
import Base.+, Base.-, Base./, Base.*, Base.^
import Base.sort
import SchumakerSpline.evaluate
abstract type UnivariateFunction end

struct Undefined_Function <: UnivariateFunction
end
Base.broadcastable(e::Undefined_Function) = Ref(e)
struct PE_Function{F<:Real,I<:Integer} <: UnivariateFunction
    # a exp(b(x-base)) (x-base)^d
    a_::F
    b_::F
    base_::F
    d_::I
    function PE_Function(a_::R,b_::S,base_::T,d_::I) where R<:Real where S<:Real where T<:Real where I<:Integer
        promo_type = promote_type(T,S,T)
        if d_ < 0
            error("Negative polynomial powers are not supported by this package")
            # These are banned due to the complications for calculus. Most
            # divisions are not allowed for the same reason.
        elseif (abs(b_) < tol) & (d_ == 0)
            return new{promo_type,I}(promo_type(a_),promo_type(0.0),promo_type(0.0),I(0))
        elseif abs(a_) < tol
            return new{promo_type,I}(promo_type(0.0),promo_type(0.0),promo_type(0.0),I(0))
        else
            return new{promo_type,I}(promo_type(a_),promo_type(b_),promo_type(base_),I(d_))
        end
    end
    function PE_Function(a_::R,b_::S,base_::Union{Date,DateTime},d_::I) where R<:Real where S<:Real where I<:Real
        new_base_ = years_from_global_base(base_)
        return PE_Function(a_, b_, new_base_, d_)
    end
    function PE_Function{Q,N}(a_::R,b_::S,base_::T,d_::I) where Q<:Real where N<:Integer where R<:Real where S<:Real where T<:Real where I<:Integer
        return new{Q,N}(Q(a_), Q(b_), Q(base_), N(d_))
    end
end
Base.broadcastable(e::PE_Function) = Ref(e)
struct Sum_Of_Functions <: UnivariateFunction
    functions_::Union{Array{PE_Function,1},Array{PE_Function{<:Real,<:Integer},1}}
    function Sum_Of_Functions(funcs::Sum_Of_Functions)
        return funcs
    end
    function Sum_Of_Functions(funcs::AbstractArray)
        undefined_funcs  = funcs[isa.(funcs, Ref(Undefined_Function))]
        if length(undefined_funcs) > 0
            return Undefined_Function()
        end
        functions_ = clean_array_of_functions(funcs)
        return new(functions_)
    end
end
Base.broadcastable(e::Sum_Of_Functions) = Ref(e)
struct Piecewise_Function <: UnivariateFunction
    starts_::Array{<:Real,1}
    functions_::Array{UnivariateFunction,1}
    function Piecewise_Function(starts::Array{<:Real,1}, functions::Array)
        if !issorted(starts)
            error("Piecewise_Function must be constructed with a sorted increasing list of starts and the corresponding functions.")
        end
        if starts[1] != -Inf
            starts    = vcat(-Inf, starts)
            functions = vcat(Undefined_Function(), functions)
        end
        len = length(starts)
        if starts[len] == Inf
            starts    = starts[1:(len-1)]
            functions = functions[1:(len-1)]
        end
        starts_without_pw_parts, functions_without_pw_parts = deal_with_piecewise_inputs(starts, functions)
        new(starts_without_pw_parts, functions_without_pw_parts)
    end
    function Piecewise_Function(starts::Union{Array{DateTime,1},Array{Date,1},Array{Union{DateTime,Date},1}}, functions::Array)
        starts_ = years_from_global_base.(starts)
        return Piecewise_Function(starts_, functions)
    end
    function Piecewise_Function(start::Real, func)
        return Piecewise_Function([start], [func])
    end
    function Piecewise_Function(start::Union{Date,DateTime}, func)
        start_float = years_from_global_base(start)
        return Piecewise_Function(start_float, func)
    end
end
Base.broadcastable(e::Piecewise_Function) = Ref(e)

function sort(funcs::Array{PE_Function,1})
    len = length(funcs)
    attributes = hcat(map(f -> f.d_, funcs), map(f -> f.base_, funcs), map(f -> f.b_, funcs) )
    ordering = convert(Array{Int}, sortslices(hcat(attributes, 1:len) , dims = 1)[:,4])
    return funcs[ordering]
end

function trim_piecewise_function(func::Piecewise_Function, left_limit::Real, right_limit::Real)
    if func.starts_[1] != -Inf
        starts_ = [-Inf, left_limit, right_limit]
        funcs_  = [Undefined_Function(), func, Undefined_Function()]
    else
        starts_ = [left_limit, right_limit]
        funcs_  = [func, Undefined_Function()]
    end
    return Piecewise_Function(starts_, funcs_)

end

function deal_with_piecewise_inputs(starts::Array{<:Real,1}, functions::Array)
    where_are_pw_bits = findall(typeof.(functions) .== UnivariateFunctions.Piecewise_Function)
    if length(where_are_pw_bits) < 1
        return starts, functions
    else
        first_pw_slice_i = where_are_pw_bits[1]
        pw = functions[first_pw_slice_i]
        from = starts[first_pw_slice_i]
        if length(starts) > first_pw_slice_i
            to = starts[first_pw_slice_i+1]
        else
            to = Inf
        end
        pw_starts, pw_funcs = take_piecewise_slice(pw.starts_, pw.functions_,from, to)
        pw_starts[1] = from
        new_starts = vcat(first_entries(starts, first_pw_slice_i-1), pw_starts, last_entries(starts, length(starts) - first_pw_slice_i))
        new_funcs = vcat(first_entries(functions, first_pw_slice_i-1), pw_funcs, last_entries(functions, length(starts) - first_pw_slice_i))
        return deal_with_piecewise_inputs(new_starts, new_funcs)
    end

end

function first_entries(vec::Array, how_many::Integer)
    if how_many < 1
        return Array{Any,1}()
    elseif how_many > length(vec)
        return vec
    else
        return vec[1:how_many]
    end
end

function last_entries(vec::Array, how_many::Integer)
    if how_many < 1
        return Array{Any,1}()
    elseif how_many > length(vec)
        return vec
    else
        leng =  length(vec)
        return vec[(leng-how_many+1):leng]
    end
end

function take_piecewise_slice(starts::Array{<:Real,1}, functions::Array, from::Real, to::Real)
    from_i = searchsortedlast(starts, from)
    to_i = searchsortedlast(starts, to)
    return starts[from_i:to_i], functions[from_i:to_i]
end

function find(a::BitArray)
    len = length(a)
    indices = Array{Int,1}()
    for i in 1:len
        if a[i]
            append!(indices, i)
        end
    end
    return indices
end

function rationalise_array_of_functions(funcs::Array{PE_Function,1})
    len = length(funcs)
    if len < 2
        return funcs
    end
    attributes = hcat(map(f -> f.b_, funcs), map(f -> f.base_, funcs), map(f -> f.d_, funcs) )
    multipliers = map(f -> f.a_, funcs)
    unique_attributes = unique(attributes, dims = 1)
    new_len = size(unique_attributes)[1]
    new_funcs = Array{PE_Function,1}(undef,new_len)
    for i in 1:new_len
        atts =  transpose(unique_attributes[i,:])
        which_multipliers = all(attributes .== atts, dims=2)
        multiplier = sum(multipliers[find(which_multipliers)])
        new_funcs[i] = PE_Function(multiplier, atts[1], atts[2], convert(Int, atts[3]))
    end
    return new_funcs
end

function clean_array_of_functions(funcs::Array)
    undefined_funcs  = funcs[isa.(funcs, Ref(Undefined_Function))]
    piecewise_funcs  = funcs[isa.(funcs, Ref(Piecewise_Function))]
    if length(undefined_funcs) > 0
        return Undefined_Function()
    elseif length(piecewise_funcs) > 0
        error("You cannot directly construct a Sum_Of_Functions with a Piecewise_Function in the input array. Instead add up the piecewise functions directly, for instance typing 'f1 + f2' for the two piecewise functions.")
    end
    pe_funcs  = funcs[isa.(funcs, Ref(PE_Function))]
    pe_funcs = pe_funcs[abs.(map( x -> x.a_, pe_funcs)) .>= eps()]
    sum_funcs = funcs[isa.(funcs, Ref(Sum_Of_Functions))]
    if length(sum_funcs) > 0
        for sf in sum_funcs
            pe_funcs_in_sf = clean_array_of_functions(sf.functions_)
            pe_funcs       = vcat(pe_funcs, pe_funcs_in_sf)
        end
    end
    if length(pe_funcs) == 0
        return [PE_Function(0.0,0.0,0.0,0)]
    end
    simplified_functions = rationalise_array_of_functions(convert(Array{PE_Function,1}, pe_funcs))
    return simplified_functions
end

function +(number::Real, f::UnivariateFunction)
    return +(f,number)
end
function -(number::Real, f::UnivariateFunction)
    return +(number, -1*f)
end
function *(number::Real, f::UnivariateFunction)
    return *(f,number)
end
function /(number::Real, f::UnivariateFunction)
    error("It is not possible yet to divide scalars by UnivariateFunctions")
end
function ^(number::Real, f::UnivariateFunction)
    error("It is not possible yet to raise to the power of a UnivariateFunctions")
end

function evaluate(f::UnivariateFunction, d::Union{Date,DateTime})
    date_in_relation_to_global_base = years_from_global_base(d)
    return evaluate(f, date_in_relation_to_global_base)
end
function evaluate(f::UnivariateFunction, x::DatePeriod)
    return evaluate(f, period_length(x))
end

function ^(f1::UnivariateFunction,num::Integer) # This will get overridden for undefined and zeros.
    if num < 0
        error("Cannot raise any univariate function to a negative power")
    elseif num == 0
        return PE_Function(1.0,0.0,0.0,0)
    elseif num == 1
        return f1
    elseif num == 2
        return f1 * f1
    else
        product = f1 * f1
        for i in 1:(num-2)
            product = product * f1
        end
        return product
    end
end

function change_base_of_PE_Function(f::PE_Function, new_base::Real)
    old_base = f.base_
    diff = new_base - old_base
    if abs(diff) < tol
        return f
    end
    # First the exponential part.
    new_a = f.a_ * exp(f.b_*diff)
    if new_a < tol
        error("Underflow problem. Changing to this base cannot be done")
    end
    # Now the polynomial part.
    if f.d_ == 0
        return PE_Function(new_a, f.b_, new_base, 0)
    else
        n = f.d_
        funcs = Array{UnivariateFunction,1}(undef,n+1)
        for r in 0:n
            binom_coeff = factorial(n) / (factorial(r) * factorial(n-r))
            new_multiplier = binom_coeff * new_a * diff^r
            new_func = PE_Function(new_multiplier, f.b_, new_base, n-r )
            funcs[r+1] = new_func
        end
        return Sum_Of_Functions(funcs)
    end
end

# Conversions for linearly rescaling inputs.
function convert_to_linearly_rescale_inputs(f::Undefined_Function, alpha::Real, beta::Real)
    return f
end
function convert_to_linearly_rescale_inputs(f::PE_Function, alpha::Real, beta::Real)
    # We want the change the function so that whenever we put in x it is like we put in alpha x + beta.
    beta = beta / alpha
    alpha = 1.0/alpha
    new_base_ = (f.base_ + beta)/alpha
    new_multiplier = f.a_ * alpha^(f.d_)
    new_power_ = f.b_ * alpha
    return PE_Function(new_multiplier, new_power_, new_base_, f.d_)
end
function convert_to_linearly_rescale_inputs(f::Sum_Of_Functions, alpha::Real, beta::Real)
    funcs = convert_to_linearly_rescale_inputs.(f.functions_, alpha,beta)
    return Sum_Of_Functions(funcs)
end
function convert_to_linearly_rescale_inputs(f::Piecewise_Function, alpha::Real, beta::Real)
    new_starts_ = (f.starts_ .* alpha) .+ beta
    funcs_ = convert_to_linearly_rescale_inputs.(f.functions_, alpha,beta)
    return Piecewise_Function(new_starts_, funcs_)
end


text:
#suite["edges_of"] = BenchmarkGroup(["operator"])
#suite["edges_of"]["graphs"] = BenchmarkGroup()
#suite["edges_of"]["digraphs"] = BenchmarkGroup()
#for (name,g) in GRAPHS
#    suite["edges_of"]["graphs"][name] = @benchmarkable edges_of($g,$1)
#end
#for (name,g) in DIGRAPHS
#    suite["edges_of"]["digraphs"][name] = @benchmarkable edges_of($g,$1)
#end


text:
mutable struct SparseSymmetricMatrix{T<:Number}
    entries::Dict{Tuple{Any, Any}, T}
    indices::Set{Any}
end

SparseSymmetricMatrix(T::Type) = SparseSymmetricMatrix{T}(Dict{Tuple{Any,Any},T}(), Set())

SparseSymmetricMatrix() = SparseSymmetricMatrix(Float64)

entries(m::SparseSymmetricMatrix) = m.entries

indices(m::SparseSymmetricMatrix) = m.indices

function setindex!(m::SparseSymmetricMatrix{Float64}, v::Float64, i, j)
    entries(m)[(i, j)] = v
    push!(indices(m), i)
    push!(indices(m), j)
end

function getindex(m::SparseSymmetricMatrix{Float64}, i, j)
    get(m.entries, (i, j), zero(Float64))
end

# #  this code is deprecated
# start(m::SparseSymmetricMatrix) = start(m.entries)
# done(m::SparseSymmetricMatrix, state) = done(m.entries, state)
# next(m::SparseSymmetricMatrix, state) = next(m.entries, state)

function size(m::SparseSymmetricMatrix)
    is = Set()
    js = Set()
    for (i, j) in keys(m.entries)
        is = union(is,i)
        js = union(js,j)
    end
    return max(length(is), length(js))
end

copy(m::SparseSymmetricMatrix) = SparseSymmetricMatrix(copy(entries(m)), copy(indices(m)))

function *(x::Float64, m::SparseSymmetricMatrix{Float64})
    p = copy(m)
    for k in keys(entries(m))
        entries(p)[k] *= x
    end
    p
end

function delete!(m::SparseSymmetricMatrix{Float64}, i, j)
    delete!(entries(m), (i, j))
    igone = true
    jgone = true
    for (ti, tj) in keys(entries(m))
        if i == ti
            igone = false
        end
        if j == tj
            jgone = false
        end
    end
    if igone
        delete!(indices(m), i)
    end
    if jgone && contains(indices(m), j)
        delete!(indices(m), j)
    end
end

function ==(A::SparseSymmetricMatrix, B::SparseSymmetricMatrix)
    A.entries == B.entries && A.indices == B.indices
end

function isapprox(A::SparseSymmetricMatrix, B::SparseSymmetricMatrix)
    if indices(A) != indices(B)
        return false
    end
    for ((ri, ci), v) in A.entries
        if !isapprox(v, B[ri, ci])
            return false
        end
    end
    for ((ri, ci), v) in B.entries
        if !isapprox(v, A[ri, ci])
            return false
        end
    end
    true
end

max_stars_repo_name: JAndriesJ/SemidefiniteProgramming

max_issues_repo_name: JAndriesJ/SemidefiniteProgramming

max_forks_repo_name: JAndriesJ/SemidefiniteProgramming


text:
module ParameterizedNotebooks

using MacroTools

export @nbarg, @nbonly, @nbreturn, ParameterizedNotebook

macro nbarg(ex)
    @capture(ex, arg_ = val_) || error("usage: @nbarg name = val")
    esc(ex)
end

macro nbonly(ex)
    esc(ex)
end

macro nbreturn(ex=nothing)
    :(ReturnValue($(esc(ex))))
end

struct ReturnValue
    val
end

struct ParameterizedNotebook
    filename :: String
    parameters :: Vector{Symbol}
    exprs :: Vector
end

function ParameterizedNotebook(filename::String)
    script = read(`jupyter nbconvert --to script $filename --stdout`, String)
    exprs = parsecode(script)
    parameters = Symbol[]
    for ex in exprs
        if @capture(ex, @nbarg name_ = val_)
            push!(parameters, name)
        end
    end
    ParameterizedNotebook(filename, parameters, exprs)
end

function Base.show(io::IO, nb::ParameterizedNotebook)
    print(io, "ParameterizedNotebook(\"", nb.filename, "\") with parameters: ")
    print(io, join(string.(nb.parameters), ", "))
end

function (nb::ParameterizedNotebook)(; kwargs...)
    kwargs = Dict(kwargs)
    for name in nb.parameters
        haskey(kwargs, name) || error("Must provide keyword argument for $name.")
    end
    for ex in nb.exprs
        if @capture(ex, @nbarg name_ = val_)
            @eval Main $name = $(kwargs[name])
        elseif @capture(ex, @nbonly _)
            continue
        else
            ans = Main.eval(ex)
            if ans isa ReturnValue
                return ans.val
            end
        end
    end
end

function parsecode(code::String)
    # https://discourse.julialang.org/t/parsing-a-julia-file/32622
    filter(
        x -> !(x isa LineNumberNode),
        Meta.parse(join(["quote", code, "end"], ";")).args[1].args
    )
end

end


text:
using CBinding
using Libdl
using Clang_jll

const LIBCLANG_PATH = get(ARGS, 1, Clang_jll.libclang_path)
const LIBLLVM_PATH = get(ARGS, 2, joinpath(dirname(Base.julia_cmd().exec[1]), Base.LIBDIR, "julia", "libLLVM.so"))

withenv("LD_LIBRARY_PATH" => get(ENV, "LD_LIBRARY_PATH", "")*":"*normpath(dirname(LIBLLVM_PATH))) do
	dlopen(LIBLLVM_PATH) do lib  # needed to hijack the libLLVM that clank will be linked to
		@eval c`-std=c99 -I$(dirname(dirname(LIBCLANG_PATH)))/include -L$(dirname(LIBCLANG_PATH)) -lclang`
		@eval c"""
			#include <time.h>
		"""s

		libclang = quote
			const basedir = "$(CBinding.dirname(CBinding.dirname(CBinding.libclang_path())))"
			const incdir  = "$(basedir)/include"
			const libdir  = "$(basedir)/$(Sys.iswindows() ? "bin" : "lib")"
			const libpath = Symbol("$(libdir)/libclang")

			const var"c\"time_t\"" = $(Cvoid)  # HACK: this is not used anyways
			const var"c\"size_t\"" = $(Csize_t)

			$(@eval @macroexpand c"""
				#include <clang-c/Documentation.h>
				#include <clang-c/Index.h>
			"""jium)
		end

		function sanitize(str::String)
			str = replace(str, r"(\W)"*dirname(dirname(LIBCLANG_PATH))*"/" => s"\1$(basedir)/")
			str = replace(str, r"(\W)"*"\$(basedir)/lib/" => s"\1$(libdir)/")
			str = replace(str, r"(\W)"*"\$(basedir)/include/" => s"\1$(incdir)/")
			str = replace(str, "Symbol(\"\$(libdir)/libclang\")," => "libpath,")
			str = replace(str, r"(\W)"*dirname(@__DIR__)*"/" => s"\1../")

			lines = []
			for line in split(str, '\n')
				l = strip(line)
				startswith(l, "#=") && endswith(l, "=#") && continue

				m = match(r"/(TEST|C|C\.jl|home|usr)/", l)
				isnothing(m) || error("Unsanitized paths still remain:  $(l)")

				push!(lines, line)
			end
			return join(lines, '\n')
		end

		libclang = sanitize(string(libclang))  # NOTE: do this _before_ we truncate the file!
		open(joinpath(@__DIR__, "libclang-$(CBinding.libclang_version(LIBCLANG_PATH)).jl"), "w+") do io
			println(io, libclang)
		end
	end
end


text:
# Make a plot of the periods of the planets versus integer values:
using PyPlot; using DelimitedFiles
using Statistics

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

data = readdlm("../../tables/times_forecast.txt",',')

# Now compute the periods and uncertainties:
period = Float64[]
sigp   = Float64[]
nplanet = 7
for i=1:nplanet
  indx = data[:,1] .== i
  nt = sum(indx)
  fn = zeros(2,nt)
  fn[1,:] .= 1.0
  fn[2,:] = data[indx,2]
  coeff,cov = regress(fn,data[indx,3],data[indx,4])
  push!(period,coeff[2])
  push!(sigp,sqrt(cov[2,2]))
end

#period = [ 1.510826, 2.421937, 4.049219, 6.101013, 9.207540, 12.352446, 18.772866]
#sigp = [0.000006, 0.000018, 0.000026, 0.000035, 0.000032, 0.000054, 0.000214]
i0 = [24,15,9,6,4,3,2]
n_i = 1 ./ period
sign = sigp .* n_i ./period
# Now, regress:
#fn = zeros(3,nplanet)
fn = zeros(2,nplanet)
fn[1,:] .= 1.0
fn[2,:] = i0
#fn[3,:] = i0.^2
coeff,cov = regress(fn,n_i,sign)
#errorbar(i0,period,sigp,fmt=".")
#errorbar(i0,n_i,sigp,fmt=".")
#errorbar(i0,n_i,sigp,fmt=".")
rc("xtick",labelsize=15)
rc("ytick",labelsize=15)
fig,axes = subplots(2,1,sharex="col")
ax = axes[1]
#ax.errorbar(i0,n_i .- coeff[2] .* i0,sigp,fmt=".")
ax.errorbar(i0,n_i,sign,fmt=".",label=L"$P^{-1}$")
#errorbar(i0,n_i .- coeff[2] .* i0 .- coeff[3] .* i0.^2,sigp,fmt=".")
#plot(i0,1 ./(coeff[1] .+ coeff[2] .* i0))
ax.plot([0;i0;25],[coeff[1];(coeff[1] .+ coeff[2] .* i0);coeff[1]+coeff[2]*25],label=L"$nP_0^{-1} - P_{TTV}^{-1}$")
#ax.plot([0,25],[coeff[1],coeff[1]])
#plot([0;i0;25],[0;coeff[2] .*i0; coeff[2]*25],linestyle=":")
#plot([0;i0;25],[0;coeff[2] .*i0; coeff[2]*25],linestyle=":")
pname = ["b","c","d","e","f","g","h"]
for i=1:nplanet
#  text(i0[i]+1,period[i]+0.1,string(pname[i],": ",i0[i]))
  ax.text(i0[i],n_i[i]-0.07,string(i0[i]))
  ax.text(i0[i],n_i[i]+0.03,string(pname[i]))
end
#ax.set_xlabel(L"n: Integer multiple of frequency $P_0^{-1}$")
ax.set_ylabel("1/P [1/days]",fontsize=15)
ax.legend()
#axis([0,25,0,20])
ax.axis([0,25,-0.1,0.8])
ax.plot([0,25],[0,0],linestyle=":")
#ax.text(2,0.6,L"$P^{-1} = n P_0^{-1} - P_{TTV}^{-1}$",fontsize=15)
#ax.text(5,0.65,L"$P_0 = 36.15007\pm 0.00013\ \mathrm{d} $",fontsize=10)
#ax.text(5,0.60,L"$P_{TTV} = 492.3\pm 0.1\ \mathrm{d} $",fontsize=10)
resid = (period .- 1 ./(coeff[1] .+ coeff[2] .* i0)) ./period .* 100
println("Residuals: ",resid," % ",std(resid))

ax = axes[2]
ax.errorbar(i0,n_i .- coeff[2] .* i0,sign,fmt=".",label=L"$P^{-1} - nP_0^{-1}$")
#ax.errorbar(i0,n_i .- coeff[2] .* i0 .- coeff[3] .* i0.^2,sign,fmt=".",label=L"$P^{-1} - nP_0^{-1}$")
ax.plot([0,25],[coeff[1],coeff[1]],label=L"$-P_{TTV}^{-1}$")
#ax.plot([0,25],[coeff[1]-sqrt(cov[1,1]),coeff[1]-sqrt(cov[1,1])],linestyle=":")
#ax.plot([0,25],[coeff[1]+sqrt(cov[1,1]),coeff[1]+sqrt(cov[1,1])],linestyle=":")
ax.set_xlabel(L"n: Integer multiple of frequency $P_0^{-1}$",fontsize=15)
ax.set_ylabel(L"1/P [1/days] - $nP_0^{-1}$",fontsize=15)
ax.legend()
ax.plot([0,25],[0,0],linestyle=":")
ax.axis([0,25,-0.003,0.001])
tight_layout()
subplots_adjust(hspace = 0)
savefig("T1_orbital_frequency_relation.png",bbox_inches="tight")


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

"""
    UniversalKriging(γ, degree, dim)
    UniversalKriging(data, γ, degree)

Universal Kriging with variogram model `γ` and polynomial
`degree` on a spatial domain of dimension `dim`.

Optionally, pass the geospatial `data` to the [`fit`](@ref) function.

### Notes

* [`OrdinaryKriging`](@ref) is recovered for 0th degree polynomial
* For non-polynomial mean, see [`ExternalDriftKriging`](@ref)
"""
struct UniversalKriging{G<:Variogram} <: KrigingEstimator
  γ::G
  degree::Int
  dim::Int
  exponents::Matrix{Int}

  function UniversalKriging{G}(γ, degree, dim) where {G<:Variogram}
    @assert degree ≥ 0 "degree must be nonnegative"
    @assert dim > 0 "dimension must be positive"
    exponents = UKexps(degree, dim)
    new(γ, degree, dim, exponents)
  end
end

UniversalKriging(γ, degree, dim) = UniversalKriging{typeof(γ)}(γ, degree, dim)

UniversalKriging(data::Data, γ, degree) =
  GeoStatsBase.fit(UniversalKriging(γ, degree, embeddim(data)), data)

function UKexps(degree::Int, dim::Int)
  # multinomial expansion
  expmats = [hcat(collect(multiexponents(dim, d))...) for d in 0:degree]
  exponents = hcat(expmats...)

  # sort expansion for better conditioned Kriging matrices
  sorted = sortperm(vec(maximum(exponents, dims=1)), rev=true)

  exponents[:,sorted]
end

nconstraints(estimator::UniversalKriging) = size(estimator.exponents, 2)

function set_constraints_lhs!(estimator::UniversalKriging, LHS::AbstractMatrix, domain)
  exponents = estimator.exponents
  nobs = nelements(domain)
  nterms = size(exponents, 2)

  # set polynomial drift blocks
  for i in 1:nobs
    x = coordinates(centroid(domain, i))
    for j in 1:nterms
      LHS[nobs+j,i] = prod(x.^exponents[:,j])
      LHS[i,nobs+j] = LHS[nobs+j,i]
    end
  end

  # set zero block
  LHS[nobs+1:end,nobs+1:end] .= zero(eltype(LHS))

  nothing
end

factorize(::UniversalKriging, LHS::AbstractMatrix) = bunchkaufman(Symmetric(LHS), check=false)

function set_constraints_rhs!(fitted::FittedKriging{<:UniversalKriging}, uₒ)
  exponents = fitted.estimator.exponents
  RHS = fitted.state.RHS
  nobs = nelements(fitted.state.data)
  nterms = size(exponents, 2)

  # set polynomial drift
  xₒ = coordinates(centroid(uₒ))
  for j in 1:nterms
    RHS[nobs+j] = prod(xₒ.^exponents[:,j])
  end

  nothing
end


text:
@testset "Parallel.Random Vertex Cover" begin

    g0 = SimpleGraph(0)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g0)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            @test isempty(c)
        end
    end

    g1 = SimpleGraph(1)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g1)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            @test isempty(c)
        end
    end

    add_edge!(g1, 1, 1)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g1)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            @test c == [1,]
        end
    end

    g3 = StarGraph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g3)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            @test (length(c)== 2 && (c[1] == 1 || c[2] == 1))
        end
    end

    g4 = CompleteGraph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g4)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            @test length(c)== 4 #All except one vertex
        end
    end

    g5 = PathGraph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g5)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            sort!(c)
            @test (c == [1, 2, 3, 4] || c == [1, 2, 4, 5] || c == [2, 3, 4, 5])
        end
    end

    add_edge!(g5, 2, 2)
    add_edge!(g5, 3, 3)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g5)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            sort!(c)
            @test (c == [1, 2, 3, 4] || c == [1, 2, 3, 4, 5] || c == [2, 3, 4] || c == [2, 3, 4, 5])
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct InboundNatRulePropertiesFormat2 <: SwaggerModel
    frontendIPConfiguration::Any # spec type: Union{ Nothing, SubResource } # spec name: frontendIPConfiguration
    backendIPConfiguration::Any # spec type: Union{ Nothing, NetworkInterfaceIPConfiguration } # spec name: backendIPConfiguration
    protocol::Any # spec type: Union{ Nothing, TransportProtocol2 } # spec name: protocol
    frontendPort::Any # spec type: Union{ Nothing, Int32 } # spec name: frontendPort
    backendPort::Any # spec type: Union{ Nothing, Int32 } # spec name: backendPort
    idleTimeoutInMinutes::Any # spec type: Union{ Nothing, Int32 } # spec name: idleTimeoutInMinutes
    enableFloatingIP::Any # spec type: Union{ Nothing, Bool } # spec name: enableFloatingIP
    enableTcpReset::Any # spec type: Union{ Nothing, Bool } # spec name: enableTcpReset
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function InboundNatRulePropertiesFormat2(;frontendIPConfiguration=nothing, backendIPConfiguration=nothing, protocol=nothing, frontendPort=nothing, backendPort=nothing, idleTimeoutInMinutes=nothing, enableFloatingIP=nothing, enableTcpReset=nothing, provisioningState=nothing)
        o = new()
        validate_property(InboundNatRulePropertiesFormat2, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        setfield!(o, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("backendIPConfiguration"), backendIPConfiguration)
        setfield!(o, Symbol("backendIPConfiguration"), backendIPConfiguration)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("frontendPort"), frontendPort)
        setfield!(o, Symbol("frontendPort"), frontendPort)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("backendPort"), backendPort)
        setfield!(o, Symbol("backendPort"), backendPort)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        setfield!(o, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("enableFloatingIP"), enableFloatingIP)
        setfield!(o, Symbol("enableFloatingIP"), enableFloatingIP)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("enableTcpReset"), enableTcpReset)
        setfield!(o, Symbol("enableTcpReset"), enableTcpReset)
        validate_property(InboundNatRulePropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type InboundNatRulePropertiesFormat2

const _property_map_InboundNatRulePropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("frontendIPConfiguration")=>Symbol("frontendIPConfiguration"), Symbol("backendIPConfiguration")=>Symbol("backendIPConfiguration"), Symbol("protocol")=>Symbol("protocol"), Symbol("frontendPort")=>Symbol("frontendPort"), Symbol("backendPort")=>Symbol("backendPort"), Symbol("idleTimeoutInMinutes")=>Symbol("idleTimeoutInMinutes"), Symbol("enableFloatingIP")=>Symbol("enableFloatingIP"), Symbol("enableTcpReset")=>Symbol("enableTcpReset"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_InboundNatRulePropertiesFormat2 = Dict{Symbol,String}(Symbol("frontendIPConfiguration")=>"SubResource", Symbol("backendIPConfiguration")=>"NetworkInterfaceIPConfiguration", Symbol("protocol")=>"TransportProtocol2", Symbol("frontendPort")=>"Int32", Symbol("backendPort")=>"Int32", Symbol("idleTimeoutInMinutes")=>"Int32", Symbol("enableFloatingIP")=>"Bool", Symbol("enableTcpReset")=>"Bool", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ InboundNatRulePropertiesFormat2 }) = collect(keys(_property_map_InboundNatRulePropertiesFormat2))
Swagger.property_type(::Type{ InboundNatRulePropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_InboundNatRulePropertiesFormat2[name]))}
Swagger.field_name(::Type{ InboundNatRulePropertiesFormat2 }, property_name::Symbol) =  _property_map_InboundNatRulePropertiesFormat2[property_name]

function check_required(o::InboundNatRulePropertiesFormat2)
    true
end

function validate_property(::Type{ InboundNatRulePropertiesFormat2 }, name::Symbol, val)
end


text:
module TensorIntegration
    export make_tensor, tensor_integrate

    function simpson(a,b,N)
        if mod(N,2) == 0
            N += 1
        end
        if N > 1
            xs = collect(range(a, b, length=N))
            ws = ones(N)
            ws[2:2:end-1] .= 4.0
            ws[3:2:end-2] .= 2.0
            ws = ws*(xs[2]-xs[1])/3.
        else
            xs = [(a+b)/2.]
            ws = [1.]
        end
        return xs, ws
    end

    function tensor_simpson(a,b,N; lowprecision = false)
        if lowprecision
            numbertype = Float32
        else
            numbertype = Float64
        end
        d = length(N)
        for i = 1:d
            if mod(N[i],2) == 0
                N[i] += 1
            end
        end

        xs = [zeros(numbertype, N[d]) for i = 1:d]
        for i in 1:d
            if N[i] > 1
                xs[i] = collect(range(a[i], b[i], length = N[i]))
            else
                xs[i] = [(a[i]+b[i])/2.]
            end
        end
        ws = [zeros(numbertype, N[d]) for i = 1:d]
        for i in 1:d
            if N[i] > 1
                ws[i] = ones(N[i])
                ws[i][2:2:end-1] .= 4.0
                ws[i][3:2:end-2] .= 2.0
                ws[i] .= ws[i].*(xs[i][2]-xs[i][1])./3.
            else
                ws[i] = [1.]
            end
        end
        xmat, w = tensor_integrate(xs, ws, N)
        return xmat, w
    end

    function tensor_integrate(xs, ws, N)
        totN = prod(N)
        d = length(N)
        xmat = zeros(eltype(xs[1]), totN, d)
        wmat = zeros(eltype(xs[1]), totN, d)
        for i in 1:d
            repout = Int(totN/prod(N[1:i]))
            repin = Int(totN/N[i]/repout)
            xmat[:, i] .= repeat(xs[i], inner=repin, outer=repout)
            wmat[:, i] .= repeat(ws[i], inner=repin, outer=repout)
        end
        w = prod(wmat, dims = 2)
        return xmat, w
    end

    function make_tensor(xs, N)
        totN = prod(N)
        d = length(N)
        xmat = zeros(eltype(xs[1]), totN, d)
        for i in 1:d
            repout = Int(totN/prod(N[1:i]))
            repin = Int(totN/N[i]/repout)
            xmat[:, i] .= repeat(xs[i], inner=repin, outer=repout)
        end
        return xmat
    end
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	mga(EP::Model, path::AbstractString, setup::Dict, inputs::Dict, outpath::AbstractString)

We have implemented an updated Modeling to Generate Alternatives (MGA) Algorithm proposed by [Evelina et al., (2017)](https://www.sciencedirect.com/science/article/pii/S0360544217304097) to generate a set of feasible, near cost-optimal technology portfolios. This algorithm was developed by [Brill Jr, E. D., 1979](https://pubsonline.informs.org/doi/abs/10.1287/mnsc.25.5.413) and introduced to energy system planning by [DeCarolia, J. F., 2011](https://www.sciencedirect.com/science/article/pii/S0140988310000721).

To create the MGA formulation, we replace the cost-minimizing objective function of GenX with a new objective function that creates multiple generation portfolios by zone. We further add a new budget constraint based on the optimal objective function value $f^*$ of the least-cost model and the user-specified value of slack $\delta$. After adding the slack constraint, the resulting MGA formulation is given as:

```math
\begin{aligned}
	\text{max/min} \quad
	&\sum_{z \in \mathcal{Z}}\sum_{r \in \mathcal{R}} \beta_{z,r}^{k}P_{z,r}\\
	\text{s.t.} \quad
	&P_{zr} = \sum_{y \in \mathcal{G}}\sum_{t \in \mathcal{T}} \omega_{t} \Theta_{y,t,z,r}  \\
	& f \leq f^* + \delta \\
	&Ax = b
\end{aligned}
```

where, $\beta_{zr}$ is a random objective fucntion coefficient betwen $[0,100]$ for MGA iteration $k$. $\Theta_{y,t,z,r}$ is a generation of technology $y$ in zone $z$ in time period $t$ that belongs to a resource type $r$. We aggregate $\Theta_{y,t,z,r}$ into a new variable $P_{z,r}$ that represents total generation from technology type $r$ in a zone $z$. In the second constraint above, $\delta$ denote the increase in budget from the least-cost solution and $f$ represents the expression for the total system cost. The constraint $Ax = b$ represents all other constraints in the power system model. We then solve the formulation with minimization and maximization objective function to explore near optimal solution space.
"""
function mga(EP::Model, path::AbstractString, setup::Dict, inputs::Dict, outpath::AbstractString)

	println("Starting Model to Generate Alternatives (MGA) Iterations")

	if setup["ModelingToGenerateAlternatives"]==1
        # Start MGA Algorithm
	    println("MGA Module")

	    # Objective function value of the least cost problem
	    Least_System_Cost = objective_value(EP)

	    # Read sets
	    dfGen = inputs["dfGen"]
	    T = inputs["T"]     # Number of time steps (hours)
	    Z = inputs["Z"]     # Number of zonests
	    G = inputs["G"]

	    # Create a set of unique technology types
	    TechTypes = unique(dfGen[dfGen[!, :MGA] .== 1, :Resource_Type])

	    # Read slack parameter representing desired increase in budget from the least cost solution
	    slack = setup["ModelingtoGenerateAlternativeSlack"]

	    ### Variables ###

	    @variable(EP, vSumvP[TechTypes = 1:length(TechTypes), z = 1:Z] >= 0) # Variable denoting total generation from eligible technology of a given type

	    ### End Variables ###

	    ### Constraints ###

	    # Constraint to set budget for MGA iterations
	    @constraint(EP, budget, EP[:eObj] <= Least_System_Cost * (1 + slack) )

        # Constraint to compute total generation in each zone from a given Technology Type
	    @constraint(EP,cGeneration[tt = 1:length(TechTypes), z = 1:Z], vSumvP[tt,z] == sum(EP[:vP][y,t] * inputs["omega"][t]
	    for y in dfGen[(dfGen[!,:Resource_Type] .== TechTypes[tt]) .& (dfGen[!,:Zone] .== z),:][!,:R_ID], t in 1:T))

	    ### End Constraints ###

	    ### Create Results Directory for MGA iterations
        outpath_max = joinpath(path, "MGAResults_max")
	    if !(isdir(outpath_max))
	    	mkdir(outpath_max)
	    end
        outpath_min = joinpath(path, "MGAResults_min")
	    if !(isdir(outpath_min))
	    	mkdir(outpath_min)
	    end

	    ### Begin MGA iterations for maximization and minimization objective ###
	    mga_start_time = time()

	    print("Starting the first MGA iteration")

	    for i in 1:setup["ModelingToGenerateAlternativeIterations"]

	    	# Create random coefficients for the generators that we want to include in the MGA run for the given budget
	    	pRand = rand(length(unique(dfGen[dfGen[!, :MGA] .== 1, :Resource_Type])),length(unique(dfGen[!,:Zone])))

	    	### Maximization objective
	    	@objective(EP, Max, sum(pRand[tt,z] * vSumvP[tt,z] for tt in 1:length(TechTypes), z in 1:Z ))

	    	# Solve Model Iteration
	    	status = optimize!(EP)

            # Create path for saving MGA iterations
	    	mgaoutpath_max = joinpath(outpath_max, string("MGA", "_", slack,"_", i))

	    	# Write results
	    	write_outputs(EP, mgaoutpath_max, setup, inputs)

	    	### Minimization objective
	    	@objective(EP, Min, sum(pRand[tt,z] * vSumvP[tt,z] for tt in 1:length(TechTypes), z in 1:Z ))

	    	# Solve Model Iteration
	    	status = optimize!(EP)

            # Create path for saving MGA iterations
	    	mgaoutpath_min = joinpath(outpath_min, string("MGA", "_", slack,"_", i))

	    	# Write results
	    	write_outputs(EP, mgaoutpath_min, setup, inputs)

	    end

	    total_time = time() - mga_start_time
	    ### End MGA Iterations ###
	end

end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
# Simple example of an eigenproblem:
#
#           y
#             |
#             |            u_y = 0
#          Ly ----------------------------------
#             |                                |
#             |                                |
#             |                                |
#       u = 0 |     u_xx + u_yy + λ u = 0      | u = 0
#             |                                |
#             |                                |
#             |                                |
#           -----------------------------------|-----  x
#                          u_y = 0              Lx
#
# Here, exact solutions are
#
#      u = sin( n π x / Lx ) cos( k π y / Ly )
#
#                      2               2
#      λ = ( n π / Lx )  + ( k π / Ly )
#
# for n = 1, 2, 3, ... and k = 0, 1, 2, ....

using FinElt
using FinElt.PlanarPoisson
using Printf
using Arpack

include("params.jl")
const nev = 4

λ = Float64[]
for n = 1:nev, k = 0:nev-1
    push!(λ, (n*pi/Lx)^2+(k*pi/Ly)^2)
end

sort!(λ)

essential_bc = [ "Left", "Right" ]

err = zeros(nev, refinements+1)
@printf("Eigenvalue errors\n\n")
@printf("%7s|", "N")
for j = 1:nev
    @printf("%16s|", "λ_$j    ")
end
@printf(" elapsed\n")
@printf("%84s\n", "-"^84)
for k = 0:refinements
    start = time()
    mesh = read_msh_file("rect$k.msh")
    ep = EigenProblem(mesh, essential_bc)
    add_bilin_form!(ep, "Omega", grad_dot_grad!,   :LHS)
    add_bilin_form!(ep, "Omega", func_times_func!, :RHS)
    A, B = assembled_eigenproblem_matrices(ep)
    d, nconv, niter, nmult, resid = eigs(A, B, nev=nev, which=:SM,
                                         ritzvec=false)
    finish = time()
    err[:,k+1] = abs.(d-λ[1:nev])
    N = size(A, 1)
    if k == 0
        elapsed = finish - start
        @printf("%7d|", N)
        for j = 1:nev
            @printf("%9.2e %6s|", err[j,k+1], "")
        end
        @printf(" %7.4f\n", elapsed)
    else
        elapsed = finish - start
        @printf("%7d|", N)
        for j = 1:nev
            rate = log2(err[j,k]/err[j,k+1])
            @printf("%9.2e %6.4f|", err[j,k+1], rate)
        end
        @printf(" %7.4f\n", elapsed)
    end
end


text:
module ImageStitch

using FileIO, Images
#using OEMTIFF #zero array
using TiffImages
using ImageView
using NRRD

using RegisterQD, CoordinateTransformations
using OffsetArrays
using LinearAlgebra
using StatsBase
using JLD2
using Statistics #for median
using Unitful

export checkimgs, merged_axes, merge_images, getOffsetOrigin
# Write your package code here.
#import Unitful: s

## Check image
function checkimgs(img1, img2)
  imgsub0 = similar(img1)
    for i in 1:size(img1, 3)
      if i%2 != 0
        imgsub0[:,:,i] = img2[:,:,i]
      else
        imgsub0[:,:,i] = img1[:,:,i-1]
      end
    end
  return(imgsub0)
end

## Merge image and save data into the hard drive
function merged_axes(img1, imgw)
  axs =
      min(minimum(axes(imgw, 1)), minimum(axes(img1, 1))):max(maximum(axes(imgw, 1)), maximum(axes(img1, 1))),
      min(minimum(axes(imgw, 2)), minimum(axes(img1, 2))):max(maximum(axes(imgw, 2)), maximum(axes(img1, 2))),
      min(minimum(axes(imgw, 3)), minimum(axes(img1, 3))):max(maximum(axes(imgw, 3)), maximum(axes(img1, 3)))
  axs
end

function merge_images(fn, img1, imgw)
  axs = merged_axes(img1, imgw)
  ss = open(fn, "w+")
  for i in axs[3]
    for j in axs[2]
      for k in axs[1]
        if i in axes(imgw,3) && j in axes(imgw,2) && k in axes(imgw,1) && imgw[k,j,i] != 0
          write(ss, UInt16(imgw[k,j,i].*typemax(UInt16)))
        elseif i in axes(img1,3) && j in axes(img1,2) && k in axes(img1,1)
          write(ss, UInt16(img1[k,j,i].*typemax(UInt16)))
        else
          write(ss, zero(UInt16))
        end
      end
    end
  end
  close(ss)
end

## calculate the origin for image stitch
function getOffsetOrigin(img, overlap_coords)
  xOrigin = -(length(axes(img,1))-length(overlap_coords[1])-1)
  yOrigin = -(length(axes(img,2))-length(overlap_coords[2])-1)
  zOrigin = -(length(axes(img,3))-length(overlap_coords[3])-1)
  return(xOrigin, yOrigin, zOrigin)
end

end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ acc1f923-ff96-4edb-9c4e-13ed4a55f871
begin
	using Pkg
	Pkg.activate("..")
	using RiceBPH.PostProcess
end

# ╔═╡ f329fe9e-cb9b-11eb-1e8b-e7b224d1d323
using JLD2, PlutoUI, CairoMakie, DataFrames, CSV, Statistics

# ╔═╡ 8a9e6ee0-967b-4366-ba93-387a912c2e42
using Clustering

# ╔═╡ 48a326b4-cac9-43a8-b0e6-b2dfadcc8be3
using DataFramesMeta

# ╔═╡ 44833546-b966-430e-a4d1-2a51150f1e8b
using Latexify

# ╔═╡ f64c2ce4-2f5a-4db9-94b4-498620d8886c
files = joinpath.("../results", readdir("../results"))

# ╔═╡ 07b2fed7-dd09-4848-ad76-24aa305247ed
function ma(X, k)
	pad = zeros(eltype(X), k ÷ 2)
	[pad; [mean(X[i:(i+k)]) for i in 1:length(X) - k]; pad]
end

# ╔═╡ 85c54663-ba49-4b4c-bb3d-aff8f7f09636
function peak_population(X; smooth = 48 * 7 ÷ 2, threshold=0.0)
	# Smooth signal
	Y = X
	Y = ma(X, smooth)
	# Normlize
	MX = maximum(Y)
	mX = minimum(Y)
	Y = (Y .- mean(Y)) ./ std(Y)

	# Find the peaks
	peaks = let r = Y .≥ threshold
		ranges = findall(isone, abs.(diff(r)))
		if isodd(length(ranges))
			push!(range, length(X))
		end
		map(Iterators.partition(ranges, 2)) do (a, b)
			_, offset = findmax(@view X[a:b])
			a + offset
		end
	end
end

# ╔═╡ d1c54595-d257-4115-91f2-d9802fb2053a
function file2peaks(file)
	peaks = jldopen(file) do f
		map(1:1000) do seed
			df = f["$seed"]
			peak_population(df.count_bph)
		end
	end
	peaks = let lens = length.(peaks)
		s = std(lens)
		m = mean(lens)
		a, b = m - 3s, m + 3s
		filter(x -> a ≤ length(x) ≤ b, peaks)
	end
	npeaks = maximum(length.(peaks))
	peaks = reduce(vcat, peaks)

	cl = kmeans(transpose(peaks), npeaks)
	peaks = map(1:npeaks) do k
		peaks[k .== cl.assignments]
	end
	mpeaks = mean.(peaks)
	perm = sortperm(mpeaks)
	peaks[perm]
end

# ╔═╡ 0356fe38-14ef-4470-9506-bbdec875c222
function parsefile(path)
	meta = jldopen(f -> f["metadata"], path)
	name = last(split(meta.envmap, r"[\\/]"))
	(meta..., envmap=name, path=path)
end

# ╔═╡ 4cd2a9b0-c599-4b50-a96b-e9a46e32d7bf
function splitpeaks(df)
	function getstat(stat, x, n)
		stat(get(x, n, [missing]))
	end
	npeaks = maximum(length.(df.peaks))
	transformation = map(1:npeaks) do n
		[:peaks => ByRow(x -> getstat(mean, x, n)) => Symbol("peak_$n"),
		:peaks => ByRow(x -> getstat(std, x, n)) => Symbol("std_peak_$n")]
	end
	df = transform(df, transformation...)
	df = select(df, Not([:peaks, :path]))
end

# ╔═╡ e09b9532-1872-4ae7-8472-f48e668b9806
peakdf = let df = DataFrame(parsefile.(files))
	df
	df = transform(df, :path => ByRow(file2peaks) => :peaks)
end

# ╔═╡ 4c735053-adb2-4b3e-b3a6-e518e54c0e30
statpeak = splitpeaks(peakdf)

# ╔═╡ 553e76b2-349a-4719-9dac-792803534c84
cleanpeak = let to_day = x -> @. round(x / 24, digits=1)
	cols = names(statpeak)
	cols = filter(x -> occursin("peak", x), names(statpeak))
	tfm = map(cols) do col
		col => to_day => col
	end
	transform(statpeak, tfm...)
end

# ╔═╡ 3cbc8ee8-d2e8-41c6-a088-e273119aad26
cleandf2 = let df = cleanpeak
	df = @where(df, (:pr_killed .== 0.15) .| (:pr_killed  .== 0))
	df = select(df,  [[3, 1, 4, 2]; 5:size(df, 2)])
	df = transform(df, :envmap => ByRow(function(x)
				if x == "no-flower.csv"
					[0, missing]
				else
					m = match(r"(\d+)-(\d+x\d+).csv", x)
					if isnothing(m)
						[0, x]
					else
						[parse(Int, m[1]), m[2]]
					end
				end
			end) => [:F, :envmap])
	df = sort(df, [:F, :envmap, :pr_killed, :init_position, :nb_bph_init], rev=true)
	df = let cols = string.([:F, :envmap, :pr_killed, :init_position, :nb_bph_init])
		dcols = setdiff(names(df), cols)
		select(df, [cols; dcols])
	end
	cols = names(df)
	cols = map(enumerate(cols)) do (i, name)
		if i == 1
			raw"$F$"
		elseif i == 2
			raw"$L$"
		elseif i == 3
			raw"$\Pr_{\text{kill}}$"
		elseif i == 4
			raw"$P_0$"
		elseif i == 5
			raw"$N$"
		else
			name = replace(raw"$" * name * raw"$", "std" => raw"\sigma")
			name = replace(name, r"sigma_peak_(\d+)" => s"sigma(t_\1)")
			name = replace(name, raw"$sigma" => raw"$\sigma")
			name = replace(name, r"peak_(\d+)" => s"mu(t_\1)")
			name = replace(name, "mu(t_" => raw"\mu(t_")
		end
	end
	df = rename(df, cols)
	n = size(df, 2)
	# df = select(df, Not([n, n-1]))
	df
end

# ╔═╡ 8c263acf-8585-462a-a6b3-617161c6c3eb
function formaltable(df)
	s = latextabular(df; latex=false)
	lines = split(s, "\n")
	tbegin = replace(lines[1], "" => "")
	# header
	thead = strip(replace(lines[2], "\\\\" => ""))
	# thead = join("\\textbf{" .* split(thead, r" *& *") .* "}", " & ")
	tbody = join(lines[3:end-2], "\n\\midrule\n")
	s = """
	$(tbegin)
	\\toprule
	$(thead) \\\\
	\\midrule
	\\midrule
	$(tbody)
	\\bottomrule
	\\end{tabular}
	"""
	s = replace(s, r"([^a-zA-Z_.])([0-9]+\.?[0-9]*)([^a-zA-Z_.])" => s"\1$\2$\3")
	# s = replace(s, r"([0-9]+\.[0-9]+)" => s"$\1$")
	s = replace(s, "1x2" => raw"\OneByTwo")
	s = replace(s, "3x3" => raw"\ThreeByThree")
	s = replace(s, "false" => raw"\Reject")
	s = replace(s, "true" => raw"\Accept")
	s = replace(s, "corner" => raw"\Corner")
	s = replace(s, "border" => raw"\Border")
	# s = replace(s, r"_" => raw"\_")
	s = replace(s, "#" => raw"\#")
	s = replace(s, "missing" => "--")
	clipboard(s)
	return s
end

# ╔═╡ b1e1589a-f3fa-4756-af4a-308d268c006b
formaltable(cleandf2)

# ╔═╡ af6b92d3-1cd7-44c3-b417-2965406badaa
peakdfcsv = let df = CSV.read("all-peak-data.csv", DataFrame)
	transform(df, :peaks => ByRow(x -> eval(Meta.parse(x))) => :peaks)
end

# ╔═╡ 4ace47db-2992-4c0e-9076-9f37d4f52b6a
_statpeak = let df = splitpeaks(peakdfcsv)
	to_day = x -> @. round(x / 24, digits=1)
	cols = names(df)
	cols = filter(x -> occursin("peak", x), names(df))
	tfm = map(cols) do col
		col => to_day => col
	end
	transform(df, tfm...)
end

# ╔═╡ f7d70d28-e1f0-4369-8319-0c2661fd348e
finalpeak = let df = _statpeak
	df = @where(df, (:pr_killed .== 0.15) .| (:pr_killed  .== 0))
	df = select(df,  [[3, 1, 4, 2]; 5:size(df, 2)])
	df = transform(df, :envmap => ByRow(function(x)
				if x == "no-flower.csv"
					[0, missing]
				else
					m = match(r"(\d+)-(\d+x\d+).csv", x)
					if isnothing(m)
						[0, x]
					else
						[parse(Int, m[1]), m[2]]
					end
				end
			end) => [:F, :envmap])
	df = sort(df, [:F, :envmap, :pr_killed, :init_position, :nb_bph_init], rev=true)
	df = let cols = string.([:F, :envmap, :pr_killed, :init_position, :nb_bph_init])
		dcols = setdiff(names(df), cols)
		select(df, [cols; dcols])
	end
	cols = names(df)
	cols = map(enumerate(cols)) do (i, name)
		if i == 1
			raw"$F$"
		elseif i == 2
			raw"$L$"
		elseif i == 3
			raw"$\Pr_{\text{kill}}$"
		elseif i == 4
			raw"$P_0$"
		elseif i == 5
			raw"$N$"
		else
			name = replace(raw"$" * name * raw"$", "std" => raw"\sigma")
			name = replace(name, r"sigma_peak_(\d+)" => s"sigma(t_\1)")
			name = replace(name, raw"$sigma" => raw"$\sigma")
			name = replace(name, r"peak_(\d+)" => s"mu(t_\1)")
			name = replace(name, "mu(t_" => raw"\mu(t_")
		end
	end
	df = rename(df, cols)
	n = size(df, 2)
	# df = select(df, Not([n, n-1]))
	df
end

# ╔═╡ 7d9f469e-c857-4589-8339-da4c0e79d597
peakto_report = let n = size(finalpeak, 2)
	df = select(finalpeak, [1:5; 6+2*3:n])
	# df = select(finalpeak, 1:n-2*3)
	df = DataFrame(filter(eachrow(df)) do r
			!all(ismissing, r[end-5:end])
	end)

	formaltable(df)
end

# ╔═╡ 5e463e13-2215-42f5-8112-6987d202e65d
testdf = PostProcess.batch_test_rice("../results", 0.01; alpha=0.05, clean=true)

# ╔═╡ 9fd36c91-e299-4662-83ef-93bb8f290692
newnames=[
	:pr_killed => raw"$\Pr_{\text{kill}}$",
	:init_position => raw"$P_0$",
	:layout => raw"$L$",
	:nb_bph_init => raw"$N$",
]

# ╔═╡ 05e937cd-9395-4bb8-bbf7-de59b3f8e04e
function clean(df; ditch, group, key, val, newnames = nothing, sortcols=nothing)
	df = select(df, Not(ditch))
	function parseenvmap(m)
		ma = match(r"(\d+)-(\dx\d).csv", m)
		parse(Int, ma[1]), ma[2]
	end
	df = transform(df, :envmap => ByRow(parseenvmap) => [:flower, :layout])
	df = select(df, Not(:envmap))
	df = [unstack(d, key, val) for d in groupby(df, group)]
	df = vcat(df...)

	if !isnothing(sortcols)
		df = sort(df, sortcols)
		sortcols = String.(sortcols)
		othercols = setdiff(names(df), sortcols)
		df = select(df, [sortcols; othercols])
	end
	if !isnothing(newnames)
		df = rename(df, newnames)
	end
	return df
end

# ╔═╡ 05f9321d-0272-4720-9e0c-1b76358a1b1c
sortcols = [:pr_killed, :layout, :init_position, :nb_bph_init]

# ╔═╡ ad29fa9e-175f-4168-acac-a74b5e7f11a0
function clean_for_k(df)
	clean(df;
		ditch=[:accept, :n],
		group=[:pr_killed, :layout],
		key=:flower,
		val=:k,
		newnames=newnames,
		sortcols = sortcols
	)
end

# ╔═╡ b4b9588e-5a68-43e6-83f4-9c937f57a445
function clean_for_h0(df)
	clean(df;
		ditch=[:k, :n],
		group=[:pr_killed, :layout],
		key=:flower,
		val=:accept,
		newnames=newnames,
		sortcols=sortcols,
	)
end

# ╔═╡ d8126628-ee50-4d73-afa1-6e6b3a355bc8
testdf2 = @where(testdf, (:envmap .!= "no-flower.csv") .& (:pr_killed .> 0));

# ╔═╡ e6a060d5-0ffc-4a71-bbe6-df8d7fb5b15a
formaltable(clean_for_k(testdf2))

# ╔═╡ Cell order:
# ╠═f329fe9e-cb9b-11eb-1e8b-e7b224d1d323
# ╠═f64c2ce4-2f5a-4db9-94b4-498620d8886c
# ╠═07b2fed7-dd09-4848-ad76-24aa305247ed
# ╠═85c54663-ba49-4b4c-bb3d-aff8f7f09636
# ╠═d1c54595-d257-4115-91f2-d9802fb2053a
# ╠═0356fe38-14ef-4470-9506-bbdec875c222
# ╠═4cd2a9b0-c599-4b50-a96b-e9a46e32d7bf
# ╠═8a9e6ee0-967b-4366-ba93-387a912c2e42
# ╠═e09b9532-1872-4ae7-8472-f48e668b9806
# ╠═4c735053-adb2-4b3e-b3a6-e518e54c0e30
# ╠═553e76b2-349a-4719-9dac-792803534c84
# ╠═48a326b4-cac9-43a8-b0e6-b2dfadcc8be3
# ╠═3cbc8ee8-d2e8-41c6-a088-e273119aad26
# ╠═44833546-b966-430e-a4d1-2a51150f1e8b
# ╠═b1e1589a-f3fa-4756-af4a-308d268c006b
# ╠═8c263acf-8585-462a-a6b3-617161c6c3eb
# ╠═af6b92d3-1cd7-44c3-b417-2965406badaa
# ╠═4ace47db-2992-4c0e-9076-9f37d4f52b6a
# ╟─f7d70d28-e1f0-4369-8319-0c2661fd348e
# ╠═7d9f469e-c857-4589-8339-da4c0e79d597
# ╠═acc1f923-ff96-4edb-9c4e-13ed4a55f871
# ╠═5e463e13-2215-42f5-8112-6987d202e65d
# ╠═9fd36c91-e299-4662-83ef-93bb8f290692
# ╠═05e937cd-9395-4bb8-bbf7-de59b3f8e04e
# ╠═05f9321d-0272-4720-9e0c-1b76358a1b1c
# ╠═ad29fa9e-175f-4168-acac-a74b5e7f11a0
# ╠═b4b9588e-5a68-43e6-83f4-9c937f57a445
# ╠═d8126628-ee50-4d73-afa1-6e6b3a355bc8
# ╠═e6a060d5-0ffc-4a71-bbe6-df8d7fb5b15a

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
function closeness(g::AbstractStaticGraph)

end

function closeness(g::AbstractEvolvingGraph)

end


text:
using Plots

x = -π:0.01:π
ys=hcat(sin.(x), cos.(x),tan.(x),sinh.(x))
plot(x,ys,layout=(4,1))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Statistics
using Printf
using Dates

Base.MainInclude.include("dataStruct.jl")
# Base.include(Main, "ConfigFileLib.jl")
Base.MainInclude.include("ConfigFileLib.jl")
const segNum =8192
const maxhistorydataCnt = 2*60 #*60

#######################################################################################################################
extensionName=""
headlen=512
cfglen=14
extensionsAll=( ".RAW1",  ".RAW3" ,".RAW4" ,".LC" ,".RAW" , )
extensions1=( ".RAW1",   ".RAW" , )
extensions2=( ".RAW2", ".FIL" , )
extensions3=(  ".RAW3" , ".LC"  , )

mutable struct DatetimeData
    dt::String
end

abstract type NULLDATA
end

###################################################################################
function readin!(rawdatafileObj,cfgdata::cfgStruct)
    global     extensionName
    # println( "readin  ",extensionName )
    if  extensionName  in ( ".RAW3"  ,".RAW4"  ,".RAW1",".RAW2")
            seek( rawdatafileObj,512)
            cfgdata.meterperSeg = read( rawdatafileObj, UInt16 )/5000
            cfgdata.chn2SegBegIdx = read( rawdatafileObj, UInt16 )
            cfgdata.reflectorFactor = read( rawdatafileObj, UInt16 )/10000
            cfgdata.attenutionFactor = read( rawdatafileObj, UInt16 )/10000
            cfgdata.scanrate = read( rawdatafileObj, UInt16 )
            cfgdata.calSamples = read( rawdatafileObj, UInt16 )
            cfgdata.SegNumber = read( rawdatafileObj, UInt16 )
    end
end

function readin(rawdatafileObj,NULLDATA)
    global     extensionName
    # println( "readin  ",extensionName )
    if  extensionName  in (   ".RAW3"  ,".RAW4"  ,".RAW1",".RAW2")
            seek( rawdatafileObj,512)
            skip( rawdatafileObj,14 )
    end
end

function readin!(rawdatafileObj,dtdata::DatetimeData)
    global extensionName

    if  extensionName  in (   ".RAW3"  ,".RAW4"  ,".RAW1",".RAW2")
                dtdata.dt=string("")
            # println("===",dtdata.dt,"------------" )
                curpos = position(rawdatafileObj )
                for i in 1:14
                    s =  read( rawdatafileObj,UInt8 )
                    dtdata.dt =string( dtdata.dt ,Char(s) )
                end
                s =  read( rawdatafileObj,UInt16)  # readin the type info
    else
        b1 =  read( rawdatafileObj,UInt8)  # readin the type info
        b2 =  read( rawdatafileObj,UInt8)  # readin the type info
        b3 =  read( rawdatafileObj,UInt8)  # readin the type info
        b4 =  read( rawdatafileObj,UInt8)  # readin the type info
        dtdata.dt =string("20201012123456" )

        df = Dates.DateFormat("y-m-d");
        jzdt = Dates.DateTime( "2010-1-1" ,df)
        a =  UInt(1)::UInt64
        a= b1*1 +b2*256 +b3*256*256 +b4*256*256*256
        jzdt1= jzdt + Dates.Second( a )
        # print( jzdt1," -----" )
        dtdata.dt =@sprintf( "%04d%02d%02d%02d%02d%02d",Dates.year(jzdt1 ) ,   Dates.month(jzdt1 ), Dates.day(jzdt1 ), Dates.hour(jzdt1 ),Dates. minute(jzdt1 ), Dates.second(jzdt1 ) )
    end
end

function readin!(rawdatafileObj,data::Array{UInt16,1})
    segNumber= length(data )
    # curpos = position(rawdatafileObj )
    for i in 1:segNumber
        s =  read( rawdatafileObj, UInt16 )
        data[i]=s
    end

end

## THE SPEED SLOW THEN ABOVE CODE
function readin1!(rawdatafileObj,data::Array{UInt16,1})
    segNumber= length(data )

    gr = reinterpret(UInt8,data)
    readbytes!( rawdatafileObj,  gr )

end

function show(dtdata::DatetimeData)
    print( dtdata.dt[1:4] , "-",dtdata.dt[5:6], "-", dtdata.dt[7:8] ,"  " ,dtdata.dt[9:10]  ,":", dtdata.dt[11:12] , ":",dtdata.dt[13:14]," ," )
end

function showRawData( ID , Dtdata::DatetimeData,RowData::Array{UInt16,1} )
    print( "ID :$ID  " )
    show(Dtdata   )
    println("data len is :$(length(RowData)  )"  )
end

function  writeTitle(CumSumDataobj)
    meter = CumSumDataobj.CfgInfoObj.meterperSeg
    chn2Pos = CumSumDataobj.CfgInfoObj.chn2SegBegIdx
    ## write the title s: 通道号下的光缆 时间、各段的米标 。
    print( "时间" ,"," )
    for segID = CumSumDataobj.AlrmProcStartID : CumSumDataobj.AlrmProcEndtID
        if segID  >= chn2Pos;
            print( round( UInt,Float64(segID-chn2Pos+1)* meter) ,"  $segID 米," )
        else
            print( round( UInt,Float64( segID  )* meter )," $segID 米," )
        end
    end
    println(   )
end

function saveLogInfo(spaceCumSumDataObj ,row  )
f=open("tem001.data","a")
for i in  row;
    print(f,i,",")
end
println(f)
close(f)
end

function readinDatainFile( filename,spaceCumSumDataObj ,Count=-1)
    global extensionName
    fn = basename( filename )
    extensionName  =uppercase( splitext( fn)[2]  )
    # println(  abspath( filename ) ,"   ",fn,"   " ,extensionName ,"   ",splitext( fn)[1]   )
    if  extensionName == ".RAW1"  ||   extensionName  == ".RAW"
        SegNumber =8192*512
    elseif  extensionName  in (   ".RAW3"  ,".RAW4"  ,".LC")
        SegNumber =8192*1
    end

    rawdatafileObj= open( filename,"r")
    # cfgdata =  cfgStruct()
    # readin!(rawdatafileObj,cfgdata)
    readin(rawdatafileObj,NULLDATA)
    # historyNum =3
    cnt=1
    dtdata = DatetimeData("")
    row = Array{UInt16,1}(undef,SegNumber)

    while( !eof(rawdatafileObj ) )
            readin!(rawdatafileObj,dtdata)
            readin!(rawdatafileObj,row)
            setDTValue( spaceCumSumDataObj , dtdata.dt)
            #@time
            # saveLogInfo(spaceCumSumDataObj ,row  )
            postProcFeatureData(spaceCumSumDataObj ,row )
            cnt= cnt+1
            if cnt > Count && Count >= 0
                break
            end
        end
    close( rawdatafileObj )
end

function getAllfiles( rootfolder ,extensions=( ".RAW1",  ".RAW3" ,".RAW4" ,".LC" ,".RAW" , ) )
    fns= []
    # rootfolder=cwd
    for  (root, dirs, files) in  walkdir(rootfolder)
        for fn  in files
               extensionName  =uppercase( splitext( fn)[2]  )
                if  extensionName in extensions
                    fullname = joinpath(root, fn )
                    push!(fns ,  fullname)
                    # println(fullname) # path to files
                end
        end
    end
    return fns
end #getAllfiles

function entryFun( fns  ,spaceCumSumDataObj,Count =-1 )
    for  file in fns
        println( "\n\nFilename: $file")
        readinDatainFile( file,spaceCumSumDataObj,Count )
        println()
    end
end

############################
function calInternal(chnid,startMeters ,endmeters,maxTimeCnt1, dataFilename ,cfgfilename ="", ProtectMeters =50)
    cfgObj = cfgStruct()
    aRNG = RNG()
    CalObj = spaceCumSumData( )
    if  cfgfilename != ""  ;
        read!(cfgObj,cfgfilename)
    else
        readin!(cfgObj,dataFilename)
    end
    show( cfgObj  )
    getRNG!( aRNG,cfgObj, startMeters ,endmeters ,ProtectMeters,chnid )
    show( aRNG)
    CalObj.CfgInfoObj = cfgObj
    InitObj(CalObj, aRNG.AlrmProcStartID,aRNG.AlrmProcEndtID1 , aRNG.startID1 ,aRNG.endID1,maxTimeCnt1 )

    if typeof(dataFilename) == String
        tfiles=[dataFilename, ]
    else
        tfiles = dataFilename
    end
    writeTitle(CalObj)
    entryFun( tfiles  ,CalObj,maxTimeCnt1 )
end

function calPT(chnid,startMeters ,endmeters,maxTimeCnt1, dataFilename , ProtectMeters =50)
    cfgfilename =""
    return calInternal(chnid,startMeters ,endmeters,maxTimeCnt1, dataFilename ,cfgfilename  , ProtectMeters  )
end

function calPTFiles(chnid,startMeters ,endmeters,maxTimeCnt1, dataFilenames , ProtectMeters =50)
    if typeof(dataFilenames) == String
        tfiles=[dataFilenames, ]
    else
        tfiles = dataFilenames
    end

    if length(tfiles ) <1
        return -1;
    end
        dataFilename = tfiles[1]
        cfgObj = cfgStruct()
        aRNG = RNG()
        CalObj = spaceCumSumData( )
        readin!(cfgObj,dataFilename)
        show( cfgObj  )
        getRNG!( aRNG,cfgObj, startMeters ,endmeters ,ProtectMeters,chnid )
        show( aRNG)
        CalObj.CfgInfoObj = cfgObj
        InitObj(CalObj, aRNG.AlrmProcStartID,aRNG.AlrmProcEndtID1 , aRNG.startID1 ,aRNG.endID1,maxTimeCnt1 )
        writeTitle(CalObj)
        # entryFun( dataFilename  ,CalObj,maxTimeCnt1 )
        for dataFilename in tfiles
            # println( "\n\nFilename: $dataFilename")
            readinDatainFile( dataFilename,CalObj,maxTimeCnt1 )
            # println()
        end
        return 0
end

function GetMachinenameFromFilename(filename)

    # tgrname :Data^Ver00^GYYS-001^20190304090024^ID00.*
    a= split(filename,"^")
    if length(a ) >=5
        return a[3]
    else
        return "OFFline-0001"
    end
end

function readinPT_1_File( filename,cfgObj::cfgStruct )
    global extensionName
    # @warn "Will Put Data into Channel....."
    fn = basename( filename )
    extensionName  =uppercase( splitext( fn)[2]  )
    println(  abspath( filename ) ,"   ",fn,"   " ,extensionName ,"   ",splitext( fn)[1]   )
    if  extensionName == ".RAW1"  ||   extensionName  == ".RAW"
        SegNumber =8192*512
    elseif  extensionName  in (   ".RAW3"  ,".RAW4"  ,".LC")
        SegNumber =8192*1

    end
    if extensionName == ".RAW3"
        rawDatatype =0x0003
    elseif extensionName == ".RAW1"
        rawDatatype =0x0001
    elseif extensionName == ".RAW4"
        rawDatatype =0x0004
    elseif extensionName == ".RAW2"
        rawDatatype =0x0002
    end
    rawdatafileObj= open( filename,"r")
    readin(rawdatafileObj,NULLDATA)
    # println("1")
    cfgArray = Array{UInt16,1}(undef,7)
    cfgArray[1] = floor(UInt16,cfgObj.meterperSeg*5000)
    cfgArray[2] = cfgObj.chn2SegBegIdx
    cfgArray[3] = 0 #reflectorFactor::Float64
    cfgArray[4] = 0 #attenutionFactor::Float64
    cfgArray[5] = cfgObj.scanrate
    cfgArray[6] = cfgObj.calSamples
    cfgArray[7] = cfgObj.SegNumber

    show( cfgObj )
    # return
    cnt=1
    dtdata = DatetimeData("")
    machineName=GetMachinenameFromFilename(fn)
    while( !eof( rawdatafileObj ) )
        # println("1")
        aobj = fileFiberRawData(SegNumber)
        dtstr = ""
        row = Array{UInt16,1}(undef,SegNumber)
            readin!(rawdatafileObj,dtdata)
            for i =1:length( dtdata.dt)
                aobj.startDTStr[i] = dtdata.dt[i]
                dtstr  = string(dtstr,Char(dtdata.dt[i]) )
            end

            if Main.startTime != "" &&  Main.endTime != ""
                # @info dtstr, Main.startTime , Main.endTime
                if  dtstr >=  Main.startTime &&  dtstr <= Main.endTime
                    ;
                else
                    # @info  "Not in time RNG "
                    skip(rawdatafileObj, length( row) *2 )
                    continue
                    # close( rawdatafileObj )
                    # return
                end
            end
            readin!(rawdatafileObj,row)
            aobj.datas .= row
            aobj.rawDatatype =  rawDatatype
            for i =1:length( machineName)
                aobj.machineName[i] = machineName[i]
            end
            len = length( machineName ) +1
            aobj.machineName[len] = 0x00

            aobj.cfgdata .= cfgArray
            # print("+")
            # global fileRawDataChannel
            # @warn "Will Put Data into Channel....."
            put!( Main.fileRawDataChannel ,aobj )
            # println(".+")
    end
    close( rawdatafileObj )
end

function fetchOfflineFiberDatas(fdobj)

    aobj = take!(  Main.fileRawDataChannel   )

    fdobj.rawDatatype[] = aobj.rawDatatype

    fdobj.machineName .= aobj.machineName
    fdobj.startDTStr .= aobj.startDTStr
    fdobj.cfgdata .= aobj.cfgdata
    # println( getStringFromArray( fdobj.machineName  ) )
    len =aobj.actualData
    @views a = fdobj.datas[1:len]
    a .= aobj.datas

    actualDatabytes = length( aobj.datas )*2
    fdobj.actualDatabytes = Ref{Cint}(actualDatabytes)
    # exit(-1)
end

function readinPTFiles( dataFilenames)
    # @warn "po？p"
        # println( dataFilenames )
        if typeof(dataFilenames) == String
            tfiles=[dataFilenames, ]
        else
            tfiles = dataFilenames
        end

        if length(tfiles ) <1
            println("No files found ,Exit out !!")
            exit(0)
            return -1;
        end

        for dataFilename in tfiles
            cfgObj = cfgStruct()
            # println( dataFilename )
            readin!(cfgObj,dataFilename)
            # println( dataFilename )

            # t = calPeriodLen(dataFilename )

            readinPT_1_File( dataFilename ,cfgObj)
            # println( dataFilename )
        end
        @warn "------------will exit---------cur time is :", Dates.now()
        println( "All File is Processed . will Exit out After 10 Seconds ...." )
        sleep(10)
        exit(0)
        return 0
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
#------------------------------
#DEFINE STATIONARY LISTENER w/ OMNIDIRECTIONAL ANTENNA

struct stationaryListenerOmnidirectional
   sourceList::Vector{Source}
   position::Vector{Float64}
end

function (listener::stationaryListenerOmnidirectional)(t₀::Float64)
sourceList = listener.sourceList
𝐩ᵣ = listener.position
   val = 0.0
   for i = 1:length(sourceList)
      val+=sourceList[i](𝐩ᵣ,t₀)
   end
   return val
end

#------------------------------
#DEFINE STATIONARY LISTENER w/ DIRECTIONAL ANTENNA and TIME-INVARIANT BEAM CENTER

struct stationaryListenerDirectional
   sourceList::Vector{Source}
   position::Vector{Float64}
   beamCenter::Vector{Float64}
   antennaGain ::Function
end

function (listener::stationaryListenerDirectional)(t₀::Float64)
   sourceList = listener.sourceList
   𝐩ᵣ = listener.position
   𝐛, G = listener.beamCenter , listener.antennaGain
      val = 0.0
      for i = 1:length(sourceList)
         if isa(sourceList[i],movingSource)
            val+=sourceList[i](𝐩ᵣ,t₀) * G( angleBetween(𝐛, 𝐩ᵣ-sourceList[i].position(t₀)) )
         else
            val+=sourceList[i](𝐩ᵣ,t₀) * G( angleBetween(𝐛, 𝐩ᵣ-sourceList[i].position) )
         end
      end
      return val
end

#------------------------------
#DEFINE MOVING LISTENER w/ OMNIDIRECTIONAL ANTENNA

struct movingListenerOmnidirectional
   sourceList::Vector{Source}
   position::Function
end

function (listener::movingListenerOmnidirectional)(t₀::Float64)
sourceList = listener.sourceList
𝐩ᵣ = listener.position
   val = 0.0
   for i = 1:length(sourceList)
      val+=sourceList[i](𝐩ᵣ(t₀),t₀)
   end
   return val
end

#------------------------------
#DEFINE MOVING LISTENER w/ DIRECTIONAL ANTENNA and TIME-VARYING BEAM CENTER

struct movingListenerDirectional
   sourceList::Vector{Source}
   position::Function
   beamCenter::Function
   antennaGain ::Function
end

function (listener::movingListenerDirectional)(t₀::Float64)
sourceList = listener.sourceList
𝐩ᵣ = listener.position
𝐛, G = listener.beamCenter , listener.antennaGain
   val = 0.0
   for i = 1:length(sourceList)
      if isa(sourceList[i],movingSource)
         val+=sourceList[i](𝐩ᵣ(t₀),t₀) * G( angleBetween(𝐛(t₀), 𝐩ᵣ(t₀)-sourceList[i].position(t₀)) )
      else
         val+=sourceList[i](𝐩ᵣ(t₀),t₀) * G( angleBetween(𝐛(t₀), 𝐩ᵣ(t₀)-sourceList[i].position) )
      end
   end
   return val
end

#------------------------------
#DEFINE GENERAL LISTENER
stationaryListener = Union{   stationaryListenerOmnidirectional,
                              stationaryListenerDirectional,
                              }

movingListener = Union{ movingListenerOmnidirectional,
                        movingListenerDirectional,
                        }

Listener = Union{ stationaryListener,
                  movingListener,
                  }

#multi-thread model evaluation over a time interval
function (z::Listener)(t::Vector{Float64})
   Z = zeros( typeof(z(0.0)), size(t))
   Threads.@threads for i = 1:length(t)
      Z[i] = z(t[i])
   end
   return Z
end

max_stars_repo_path: src/Basic LTI Sources/listeners.jl

max_issues_repo_path: src/Basic LTI Sources/listeners.jl

max_forks_repo_path: src/Basic LTI Sources/listeners.jl


text:
module Klusters
using ..Common, LightXML
export readklusters

"""
readklusters(path::AbstractString) -> GenericSpikeChannel

Read spike data from Klusters.
"""
function readklusters(path::AbstractString)
    dotindex = rsearch(path, '.')
    stem = path[1:rsearch(path, '.', dotindex-1)-1]
    chnumber = parse(Int, path[dotindex+1:end])

    xdoc = parse_file(stem*".xml")
    xroot = root(xdoc)
    acquisitionSystem = get_elements_by_tagname(xroot, "acquisitionSystem")[1]
    sample_rate = parse(Float64, content(get_elements_by_tagname(acquisitionSystem, "samplingRate")[1]))
    nbits = parse(Int, content(get_elements_by_tagname(acquisitionSystem, "nBits")[1]))
    voltage_multiplier = parse(Float64, content(get_elements_by_tagname(acquisitionSystem, "voltageRange")[1]))/2^nbits
    nbits == 16 || nbits == 32 || throw(ArgumentError("nBits == $nbits unsupported"))

    spikeDetection = get_elements_by_tagname(xroot, "spikeDetection")[1]
    groups = get_elements_by_tagname(get_elements_by_tagname(spikeDetection, "channelGroups")[1], "group")
    nsamples = 0
    for group in groups
        has_channel = false
        for channel in get_elements_by_tagname(get_elements_by_tagname(group, "channels")[1], "channel")
            if parse(Int, content(channel)) == chnumber
                has_channel = true
                break
            end
        end
        !has_channel && continue
        nsamples = parse(Int, content(get_elements_by_tagname(group, "nSamples")[1]))
    end
    nsamples == 0 && throw(ArgumentError("number of samples for channel $chnumber not specified in $stem.xml"))

    times = open("$stem.fet.$chnumber") do f
        readline(f)
        out = Float64[]
        for line in eachline(f)
            line == "\n" || line == "\r\n" && continue
            push!(out, parse(Int, line[rsearch(line, ' ')+1:end]))
        end
        out
    end
    times /= sample_rate
    if isfile("$stem.clu.$chnumber")
        unit_numbers = open("$stem.clu.$chnumber") do f
            readline(f)
            out = Int16[]
            for line in eachline(f)
                line == "\n" || line == "\r\n" && continue
                push!(out, parse(Int, line))
            end
            out
        end
    else
        unit_numbers = zeros(Int16, length(times))
    end
    data = open(f->read(f, nbits == 16 ? Int16 : Int32, (nsamples, div(filesize(f), nsamples*nbits÷8))), "$stem.spk.$chnumber")

    length(times) == length(unit_numbers) == size(data, 2) ||
        throw(ArgumentError("number of waveforms, times, and unit numbers do not match"))

    GenericSpikeChannel("ch$chnumber", sample_rate, voltage_multiplier, times, unit_numbers, data)
end
end


text:
# CuNNFTPlan2d
mutable struct CuNFFTPlan{D,T} <: AbstractNFFTPlan{D,0,T}
  N::NTuple{D,Int64}
  M::Int64
  x::Array{T,2}
  m::Int64
  sigma::T
  n::NTuple{D,Int64}
  K::Int64
  windowHatInvLUT::CuArray{Complex{T},D}
  forwardFFT::CUDA.CUFFT.cCuFFTPlan{Complex{T},-1,true,D}
  backwardFFT::CUDA.CUFFT.cCuFFTPlan{Complex{T},1,true,D}
  tmpVec::CuArray{Complex{T},D}
  tmpVec2::CuArray{Complex{T},D}
  apodIdx::CuArray{Int64,1}
  B::CuSparseMatrixCSC{Complex{T}}
end

# constructor for CuNFFTPlan2d
function CuNFFTPlan(x::Matrix{T}, N::NTuple{D,Int}, m=4, sigma=2.0,
  window=:kaiser_bessel, K=2000; kwargs...) where {T,D}

  n = ntuple(d->round(Int,sigma*N[d]), D)

  tmpVec = CuArray(zeros(Complex{T},n))
  tmpVec2 = CuArray(zeros(Complex{T},N))
  M = size(x,2)

  # create FFT-plans
  dims = ntuple(d->d, D)
  FP = plan_fft!(tmpVec, dims)
  BP = plan_bfft!(tmpVec, dims)

  # Create lookup table for 1d interpolation kernels
  win, win_hat = NFFT.getWindow(window)
  windowHatInvLUT  = precomp_windowHatInvLUT(T, win_hat, N, sigma, m)
  windowHatInvLUT_d = CuArray(windowHatInvLUT)

  # compute (linear) indices for apodization (mapping from regular to oversampled grid)
  apodIdx = precomp_apodIdx(N,n)
  apodIdx_d = CuArray(apodIdx)

  # precompute interpolation matrix
  U1 = ntuple(d-> (d==1) ? 1 : (2*m+1)^(d-1), D)
  U2 = ntuple(d-> (d==1) ? 1 : prod(n[1:(d-1)]), D)
  B = precomputeB(win, x, n, m, M, sigma, T, U1, U2)
  B_d = CuSparseMatrixCSC(Complex{T}.(B))

  return CuNFFTPlan{D,T}(N, M, x, m, sigma, n, K, windowHatInvLUT_d,
                  FP, BP, tmpVec, tmpVec2, apodIdx_d, B_d )
end

size(p::CuNFFTPlan) = p.N
numFourierSamples(p::CuNFFTPlan) = p.M

"""  in-place NFFT on the GPU"""
function nfft!(p::CuNFFTPlan{D,T}, f::CuArray, fHat::CuArray) where {D,T}
  consistencyCheck(p, f, fHat)

  # apodization
  fill!(p.tmpVec, zero(Complex{T}))
  p.tmpVec2 .= f.*p.windowHatInvLUT
  p.tmpVec[p.apodIdx] = p.tmpVec2

  # FFT
  p.forwardFFT * p.tmpVec # equivalent to fft!(p.tmpVec)

  # convolution
  mul!(fHat, transpose(p.B), vec(p.tmpVec))

  return fHat
end

"""  in-place adjoint NFFT on the GPU"""
function nfft_adjoint!(p::CuNFFTPlan{D,T}, fHat::CuArray, f::CuArray) where {D,T}
  consistencyCheck(p, f, fHat)

  # adjoint convolution
  mul!(vec(p.tmpVec), p.B, fHat)

  # FFT
  p.backwardFFT * p.tmpVec # bfft!(p.tmpVec)

  # adjoint apodization
  p.tmpVec2[:] = p.tmpVec[p.apodIdx]
  f .= p.tmpVec2.*p.windowHatInvLUT

  return f
end


text:
using ONS, Test

@test typeof(list_datasets()) == ONS.Records
@test typeof(list_timeseries()) == ONS.Records
@test typeof(search_datasets("RPI")) == ONS.Records
@test typeof(get_dataset("MM23")) == ONS.Records
@test typeof(get_data("MM23","CHAW")) == ONS.Data