text:
########################################################################
# Training
########################################################################

# Implementation of the FedAvg and FedProx algorithm
function fedAvgAndProx(
    server::Union{FedProxServer, FedProxServerNN},
    clients::Union{Vector{FedProxClient}, Vector{FedProxClientNN} },
    numRounds::Int64
)
    # Connect clients with server
    server.clients = clients
    # Training process
    objList = zeros(Float64, 0)
    testAccList = zeros(Float64, 0)
    startTime = time()
    @printf("Start training!\n")
    for t = 1:numRounds
        select!(server)
        sendModel!(server)
        # objValue = getObjValue(server)
        # @printf("Round : %4d, obj: %6.4e\n", t, objValue)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            update!(client)
        end
        @printf("aggregate\n")
        aggregate!(server)
        # Print log
        sendModelToAllClients!(server)
        @printf("getObj\n")
        objValue = getObjValue(server)
        acc = accuracy(server.Xtest, server.Ytest, server.W)
        @printf("Round : %4d, obj: %6.4e, acc: % 3.2f %%, time: %4.3f s\n", t, objValue, acc*100, time()-startTime)
        push!(objList, objValue)
        push!(testAccList, acc)
        if isnan(objValue)
            break
        end
    end
    endTime = time()
    @printf("Finished training, time elapsed: %.4e\n", endTime - startTime)
    return server.W, objList, testAccList
end

# Implementation of the Scaffold algorithm.
function Scaffold(
    server::Union{ScaffoldServer, ScaffoldServerNN},
    clients::Union{Vector{ScaffoldClient}, Vector{ScaffoldClientNN} },
    numRounds::Int64
)
    # Connect clients with server
    server.clients = clients
    # Training process
    objList = zeros(Float64, 0)
    testAccList = zeros(Float64, 0)
    startTime = time()
    @printf("Start training!\n")
    for t = 1:numRounds
        select!(server)
        sendModel!(server)
        # objValue = getObjValue(server)
        # @printf("Round : %4d, obj: %6.4e\n", t, objValue)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            update!(client)
        end
        aggregate!(server)
        # Print log
        sendModelToAllClients!(server)
        objValue = getObjValue(server)
        acc = accuracy(server.Xtest, server.Ytest, server.W)
        @printf("Round : %4d, obj: %6.4e, acc: % 3.2f %%, time: %4.3f s\n", t, objValue, acc*100, time()-startTime)
        push!(objList, objValue)
        push!(testAccList, acc)
        if isnan(objValue)
            break
        end
    end
    endTime = time()
    @printf("Finished training, time elapsed: %.4e\n", endTime - startTime)
    return server.W, objList, testAccList
end

# Implementation of the FedDCD algorithm (both exact and inexact)
function fedDCD(
    server::FedDCDServer,
    clients::Vector{FedDCDClient},
    numRounds::Int64
)
    # Connect clients with server
    server.clients = clients
    # Training process
    objList = zeros(Float64, 0)
    testAccList = zeros(Float64, 0)
    startTime = time()
    @printf("Start training!\n")
    for t = 1:numRounds
        # @printf("Round %d\n", t)
        select!(server)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            update!(client)
        end
        aggregate!(server)
        sendModel!(server)
        # Print log
        objValue = getObjValue(server)
        acc = accuracy(server.Xtest, server.Ytest, server.W)
        @printf("Round : %4d, obj: %6.4e, acc: % 3.2f %%, time: %4.3f s\n", t, objValue, acc*100, time()-startTime)
        push!(objList, objValue)
        push!(testAccList, acc)
        if isnan(objValue)
            break
        end
    end
    endTime = time()
    @printf("Finished training, time elapsed: %.4e\n", endTime - startTime)
    return server.W, objList, testAccList
end

function fedDCD(
    server::FedDCDServerNN,
    clients::Vector{FedDCDClientNN},
    numRounds::Int64
)
    # Connect clients with server
    server.clients = clients
    # Training process
    objList = zeros(Float64, 0)
    testAccList = zeros(Float64, 0)
    startTime = time()
    @printf("Start training!\n")
    for t = 1:numRounds
        # @printf("Round %d\n", t)
        select!(server)
        synchronize!(server)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            update!(client)
        end
        aggregate!(server)
        sendModel!(server)
        decay_lr!(server)
        # Print log
        objValue = getObjValue(server)
        # acc = 0.0
        # num_data = 0
        # for c in server.clients
        #     n = size(c.Xtrain, 1)
        #     num_data += n
        #     acc += n*accuracy(c.XtrainT, c.Ytrain, server.W)
        # end
        # acc /= num_data
        acc = accuracy(server.Xtest, server.Ytest, server.W)
        @printf("Round : %4d, obj: %6.4e, acc: % 3.2f %%, time: %4.3f s\n", t, objValue, acc*100, time()-startTime)
        push!(objList, objValue)
        push!(testAccList, acc)
        if isnan(objValue)
            break
        end
    end
    endTime = time()
    @printf("Finished training, time elapsed: %.4e\n", endTime - startTime)
    return server.W, objList, testAccList
end

# Implementation of the accelerated FedDCD algorithm
function accfedDCD(
    server::AccFedDCDServer,
    clients::Vector{AccFedDCDClient},
    numRounds::Int64;
    objMin::Float64 = 0.0
)
    # Connect clients with server
    server.clients = clients
    # Training process
    objList = zeros(Float64, 0)
    testAccList = zeros(Float64, 0)
    startTime = time()
    @printf("Start training!\n")
    for t = 1:numRounds
        # @printf("Round %d\n", t)
        # update v
        Threads.@threads for client in server.clients
            updatev!(client)
        end
        # first inner round
        select!(server)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            updateW!(client)
        end
        aggregate!(server)
        sendModel!(server, 1)
        # update u
        Threads.@threads for client in server.clients
            updateu!(client)
        end
        # second inner round
        select!(server)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            updateW!(client)
        end
        aggregate!(server)
        sendModel!(server, 2)
        # Print log
        objValue = getObjValue(server)
        acc = accuracy(server.Xtest, server.Ytest, server.W)
        @printf("Round : %4d, obj: %6.4e, acc: % 3.2f %%, time: %4.3f s\n", t, objValue, acc*100, time()-startTime)
        push!(objList, objValue)
        push!(testAccList, acc)
        if objValue < objMin
            break
        end
        if isnan(objValue)
            break
        end
    end
    endTime = time()
    @printf("Finished training, time elapsed: %.4e\n", endTime - startTime)
    return server.W, objList, testAccList
end

function accfedDCD(
    server::AccFedDCDServerNN,
    clients::Vector{AccFedDCDClientNN},
    numRounds::Int64;
    objMin::Float64 = 0.0
)
    # Connect clients with server
    server.clients = clients
    # Training process
    objList = zeros(Float64, 0)
    testAccList = zeros(Float64, 0)
    startTime = time()
    @printf("Start training!\n")
    for t = 1:numRounds
        # update v
        Threads.@threads for client in server.clients
            updatev!(client)
        end
        # first inner round
        select!(server)
        synchronize!(server)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            updateW!(client)
        end
        aggregate!(server)
        sendModel!(server, 1)
        # update u
        Threads.@threads for client in server.clients
            updateu!(client)
        end
        # second inner round
        select!(server)
        synchronize!(server)
        Threads.@threads for idx in server.selectedIndices
            client = server.clients[idx]
            updateW!(client)
        end
        aggregate!(server)
        sendModel!(server, 2)
        # Print log
        objValue = getObjValue(server)
        acc = accuracy(server.Xtest, server.Ytest, server.W)
        @printf("Round : %4d, obj: %6.4e, acc: % 3.2f %%, time: %4.3f s\n", t, objValue, acc*100, time()-startTime)
        push!(objList, objValue)
        push!(testAccList, acc)
        if objValue < objMin
            break
        end
        if isnan(objValue)
            break
        end
    end
    endTime = time()
    @printf("Finished training, time elapsed: %.4e\n", endTime - startTime)
    return server.W, objList, testAccList
end


text:
f(x) = (x - 3.0)^2.0
f2(x) = 2.0* (x -3.0)

x0 = 1

for i = 1:40
	@printf("%i\t%.4e\n",i,3.0 - x0)
	x0 = x0 - f(x0) / f2(x0)
end

ext: jl

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
using MyDG
using Test

@testset "MyDG.jl" begin
    # Write your tests here.
end


text:
__precompile__(true)
module DiDySim
    using LightGraphs, SparseArrays, IterTools, StatsBase, GZip, GraphIO, Combinatorics, Formatting, Printf
    include("types.jl")
    #include("NetWiring.jl")
    include("NetWiring2.jl")
    include("NetDynamics.jl")
    #include("NetAnalysis.jl")
    include("NetAnalysis2.jl")
    export BNode,Net,BNode2,Net2,reg_function
    export init_node!, init_net!, rand_init_net!, evolve_node,evolve_node!, evolve_net!, obs_rand_cond!
    export derrida_map!, derrida_map2!,derrida_map3!, full_derrida_map!, derrida_map_mod!, derrida_map_mod2!, derrida_map_mod3!
    export knock_out!,knock_down!,fix_state!, get_basins!,find_attractor!,find_attractor2!,derrida_map_att3!
    export generate_RBN,generate_RBN2,load_net2, load_net, to_boolean!, read_metadata
end


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb05od(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    LDA = NMAX
    NDIAG = 9
    LDWORK = NMAX*( 2*NMAX+NDIAG+1 )+NDIAG
    A = Array{Float64,2}(undef, LDA,NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    IWORK = Array{BlasInt,1}(undef, NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    DELTA = parse(Float64, replace(vs[2],'D'=>'E'))
    BALANC =  vs[3][1]
if ( N<=0 || N>NMAX )
else

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    MDIG, IDIG, INFO, IWARN = SLICOT.mb05od!(BALANC, N, NDIAG, DELTA, A, LDWORK)
    println(io, "MDIG = $MDIG")
    println(io, "IDIG = $IDIG")
    @test INFO == 0
    INFO == 0 || return
    println(io, "IWARN = $IWARN")

if ( INFO!=0 )
else
# interp output 1
    println(io,"A:")
    _nc = N
    _nr = N
    show(io,"text/plain",A[1:_nr,1:_nc])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


text:
using MatpowerCases, Base.Test

println("- extract_case check")
for c in casenames()
    mpc = loadcase(c; describe=false)
    @assert isa(extract_case(mpc), MatpowerCases.Case)
    @assert isa(extract_case(c), MatpowerCases.Case)
end

c = joinpath(Pkg.dir(),"MatpowerCases","data","case118.mat")
loadcase(c; describe=false)
@assert isa(extract_case(c), MatpowerCases.Case)


text:
using DEBmicroTrait
using CSV, DataFrames
using HypothesisTests, Statistics

df_isolates     = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/isolates2traits.csv", DataFrame, missingstring="N/A")
df_metabolites  = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/exudation_properties.csv", DataFrame, missingstring="N/A")

Min_gen_time    = df_isolates.Min_gen_time
gmax = log(2)./Min_gen_time
Gram_stain      = convert(Array{String,1}, df_isolates.gram_stain)

df_isolates.SAV = DEBmicroTrait.sav_from_growth_rate(gmax)
df_isolates.SAV = DEBmicroTrait.surface_area_volume_ratio(V_cell)
kcat = 180*60^2
df_isolates.sigma_p = @. gmax/(df_isolates.SAV*1e-12*kcat)

z_sugars        = reshape(convert(Array{Float64,1}, df_isolates.z_sugars./df_isolates.Genome_size*1e6),1,39)
z_organics      = reshape(convert(Array{Float64,1}, df_isolates.z_organic_acids./df_isolates.Genome_size*1e6),1,39)
z_aminos        = reshape(convert(Array{Float64,1}, df_isolates.z_amino_acids./df_isolates.Genome_size*1e6), 1,39)
z_fattys        = reshape(convert(Array{Float64,1}, df_isolates.z_fatty_acids./df_isolates.Genome_size*1e6),1,39)
z_nucleos       = reshape(convert(Array{Float64,1}, df_isolates.z_nucleotides./df_isolates.Genome_size*1e6),1,39)
z_auxins        = reshape(convert(Array{Float64,1}, df_isolates.z_auxins./df_isolates.Genome_size*1e6),1,39)
genome_distr    = vcat(z_sugars, z_organics, z_aminos, z_fattys, z_nucleos, z_auxins)

ρ_ps            = zeros(size(df_metabolites.Name,1), size(V_cell,1))

for j in 1:size(df_metabolites.Name,1)
    if df_metabolites.Ontology[j] == "Sugars"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            ρ_p = df_isolates.sigma_p[i]
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[1]
        end
    elseif df_metabolites.Ontology[j] == "Organic acids"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            ρ_p = df_isolates.sigma_p[i]
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[2]
        end
    elseif df_metabolites.Ontology[j] == "Amino acids"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            ρ_p = df_isolates.sigma_p[i]
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[3]
        end
    elseif df_metabolites.Ontology[j] == "Fatty acids"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            ρ_p = df_isolates.sigma_p[i]
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[4]
        end
    elseif df_metabolites.Ontology[j] == "Nucleotides"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            ρ_p = df_isolates.sigma_p[i]
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[5]
        end
    else
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            ρ_p = df_isolates.sigma_p[i]
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[6]
        end
    end
end

ρ_ps[ρ_ps.==0.0] .= 1e-12

N_C = zeros(size(df_metabolites.Name,1))
for i in 1:size(df_metabolites.Name,1)
    elementstring = convert(String, df_metabolites.Formula[i])
    N_C[i] = DEBmicroTrait.extract_composition(elementstring)[1]
end

N_SB              = DEBmicroTrait.transporter_density_to_monomer_uptake_sites(V_cell, ρ_ps, Min_gen_time, Gram_stain)
Vmax              = @. 180.0*60^2*N_SB.*N_C

D_S               = DEBmicroTrait.aqueous_diffusivity(df_metabolites.Molecular_weight)
K_D               = DEBmicroTrait.specific_reference_affinity(V_cell, ρ_ps, D_S)

a_s               = Vmax./K_D

df_assimilation = DataFrame()
df_assimilation.transporter_density = vec(ρ_ps)
df_assimilation.Vmax = vec(Vmax)
df_assimilation.KD = vec(K_D)
df_assimilation.affinity = vec(Vmax./K_D)

df_assimilation.ontology = repeat(df_metabolites.Ontology, size(V_cell,1))
response = Array{String}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     response[:,i] .= df_isolates.Rhizosphere_response[i]
 end

df_assimilation.response = vec(response)
genome_size = Array{Int}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     genome_size[:,i] .= df_isolates.Genome_size[i]
end
df_assimilation.genome_size = vec(genome_size)

cell_size = Array{Float64}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     cell_size[:,i] .= V_cell[i]
end
df_assimilation.cell_size = vec(cell_size)
min_gt = Array{Float64}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     min_gt[:,i] .= df_isolates.Min_gen_time[i]
end
df_assimilation.mingt = vec(min_gt)

SAV = Array{Float64}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     SAV[:,i] .= df_isolates.SAV[i]
end
df_assimilation.sav = vec(SAV)

df_assimilation.yield = vec(abs.(y_DEs))

CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/isolates_assimilation_sav.csv", df_assimilation)

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using Flux, Plots, Statistics, Distributions
using LaTeXStrings

x = range(-2, stop = 3, length = 50)
y = hcat(2.5 .* sin.(collect(x)), 2.5 .* cos.(collect(x)))
soft = [softmax(y[i, :])[1] for i=1:size(y,1)]

y_draws = Array{Float64}(50, 2, 1000)
for i in 1:1000
    #e = rand(Normal(0,1))
    y_draws[:, :, i] = hcat(2.5 .* sin.(collect(x)) .+ rand(Normal(0,1)) , 2.5 .* cos.(collect(x)) .+ rand(Normal(0,1)))
end
soft_draws = Array{Float64}(50, 1000)
for i in 1:50
    for j in 1:1000
        soft_draws[i, j] = softmax(y_draws[i, :, j])[1]
    end
end

soft_5 = [quantile(soft_draws[i, :], .025) for i=1:size(soft_draws, 1)]
soft_95 = [quantile(soft_draws[i, :], .975) for i=1:size(soft_draws, 1)]
soft_mean = [mean(soft_draws[i, :]) for i=1:size(soft_draws, 1)]

y_5 = hcat([quantile(y_draws[i, 1, :], .025) for i=1:size(y_draws, 1)],
        [quantile(y_draws[i, 2, :], .025) for i=1:size(y_draws, 1)])
y_95 = hcat([quantile(y_draws[i, 1, :], .975) for i=1:size(y_draws, 1)],
        [quantile(y_draws[i, 2, :], .975) for i=1:size(y_draws, 1)])
y_mean = hcat([mean(y_draws[i, 1, :]) for i=1:size(y_draws, 1)],
        [mean(y_draws[i, 2, :]) for i=1:size(y_draws, 1)])

plot(x, y[:, 1], label = L"\sigma(f(x))", legend=:bottomright,
    ribbon = (y[:, 1]-y_5[:, 1], y_95[:, 1]-y[:, 1]),
    color = :black,
    fillcolor = :grey)

plot(x, soft, label = L"\sigma(f(x))", legend=:bottomright,
    color = :red,
    fillcolor = :red)
savefig("uncertainty_example_post_softmax_slide1.png")

plot(x, soft, label = L"\sigma(f(x))", legend=:bottomright,
    ribbon = (soft-soft_5, soft_95-soft),
    color = :red,
    fillcolor = :red)
savefig("uncertainty_example_post_softmax_slide2.png")

plot(x, y[:, 1], label = L"f(x)", legend=:bottomright,
    color = :red,
    fillcolor = :red)
plot!(x, y[:, 2], label = L"\tilde{f}(x)", legend=:bottomright,
    color = :blue,
    fillcolor = :blue)
savefig("uncertainty_example_pre_softmax_slide1")

plot(x, y[:, 1], label = L"\sigma(f(x))", legend=:bottomright,
    color = :red,
    fillcolor = :red,
    ribbon = (y[:, 1]-y_5[:, 1], y_95[:, 1]-y[:, 1]))
plot!(x, y[:, 2], label = L"\tilde{f}(x)", legend=:bottomright,
        color = :blue,
        fillcolor = :blue,
        ribbon = (y[:, 2]-y_5[:, 2], y_95[:, 2]-y[:, 2]))
savefig("uncertainty_example_pre_softmax_slide2")

max_stars_repo_name: ellisbrown/BNN-Uncertainty

max_issues_repo_name: ellisbrown/BNN-Uncertainty

max_forks_repo_name: ellisbrown/BNN-Uncertainty


text:
using LinearAlgebra
using Polynomials
using Optim
using Distributions
using LsqFit
using Plots
using CSV
using DifferentialEquations
using DelimitedFiles
using Dates
using DataFrames

datas=CSV.File("../data/freq/Data_cleansed/TransnetBW/2019.csv",header=0)
fs=(datas[:,2].-50).*2*pi
fsbeg=Array{Float64,1}(undef,90)
smoothfs=Array{Float64,1}(undef,900)
smoothfs2=Array{Float64,1}(undef,900)
i_smoothfs=Array{Float64,1}(undef,300)
i_smoothfs2=Array{Float64,1}(undef,300)
i_smoothfs3=Array{Float64,1}(undef,300)
posOfdbcrossing=-1
var_ww_0 = 0.0
var_aa_0 = 0.0
rho_aw_0 = 0.0
D_guess = 0.0
minpos2=1
minpos=1
solRev=Array{Float64,2}(undef,900,3)
varOme=Array{Float64,1}(undef,900)
parguessRev=[0.02,0.0005,D_guess,var_ww_0,var_aa_0,rho_aw_0]
LB = [0.001,0.0,0.0,(var_ww_0/500), 0.0, -20.0]
UB = [0.04,1.0e-3,1.0,(var_ww_0*500), +120.0, +20.0]
parguessRes=[0.02,0.0005,D_guess,var_ww_0,var_aa_0,rho_aw_0]
q = 0.0
omega_0 = 0.0
par_guess1=[q,(q/900),omega_0,0.0,0.0,0.0]
LB1 = [q*1.2,0.0,omega_0-0.075,-320.0,0.0,0.0]
UB1 = [q*0.83,1.0,omega_0+0.075,320.0,0.0,0.0]
meanOme=Array{Float64,1}(undef,900)
sol=Array{Float64,2}(undef,900,2)
par_guess1Res=[q,(q/900),omega_0,0.0,0.0,0.0]

MLE = Array{Float64,1}(undef,10) #[gamma1,gamma2,q,r,D,mu_w_0,mu_a_0,var_ww_0,var_aa_0,rho_aw_0]
LB2 = [0.0, 0.0, -1.0,0.0,-1.0,-2.0,0.0, 50.0, 0.0, -20.0]
UB2 = [0.45,1.0e-3,0.0,1.0,1.0, 2.0,0.0, 50.0, +120.0, +20.0]
MLENoC = Array{Float64,1}(undef,8) #[gamma1,gamma2,q,r,D,mu_w_0,mu_a_0,var_ww_0,var_aa_0,rho_aw_0]
LB2NoC = [-1.0,0.0,-1.0,-2.0,0.0, 50.0, 0.0, -20.0]
UB2NoC = [0.0,1.0,1.0, 2.0,0.0, 50.0, +120.0, +20.0]
MS=Array{Float64,2}(undef,900,2)
sol4=Array{Float64,2}(undef,900,5)
sol5=Array{Float64,2}(undef,900,2)
MLEres = Array{Float64,1}(undef,10)
info=Array{Float64,2}(undef,8784,16)
t0=1000
turnOn=20
rocof=0

function dyn_moments1Rev!(dy,y,par,t)  #par will look like: [gamma1,gamma2,D]
    dy[1] = par[3]^2-(2*par[1]*y[1])-(2*par[2]*y[3])
    dy[2] = 2*y[3]
    dy[3] = -par[1]*y[3]-par[2]*y[2]+y[1]
end
function dyn_moments1_2Rev!(dy,y,par,t)  #par will look like: [gamma1,gamma2,D]
    dy[1] = par[3]^2-(2*par[1]*((t-posOfdbcrossing+1)/turnOn)*y[1])
    dy[2] = 2*y[3]
    dy[3] = -par[1]*((t-posOfdbcrossing+1)/turnOn)*y[3]+y[1]
end
function dyn_moments0Rev!(dy,y,par,t)  #par will look like: [gamma1,gamma2,D]
    dy[1] = par[3]^2
    dy[2] = 2*y[3]
    dy[3] = y[1]
end
function halfDetTimeseries_fpeRev(par) #par: [gamma1,gamma2,D,var_ww_0,var_aa_0,cov_wa_0]
    if posOfdbcrossing==0       #always in db-->control off
        so=solve(ODEProblem(dyn_moments0Rev!,[par[4],par[5],par[6]],(0.0,899.0),[par[1],par[2],par[3]]))
        for j in 1:900
            global solRev[j,:]=so(j-1.0)
        end
    else                             # posOfdbcrossing>0 (first inside db, then outside)
        so=solve(ODEProblem(dyn_moments0Rev!,[par[4],par[5],par[6]],(0.0,posOfdbcrossing-1.0),[par[1],par[2],par[3]]))
        for j in 1:posOfdbcrossing
            global solRev[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments1_2Rev!,[solRev[posOfdbcrossing,1],solRev[posOfdbcrossing,2],solRev[posOfdbcrossing,3]],(posOfdbcrossing-1.0,posOfdbcrossing+(turnOn-1.0)),[par[1],par[2],par[3]]))
        for j in posOfdbcrossing:(posOfdbcrossing+turnOn)
            global solRev[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments1Rev!,[solRev[(posOfdbcrossing+turnOn),1],solRev[(posOfdbcrossing+turnOn),2],solRev[(posOfdbcrossing+turnOn),3]],(posOfdbcrossing+(turnOn-1.0),899.0),[par[1],par[2],par[3]]))
        for j in (posOfdbcrossing+turnOn):900
            global solRev[j,:]=so(j-1.0)
        end
    end
    return solRev[:,1]  #out: array containing var(w)
end
function dyn_ode1!(dy,y,par,t) #par=[gamma1,gamma2,q,r]  y=[w,a]
    dy[1] = par[3]+(par[4]*t)-par[1]*y[1]-par[2]*y[2]
    dy[2] = y[1]
end
function dyn_ode1_2!(dy,y,par,t) #par=[gamma1,gamma2,q,r]  y=[w,a]
    dy[1] = par[3]+(par[4]*t)-par[1]*((t-posOfdbcrossing+1)/turnOn)*y[1]
    dy[2] = y[1]
end
function dyn_ode0!(dy,y,par,t) #par=[gamma1,gamma2,q,r]  y=[w,a]
    dy[1] = par[3]+(par[4]*t)
    dy[2] = y[1]
end
function timeseries_ode(par)
    if posOfdbcrossing==0       #always in db-->control off
        so=solve(ODEProblem(dyn_ode0!,[par[3],par[4]],(0.0,899.0),[par[5],par[6],par[1],par[2]]))
        for j in 1:900
            global sol[j,:]=so(j-1.0)
        end
    else                             # posOfdbcrossing>0 (first inside db, then outside)
        so=solve(ODEProblem(dyn_ode0!,[par[3],par[4]],(0.0,posOfdbcrossing-1),[par[5],par[6],par[1],par[2]]))
        for j in 1:posOfdbcrossing
            global sol[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_ode1_2!,[sol[posOfdbcrossing,1],sol[posOfdbcrossing,2]],(posOfdbcrossing-1.0,posOfdbcrossing+(turnOn-1.0)),[par[5],par[6],par[1],par[2]]))
        for j in posOfdbcrossing:(posOfdbcrossing+turnOn)
            global sol[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_ode1!,[sol[(posOfdbcrossing+turnOn),1],sol[(posOfdbcrossing+turnOn),2]],(posOfdbcrossing+(turnOn-1.0),899.0),[par[5],par[6],par[1],par[2]]))
        for j in (posOfdbcrossing+turnOn):900
            global sol[j,:]=so(j-1.0)
        end
    end
    return sol[:,1]
end
function dyn_moments12!(dy,y,par,t)  #par will look like: [gamma1,gamma2,q,r,D]
    dy[1] = par[3]+(par[4]*t)-par[1]*y[1]-par[2]*y[2]
    dy[2] = y[1]
    dy[3] = par[5]^2-(2*par[1]*y[3])-(2*par[2]*y[5])
    dy[4] = 2*y[5]
    dy[5] = -par[1]*y[5]-par[2]*y[4]+y[3]
end
function dyn_moments1_22!(dy,y,par,t)  #par will look like: [gamma1,gamma2,q,r,D]
    dy[1] = par[3]+(par[4]*t)-par[1]*((t-posOfdbcrossing+1)/turnOn)*y[1]
    dy[2] = y[1]
    dy[3] = par[5]^2-(2*par[1]*((t-posOfdbcrossing+1)/turnOn)*y[3])
    dy[4] = 2*y[5]
    dy[5] = -par[1]*((t-posOfdbcrossing+1)/turnOn)*y[5]+y[3]
end
function dyn_moments02!(dy,y,par,t)  #par will look like: [gamma1,gamma2,q,r,D]
    dy[1] = par[3]+(par[4]*t)
    dy[2] = y[1]
    dy[3] = par[5]^2
    dy[4] = 2*y[5]
    dy[5] = y[3]
end
function timeseries_fpe(par)
    if posOfdbcrossing==0       #always in db-->control off
        so=solve(ODEProblem(dyn_moments02!,[par[6],par[7],par[8],par[9],par[10]],(0.0,899.0),par[1:5]))
        for j in 1:900
            global sol4[j,:]=so(j-1.0)
        end
    else                             # posOfdbcrossing>0 (first inside db, then outside)
        so=solve(ODEProblem(dyn_moments02!,[par[6],par[7],par[8],par[9],par[10]],(0.0,posOfdbcrossing-1),par[1:5]))
        for j in 1:posOfdbcrossing
            global sol4[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments1_22!,[sol4[posOfdbcrossing,1],sol4[posOfdbcrossing,2],sol4[posOfdbcrossing,3],sol4[posOfdbcrossing,4],sol4[posOfdbcrossing,5]],(posOfdbcrossing-1.0,posOfdbcrossing+(turnOn-1.0)),par[1:5]))
        for j in posOfdbcrossing:(posOfdbcrossing+turnOn)
            global sol4[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments12!,[sol4[(posOfdbcrossing+turnOn),1],sol4[(posOfdbcrossing+turnOn),2],sol4[(posOfdbcrossing+turnOn),3],sol4[(posOfdbcrossing+turnOn),4],sol4[(posOfdbcrossing+turnOn),5]],(posOfdbcrossing+(turnOn-1.0),899.0),par[1:5]))
        for j in (posOfdbcrossing+turnOn):900
            global sol4[j,:]=so(j-1.0)
        end
    end
    global sol5[:,1]=sol4[:,1]
    global sol5[:,2]=sol4[:,3]
    return sol5
end
function dyn_moments02NoC!(dy,y,par,t)  #par will look like: [q,r,D] WHY THIS FUNCTION? if we do not have control
    dy[1] = par[1]+(par[2]*t)
    dy[2] = y[1]
    dy[3] = par[3]^2
    dy[4] = 2*y[5]
    dy[5] = y[3]
end
function timeseries_fpeNoC(par)
    #always in db-->control off
    so=solve(ODEProblem(dyn_moments02NoC!,[par[4],par[5],par[6],par[7],par[8]],(0.0,899.0),par[1:3]))
    for j in 1:900
        global sol4[j,:]=so(j-1.0)
    end
    global sol5[:,1]=sol4[:,1]
    global sol5[:,2]=sol4[:,3]
    return sol5
end

# how to find the right starting point?

for h in 1:500
    if Second(DateTime(datas[(3600*(h-1)+1),1][1:10]*"T"*datas[(3600*(h-1)+1),1][12:19]))!=Second(0)
        ch=0
        for i in -20:20
            if Second(DateTime(datas[(3600*(h-1)+1+i),1][1:10]*"T"*datas[(3600*(h-1)+1+i),1][12:19]))==Second(0)
                ch=i
            end
        end
        global fs[(3600*(h-1)+1):(length(fs)-maximum([ch,0]))]=fs[(3600*(h-1)+1+ch):length(fs)+minimum([ch,0])]
        global datas[(3600*(h-1)+1):(length(fs)-maximum([ch,0])),:]=datas[(3600*(h-1)+1+ch):length(fs)+minimum([ch,0]),:]
    end
    badData=0
    for i in (3600*(h-1)+1):(3600*(h-1)+900)
        if (typeof(fs[i])==Float64 && isnan(fs[i])==false)==false # if data is missing or NaN set badData to 1
            badData=1
            break
        end
    end
    if badData!=1
        function halfDetLoglikeRev(par::AbstractArray) #par: [gamma1,gamma2,D,var_ww_0,var_aa_0,cov_wa_0]
            global varOme=halfDetTimeseries_fpeRev(par)
            global lnp = -sum((-1/2)*log.(2*pi*abs.(varOme))-(((fs[(3600*(h-1)+1):(3600*(h-1)+900)]-smoothfs2).^2)./(2*abs.(varOme))))
            if minimum(varOme)<1.0e-10
                global lnp=10000000000000000000000000
            end
            return lnp
        end
        function halfDetLoglikeODE(par::AbstractArray) #[q,r,omega_0,theta_0,gamma1,gamma2]
            global meanOme=timeseries_ode(par)
            global varOme=halfDetTimeseries_fpeRev(parguessRes)
            global lnp = -sum((-1/2)*log.(2*pi*abs.(varOme))-(((smoothfs2-meanOme).^2)./(2*abs.(varOme))))
            if minimum(varOme)<1.0e-10
                global lnp=10000000000000000000000000
            end
            return lnp
        end
        function loglike(par::AbstractArray)
            global MS=timeseries_fpe(par)
            global lnp = -sum((-1/2)*log.(2*pi*abs.(MS[:,2]))-(((fs[(3600*(h-1)+1):(3600*(h-1)+900)]-MS[:,1]).^2)./(2*abs.(MS[:,2]))))
            if minimum(MS[:,2])<1.0e-10
                global lnp=10000000000000000000000000
            end
            return lnp
        end
        function loglikeNoC(par::AbstractArray)
            global MS=timeseries_fpeNoC(par)
            global lnp = -sum((-1/2)*log.(2*pi*abs.(MS[:,2]))-(((fs[(3600*(h-1)+1):(3600*(h-1)+900)]-MS[:,1]).^2)./(2*abs.(MS[:,2]))))
            if minimum(MS[:,2])<1.0e-10
                global lnp=10000000000000000000000000
            end
            return lnp
        end
        try
            for i in 1:70
                global smoothfs[i]=sum(fs[(((h-1)*3600)+1):(((h-1)*3600)+(2*i))])/(2*i)
            end
            for i in 71:830
                global smoothfs[i]=sum(fs[(((h-1)*3600)+(i-70)):(((h-1)*3600)+(i+70))])/140
            end
            for i in 831:900
                global smoothfs[i]=sum(fs[(((h-1)*3600)+900-(2*(900-i))):(((h-1)*3600)+900)])/(2*(900-i)+1)
            end
            for i in 880:900
                global smoothfs[i]=smoothfs[879]
            end
            for i in 1:40
                global smoothfs2[i]=sum(smoothfs[1:(2*i)])/(2*i)
            end
            for i in 41:860
                global smoothfs2[i]=sum(smoothfs[(i-40):(i+40)])/80
            end
            for i in 861:900
                global smoothfs2[i]=sum(smoothfs[(900-(2*(900-i))):900])/(2*(900-i))
            end
            for i in 880:900
                global smoothfs2[i]=smoothfs2[879]
            end
            global fsbeg[1:30].=fs[((h-1)*3600)+1]
            global fsbeg[31:90]=fs[((h-1)*3600+1):((h-1)*3600+60)]
            for i in 1:30
                global i_smoothfs2[i]=sum(fsbeg[i:(i+60)])/61
                global i_smoothfs[i]=sum(fsbeg[26+i:34+i])/9
            end
            for i in 31:300
                global i_smoothfs[i]=sum(fs[(((h-1)*3600)+(i-4)):(((h-1)*3600)+(i+4))])/9
                global i_smoothfs2[i]=sum(fs[(((h-1)*3600)+(i-30)):(((h-1)*3600)+(i+30))])/61
            end
            for i in 11:290
                i_smoothfs3[i]=sum(i_smoothfs2[(i-10):(i+10)])/21
            end
            global i_smoothfs3[1:10]=i_smoothfs2[1:10]
            global i_smoothfs3[291:300].=i_smoothfs3[290]
            if polyfit(5:15,i_smoothfs3[5:15],1)[1]>0
                for i in 11:299
                    if i_smoothfs3[i+1]<i_smoothfs3[i]
                        global minpos2=i
                        break
                    end
                end
                for i in maximum([(minpos2-25),1]):minimum([(minpos2+25),300])
                    if i_smoothfs[i+1]<i_smoothfs[i]
                        global minpos=i
                        break
                    end
                end
            else
                for i in 11:299
                    if i_smoothfs3[i+1]>i_smoothfs3[i]
                        global minpos2=i
                        break
                    end
                end
                for i in maximum([(minpos2-25),1]):minimum([(minpos2+25),300])
                    if i_smoothfs[i+1]>i_smoothfs[i]
                        global minpos=i
                        break
                    end
                end
            end
            if minpos>14
                global rocof=polyfit(1:minpos,i_smoothfs[1:minpos],1)[1]
                global turnOn=minimum([30,convert(Int64,abs(round((20/0.0045)*rocof)))])
            end
            #plot(-100:899,fs[((h-1)*3600+1-100):((h-1)*3600+900)])
            #plot!(0:299,i_smoothfs3)
            #plot!(0:299,i_smoothfs2)
            #plot!(0:299,i_smoothfs)
            #plot!(0:899,smoothfs2)
            #plot!(0:minpos-1,rocof(1:minpos))
            global posOfdbcrossing=maximum([1,convert(Int64,round(minpos-(turnOn/2)))])
            if turnOn<5
                global posOfdbcrossing = 0
            elseif minpos<15
                global posOfdbcrossing = 0
            elseif maximum(abs.(i_smoothfs[(minpos-10):(minpos+10)]))<0.13
                global posOfdbcrossing = 0
            end
            global var_ww_0 = (std(fs[((h-1)*3600+1):((h-1)*3600+900)]-smoothfs2))^2
            global var_aa_0 = var_ww_0/10
            global rho_aw_0 = var_ww_0
            global D_guess = sqrt(var_ww_0*2*0.02)
            if posOfdbcrossing==0
                global omega_0=i_smoothfs3[1]
                global MLENoC=[0.0,0.0,D_guess,omega_0,0.0,var_ww_0,var_aa_0,rho_aw_0]
                global LB2NoC = [-0.5,-0.05,0.0,omega_0-0.05,-10.0,(var_ww_0/500), 0.0, -20.0]
                global UB2NoC = [0.5, 0.05, 1.0,omega_0+0.05, 10.0, (var_ww_0*500), +120.0, +20.0]
                global MLEresNoC = optimize(loglikeNoC, LB2NoC, UB2NoC, MLENoC, Fminbox(LBFGS()),).minimizer
                global info[h,:]=[0.0,0.0,0.0,0.0,MLEresNoC[1],MLEresNoC[2],MLEresNoC[3],MLEresNoC[4],MLEresNoC[5],MLEresNoC[6],MLEresNoC[7],MLEresNoC[8],0,maximum(abs.(fs[(3600*(h-1)+1):(3600*(h-1)+900)])),maximum(abs.(smoothfs[1:450])),-loglikeNoC(MLEresNoC)]
                #[posOfdbcrossing,MLE,RoCof,nadir,MeanDevInFirstHalf,Loglike]
            else
                global parguessRev=[0.003,0.00005,D_guess,var_ww_0,var_aa_0,rho_aw_0]
                global LB = [0.001,0.0,0.0,(var_ww_0/500), 0.0, -20.0]
                global UB = [0.04,1.0e-3,1.0,(var_ww_0*500), +120.0, +20.0]
                global parguessRes = optimize(halfDetLoglikeRev, LB, UB, parguessRev, Fminbox(BFGS())).minimizer #Optim.Options(show_trace=true,iterations=100)
                #par: [gamma1,gamma2,D,var_ww_0,var_aa_0,cov_wa_0]
                #plot(0:899,fs[((h-1)*3600+1):((h-1)*3600+900)])
                #plot!(0:899,smoothfs2)
                #plot!(0:899,(smoothfs2+sqrt.(halfDetTimeseries_fpeRev(parguessRes))))
                #plot!(0:899,(smoothfs2-sqrt.(halfDetTimeseries_fpeRev(parguessRes))))
                par_temp = polyfit(1:posOfdbcrossing,i_smoothfs[1:posOfdbcrossing],1)
                #plot!(1:posOfdbcrossing,par_temp(1:posOfdbcrossing))
                global q = par_temp[1]
                global omega_0 = par_temp[0]
                global par_guess1=[q,(-q/600),omega_0,0.0,parguessRes[1],parguessRes[2]] #[q,r,omega_0,theta_0,gamma1,gamma2]
                if q<0.0
                    global LB1 = [q*1.01,-q*0.83/3000,omega_0-0.025,-10.0,0.0001,0.0]
                    global UB1 = [q*0.99,-q*1.2/300,omega_0+0.025,10.0,0.04,1.0e-3]
                else
                    global LB1 = [q*0.99,-q*1.2/300,omega_0-0.025,-10.0,0.0001,0.0]
                    global UB1 = [q*1.01,-q*0.83/3000,omega_0+0.025,10.0,0.04,1.0e-3]
                end
                global par_guess1Res = optimize(halfDetLoglikeODE, LB1, UB1, par_guess1, Fminbox(LBFGS())).minimizer #Optim.Options(show_trace=true,iterations=100)
                #plot(0:899,fs[((h-1)*3600+1):((h-1)*3600+900)])
                #plot!(0:899,smoothfs2)
                #plot!(0:899,timeseries_ode(par_guess1Res))
                #plot!(0:899,(timeseries_ode(par_guess1Res)+sqrt.(halfDetTimeseries_fpeRev(parguessRes))))
                #plot!(0:899,(timeseries_ode(par_guess1Res)-sqrt.(halfDetTimeseries_fpeRev(parguessRes))))
                global MLE=[par_guess1Res[5],par_guess1Res[6],par_guess1Res[1],par_guess1Res[2],parguessRes[3],par_guess1Res[3],par_guess1Res[4],parguessRes[4],parguessRes[5],parguessRes[6]]
                if q<0.0
                    global LB2 = [0.0001, 0.0, q*1.05,-q*0.83/3000,0.0,omega_0-0.025,-10.0, (var_ww_0/500), 0.0, -20.0]
                    global UB2 = [0.04,1.0e-3,q*0.95,-q*1.2/300,1.0, omega_0+0.025,10.0, (var_ww_0*500), +120.0, 20.0]
                else
                    global LB2 = [0.0001, 0.0, q*0.95,-q*1.2/300,0.0,omega_0-0.025,-10.0, (var_ww_0/500), 0.0, -20.0]
                    global UB2 = [0.04,1.0e-3,q*1.05,-q*0.83/3000,1.0, omega_0+0.025,10.0,(var_ww_0*500), +120.0, +20.0]
                end
                global MLEres = optimize(loglike, LB2, UB2, MLE, Fminbox(LBFGS()),).minimizer
                #[posOfdbcrossing,MLE,RoCof,nadir,MeanDevInFirstHalf,Loglike]
                global info[h,:]=[turnOn,posOfdbcrossing,MLEres[1],MLEres[2],MLEres[3],MLEres[4],MLEres[5],MLEres[6],MLEres[7],MLEres[8],MLEres[9],MLEres[10],q,maximum(abs.(fs[(3600*(h-1)+1):(3600*(h-1)+900)])),maximum(abs.(smoothfs[1:450])),-loglike(MLEres)]
            end
        catch y
            global info[h,:]=[NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,maximum(abs.(fs[(3600*(h-1)+1):(3600*(h-1)+900)])),NaN,NaN]
        end
    else
        global info[h,:]=[NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN]
    end
    println(h)
end
info[1,:]
h=1
plot(0:899,fs[((h-1)*3600+1):((h-1)*3600+900)])
plot!(0:899,timeseries_fpeNoC(info[1,5:12])[:,1])
#######################################
#######################################
#######################################

# if i find that the secnd does not match:
# move all the data such that it fits from nov on:
# instead of 00:00:00 we find 23:59:59 --> "overwrite one entry" by: global fs[(3600*(h-1)+1):end-1]=fs[(3600*(h-1)+2):end]
# instead of 00:00:00 we find 00:00:01 --> "delete one entry" by: global fs[(3600*(h-1)+1):end]=fs[(3600*(h-1)):end-1]   and same for the DateTime

h=1
for i in 1:70
    global smoothfs[i]=sum(fs[(((h-1)*3600)+1):(((h-1)*3600)+(2*i))])/(2*i)
end
for i in 71:830
    global smoothfs[i]=sum(fs[(((h-1)*3600)+(i-70)):(((h-1)*3600)+(i+70))])/140
end
for i in 831:900
    global smoothfs[i]=sum(fs[(((h-1)*3600)+900-(2*(900-i))):(((h-1)*3600)+900)])/(2*(900-i)+1)
end
for i in 880:900
    global smoothfs[i]=smoothfs[879]
end
for i in 1:40
    global smoothfs2[i]=sum(smoothfs[1:(2*i)])/(2*i)
end
for i in 41:860
    global smoothfs2[i]=sum(smoothfs[(i-40):(i+40)])/80
end
for i in 861:900
    global smoothfs2[i]=sum(smoothfs[(900-(2*(900-i))):900])/(2*(900-i))
end
for i in 880:900
    global smoothfs2[i]=smoothfs2[879]
end
global fsbeg[1:30].=fs[((h-1)*3600)+1]
global fsbeg[31:90]=fs[((h-1)*3600+1):((h-1)*3600+60)]
for i in 1:30
    global i_smoothfs2[i]=sum(fsbeg[i:(i+60)])/61
    global i_smoothfs[i]=sum(fsbeg[26+i:34+i])/9
end
for i in 31:300
    global i_smoothfs[i]=sum(fs[(((h-1)*3600)+(i-4)):(((h-1)*3600)+(i+4))])/9
    global i_smoothfs2[i]=sum(fs[(((h-1)*3600)+(i-30)):(((h-1)*3600)+(i+30))])/61
end
for i in 11:290
    i_smoothfs3[i]=sum(i_smoothfs2[(i-10):(i+10)])/21
end
global i_smoothfs3[1:10]=i_smoothfs2[1:10]
global i_smoothfs3[291:300].=i_smoothfs3[290]
if polyfit(5:15,i_smoothfs3[5:15],1)[1]>0
    for i in 11:299
        if i_smoothfs3[i+1]<i_smoothfs3[i]
            global minpos2=i
            break
        end
    end
    for i in maximum([(minpos2-25),1]):minimum([(minpos2+25),300])
        if i_smoothfs[i+1]<i_smoothfs[i]
            global minpos=i
            break
        end
    end
else
    for i in 11:299
        if i_smoothfs3[i+1]>i_smoothfs3[i]
            global minpos2=i
            break
        end
    end
    for i in maximum([(minpos2-25),1]):minimum([(minpos2+25),300])
        if i_smoothfs[i+1]>i_smoothfs[i]
            global minpos=i
            break
        end
    end
end
if minpos>14
    global rocof=polyfit(1:minpos,i_smoothfs[1:minpos],1)[1]
    global turnOn=minimum([30,convert(Int64,abs(round((20/0.0045)*rocof)))])
end
plot(0:899,fs[((h-1)*3600+1):((h-1)*3600+900)])
plot!(0:299,i_smoothfs3)
plot!(0:299,i_smoothfs2)
plot!(0:299,i_smoothfs)
#plot!(0:899,smoothfs2)
#plot!(0:minpos-1,rocof(1:minpos))
global posOfdbcrossing=maximum([1,convert(Int64,round(minpos-(turnOn/2)))])
if turnOn<5
    global posOfdbcrossing = 0
elseif minpos<15
    global posOfdbcrossing = 0
elseif maximum(abs.(i_smoothfs[(minpos-10):(minpos+10)]))<0.13
    global posOfdbcrossing = 0
end


text:
"""
    WeightedSample <: AbstractActionSelector
    WeightedSample()
"""
struct WeightedSample <: AbstractActionSelector end

"""
    (p::WeightedSample)(values::AbstractArray)

!!! note
    Action `values` are normalized to have a sum of 1.0
    and then used as the probability to sample a random action.
"""
(p::WeightedSample)(values::AbstractArray) = sample(weights(values))


text:
"""
    LocalGeometry

The necessary local metric information defined at each node.
"""
struct LocalGeometry{I, C <: AbstractPoint, FT, S}
    "Coordinates of the current point"
    coordinates::C
    "Jacobian determinant of the transformation `ξ` to `x`"
    J::FT
    "Metric terms: `J` multiplied by the quadrature weights"
    WJ::FT
    "Partial derivatives of the map from `ξ` to `x`: `∂x∂ξ[i,j]` is ∂xⁱ/∂ξʲ"
    ∂x∂ξ::Axis2Tensor{FT, Tuple{LocalAxis{I}, CovariantAxis{I}}, S}
    "Partial derivatives of the map from `x` to `ξ`: `∂ξ∂x[i,j]` is ∂ξⁱ/∂xʲ"
    ∂ξ∂x::Axis2Tensor{FT, Tuple{ContravariantAxis{I}, LocalAxis{I}}, S}
end

LocalGeometry(coordinates, J, WJ, ∂x∂ξ) =
    LocalGeometry(coordinates, J, WJ, ∂x∂ξ, inv(∂x∂ξ))

"""
    SurfaceGeometry

The necessary local metric information defined at each node on each surface.
"""
struct SurfaceGeometry{FT, N}
    "surface Jacobian determinant, multiplied by the surface quadrature weight"
    sWJ::FT
    "surface outward pointing normal vector"
    normal::N
end

undertype(::Type{LocalGeometry{I, C, FT, S}}) where {I, C, FT, S} = FT
undertype(::Type{SurfaceGeometry{FT, N}}) where {FT, N} = FT


text:
module TestBetaDivOperations
  using Test
  using EcologicalNetworks
  using LinearAlgebra

  A = UnipartiteNetwork(rand(Bool, (10, 10)))
  B = UnipartiteNetwork(rand(Bool, (10, 10)))
  U = union(A,B)
  NI = intersect(A,B)
  @test typeof(U) == typeof(B)
  @test links(U) >= links(A)
  @test links(U) >= links(B)
  @test links(NI) <= links(B)
  @test links(NI) <= links(A)
  @test richness(NI) == richness(B)

  A = BipartiteNetwork(rand(Bool, (10, 10)))
  B = BipartiteNetwork(rand(Bool, (20, 20)))
  U = union(A,B)
  NI = intersect(A,B)
  @test typeof(U) == typeof(B)
  @test links(U) >= links(A)
  @test links(U) >= links(B)
  @test richness(U) == richness(B)
  @test links(NI) <= links(B)
  @test links(NI) <= links(A)
  @test richness(NI) == richness(A)

  # No common species
  A = UnipartiteNetwork(Matrix(I, (3,3)), [:a, :b, :c])
  B = UnipartiteNetwork(Matrix(I, (3,3)), [:d, :e, :f])

  @test species(setdiff(B, A)) == species(B)
  @test species(setdiff(A, B)) == species(A)
  @test richness(intersect(A, B)) == 0
  @test links(intersect(A, B)) == 0

  A = BipartiteNetwork(Matrix(I, (3,3)), [:a, :b, :c], [:r, :t, :u])
  B = BipartiteNetwork(Matrix(I, (3,3)), [:d, :e, :f], [:j, :k, :l])

  @test species(setdiff(B, A)) == species(B)
  @test species(setdiff(A, B)) == species(A)
  @test richness(intersect(A, B)) == 0
  @test links(intersect(A, B)) == 0

  # Measures
  for m in [KGL01, KGL02, KGL03, KGL04, KGL05, KGL06, KGL07, KGL08, KGL09, KGL10,
      KGL11, KGL12, KGL13, KGL14, KGL15, KGL16, KGL17, KGL18, KGL19, KGL20, KGL21,
      KGL22, KGL23, KGL24]
      m(βos(A, B))
  end

end


text:
using Test
@time @testset "All Tests" begin
    include("test_all.jl")
end


text:
@testset "Moments" begin
    f = zeros(5,5,9)
    f1D = zeros(30,3)
    height = zeros(5,5)
    hei = zeros(30)
    velx = zeros(5,5)
    vely = zeros(5,5)
    vel = zeros(30)
    state = Swalbe.Sys(Swalbe.SysConst(Lx=5, Ly=5), "CPU")
    state1D = Swalbe.Sys(Swalbe.SysConst_1D(L=30))

    @testset "No velocity" begin
        f[:,:,1] .= 1.0
        state.fout[:,:,1] .= 1.0
        Swalbe.moments!(height, velx, vely, f)
        Swalbe.moments!(state)
        @test all(height .== 1.0)
        @test state.height[1,1] ≈ 1.0
        f1D[:,1] .= 1.0
        state1D.fout[:,1] .= 1.0
        Swalbe.moments!(hei, vel, f1D)
        Swalbe.moments!(state1D)
        @test all(hei .== 1.0)
        @test all(state1D.height .== 1.0)
    end
    @testset "Velocity" begin
        # 2D
        f[:,:,1] .= 1.0
        f[:,:,2] .= 0.1
        state.fout[:,:,2] .= 0.1
        Swalbe.moments!(height, velx, vely, f)
        Swalbe.moments!(state)
        @test all(height .== 1.1)
        @test all(velx .== 0.1/1.1)
        @test all(state.height .== 1.1)
        @test all(state.velx .== 0.1/1.1)
        # 1D
        f1D[:,1] .= 1.0
        f1D[:,2] .= 0.1
        state1D.fout[:,2] .= 0.1
        Swalbe.moments!(hei, vel, f1D)
        @test all(hei .== 1.1)
        @test all(vel .== 0.1/1.1)
        Swalbe.moments!(state1D)
        @test all(state1D.height .== 1.1)
        @test all(state1D.vel .== 0.1/1.1)
        # 2D
        f[:,:,3] .= 0.2
        state.fout[:,:,3] .= 0.2
        Swalbe.moments!(height, velx, vely, f)
        Swalbe.moments!(state)
        @test all(height .== 1.3)
        @test all(velx .== 0.1/1.3)
        @test all(vely .== 0.2/1.3)
        @test all(state.height .== 1.3)
        @test all(state.velx .== 0.1/1.3)
        @test all(state.vely .== 0.2/1.3)
        # 1D
        f1D[:,3] .= 0.2
        state1D.fout[:,3] .= 0.2
        Swalbe.moments!(hei, vel, f1D)
        Swalbe.moments!(state1D)
        @test all(hei .== 1.3)
        @test all(vel .== -0.1/1.3)
        @test state1D.height[1] ≈ 1.3
        @test state1D.vel[1] ≈ -0.1/1.3
        # 2D
        f[:,:,4] .= -0.2
        state.fout[:,:,4] .= -0.2
        Swalbe.moments!(height, velx, vely, f)
        Swalbe.moments!(state)
        @test all(height .== 1.1)
        @test all(velx .≈ 0.3/1.1)
        @test all(vely .== 0.2/1.1)
        @test all(state.height .== 1.1)
        @test all(state.velx .≈ 0.3/1.1)
        @test all(state.vely .== 0.2/1.1)
        f[:,:,5] .= -0.1
        state.fout[:,:,5] .= -0.1
        Swalbe.moments!(height, velx, vely, f)
        Swalbe.moments!(state)
        @test all(height .== 1.0)
        @test all(velx .≈ 0.3/1.0)
        @test all(vely .≈ 0.3/1.0)
        @test all(state.height .== 1.0)
        @test all(state.velx .≈ 0.3/1.0)
        @test all(state.vely .≈ 0.3/1.0)
        f[:,:,6] .= 0.1
        state.fout[:,:,6] .= 0.1
        Swalbe.moments!(height, velx, vely, f)
        Swalbe.moments!(state)
        @test all(height .== 1.1)
        @test all(velx .≈ 0.4/1.1)
        @test all(vely .≈ 0.4/1.1)
        @test all(state.height .== 1.1)
        @test all(state.velx .≈ 0.4/1.1)
        @test all(state.vely .≈ 0.4/1.1)
    end
end


text:
using Base.Test
using AutoEnvs

function test_index_ngsim_trajectory()
    filepath = Pkg.dir("NGSIM", "data", "trajdata_i101_trajectories-0750am-0805am.txt")
    index = index_ngsim_trajectory(filepath, verbose=1)
    @test length(index) > 2000
    for (k,v) in index
        @test in("ts", keys(v))
        @test in("te", keys(v))
    end
end

function test_sample_trajdata_vehicle()
    trajinfos = [Dict(1=>Dict("ts"=>1, "te"=>2))]
    traj_idx, vehid, ts, te = sample_trajdata_vehicle(trajinfos)
    @test traj_idx == 1
    @test vehid == 1
    @test ts >= 1
    @test te == 2
end

function test_sample_multiple_trajdata_vehicle_simple()

    # test that samples the valid ids
    trajinfos = [Dict(
        1=>Dict("ts"=>1, "te"=>100),
        2=>Dict("ts"=>1, "te"=>100),
        3=>Dict("ts"=>1, "te"=>100),
    )]
    offset = 50
    traj_idx, vehids, ts, te = sample_multiple_trajdata_vehicle(3, trajinfos, offset)
    @test sort!(vehids) == [1,2,3]

    # test that samples correct other id
    trajinfos = [Dict(
        1=>Dict("ts"=>1, "te"=>100),
        2=>Dict("ts"=>1, "te"=>100),
        3=>Dict("ts"=>70, "te"=>100),
    )]
    offset = 50
    traj_idx, vehids, ts, te = sample_multiple_trajdata_vehicle(2, trajinfos, offset, egoid=1, traj_idx=1)
    @test sort!(vehids) == [1,2]

    # test that it raises an error
    trajinfos = [Dict(
        1=>Dict("ts"=>1, "te"=>100),
        2=>Dict("ts"=>1, "te"=>100),
        3=>Dict("ts"=>1, "te"=>100),
    )]
    offset = 50
    @test_throws ErrorException sample_multiple_trajdata_vehicle(4, trajinfos, offset, verbose = false)

end

function test_sample_multiple_trajdata_vehicle_real()
    offset = 250
    n_veh = 100
    filepath = Pkg.dir("NGSIM", "data", "trajdata_i101_trajectories-0750am-0805am.txt")
    trajinfos = [index_ngsim_trajectory(filepath, minlength=offset, verbose=0)]
    traj_idx, vehids, ts, te = sample_multiple_trajdata_vehicle(n_veh, trajinfos, offset)
    @test length(vehids) == n_veh
    @test length(Set(vehids)) == n_veh
end

function test_build_feature_extractor()
    ext = build_feature_extractor()
end

function test_load_ngsim_trajdatas()
    filenames = [
        "trajdata_i80_trajectories-0400-0415.txt",
        "trajdata_i80_trajectories-0500-0515.txt",
        "trajdata_i80_trajectories-0515-0530.txt",
        "trajdata_i101_trajectories-0805am-0820am.txt",
        "trajdata_i101_trajectories-0820am-0835am.txt",
        "trajdata_i101_trajectories-0750am-0805am.txt"
    ]
    filepaths = [Pkg.dir("NGSIM", "data", fn) for fn in filenames]
    load_ngsim_trajdatas(filepaths, minlength=250)

end

@time test_index_ngsim_trajectory()
@time test_sample_trajdata_vehicle()
@time test_sample_multiple_trajdata_vehicle_simple()
@time test_sample_multiple_trajdata_vehicle_real()
@time test_build_feature_extractor()
@time test_load_ngsim_trajdatas()


text:
import Zygote
import JuMP

using Test: @test, @testset
using PartiallyObservedInverseGames.ForwardOptimalControl: solve_lqr
using PartiallyObservedInverseGames.InverseOptimalControl:
    solve_inverse_lqr, lqr_lagrangian, lqr_lagrangian_grad_x, lqr_lagrangian_grad_u
using LinearAlgebra: I

#=========================================== Forward LQR ===========================================#

T = 100
A = repeat([[
    1 1
    0 1
]], T)
B = repeat([[0, 1][:, :]], T)
Q = I
R = 100I
x0 = [10.0, 10.0]

forward_converged, forward_solution, forward_model = solve_lqr(A, B, Q, R, x0, T)

@testset "Forward LQR" begin
    @test forward_converged
end

#====================== Inverse LQR as nested constrained optimization problem =====================#

@testset "Inverse LQR" begin
    Q̃ = [
        [
            1//3 0
            0 0
        ],
        [
            0 0
            0 2//3
        ],
    ]
    R̃ = [R]
    inverse_converged, inverse_solution, inverse_model, Q_est, R_est =
        solve_inverse_lqr(forward_solution.x, Q̃, R̃; A, B)

    @testset "Gradient Check" begin
        grad_args = (inverse_solution.x, inverse_solution.u, inverse_solution.λ)
        grad_kwargs = (; Q = Q_est, R = R_est, A, B)

        ∇ₓL_ad, ∇ᵤL_ad = Zygote.gradient(
            (x, u) -> lqr_lagrangian(x, u, inverse_solution.λ; grad_kwargs...),
            inverse_solution.x,
            inverse_solution.u,
        )
        ∇ₓL_manual = lqr_lagrangian_grad_x(grad_args...; grad_kwargs...)
        ∇ᵤL_manual = lqr_lagrangian_grad_u(grad_args...; grad_kwargs...)
        atol = 1e-10

        @test isapprox(∇ₓL_ad, ∇ₓL_manual; atol = atol)
        @test isapprox(∇ₓL_ad[:, 2:end], inverse_solution.∇ₓL; atol = atol)
        @test isapprox(∇ᵤL_ad, ∇ᵤL_manual; atol = atol)
        @test isapprox(∇ᵤL_ad, inverse_solution.∇ᵤL; atol = atol)
    end

    @testset "Solution Sanity" begin
        @test inverse_converged
        @test Q_est[1, 1] / Q_est[2, 2] ≈ Q[1, 1] / Q[2, 2]
        @test Q_est[1, 1] / R_est[1, 1] ≈ Q[1, 1] / R[1, 1]
        @test Q_est[2, 2] / R_est[1, 1] ≈ Q[2, 2] / R[1, 1]
        @test inverse_solution.x ≈ solve_lqr(A, B, Q_est, R_est, x0, T)[2].x
    end
end


text:
using Logging
using Printf
using Statistics
using NCDatasets
using CUDA

using Oceananigans
using Oceananigans.Grids
using Oceananigans.Advection
using Oceananigans.OutputWriters
using Oceananigans.Utils
using JULES

using Oceananigans.Fields: cpudata

Logging.global_logger(OceananigansLogger())

const km = kilometers
const hPa = 100.0

#####
##### Hydrostatic base state
#####

pₛ = 1000hPa
Tₛ = 300.0

θ₀(x, y, z) = Tₛ
p₀(x, y, z) = pₛ * (1 - g*z / (cₚ*Tₛ))^(cₚ/R)
T₀(x, y, z) = Tₛ * (p₀(x, y, z)/pₛ)^(R/cₚ)
ρ₀(x, y, z) = p₀(x, y, z) / (R*T₀(x, y, z))

#####
##### Model setup
#####

L = 2.56km
Δ = 40meters
N = Int(L/Δ)

topo = (Periodic, Periodic, Bounded)
domain = (x=(0, L), y=(0, L), z=(0, L))
grid = RegularCartesianGrid(topology=topo, size=(N, N, N), halo=(3, 3, 3); domain...)

gas = DryEarth()

# Sponge layer at the top with exponential decay length scale γ.
@inline damping_profile(z, L, γ) = exp(- (L - z) / γ)

forcing_params = (λ = 1/(10seconds), L=grid.Lz, γ=0.2grid.Lz)

@inline rayleigh_damping_ρu(i, j, k, grid, clock, fields, p) =
    @inbounds - p.λ * damping_profile(grid.zC[k], p.L, p.γ) * fields.ρu[i, j, k]

@inline rayleigh_damping_ρv(i, j, k, grid, clock, fields, p) =
    @inbounds - p.λ * damping_profile(grid.zC[k], p.L, p.γ) * fields.ρv[i, j, k]

@inline rayleigh_damping_ρw(i, j, k, grid, clock, fields, p) =
    @inbounds - p.λ * damping_profile(grid.zC[k], p.L, p.γ) * fields.ρw[i, j, k]

ρu_forcing = Forcing(rayleigh_damping_ρu, discrete_form=true, parameters=forcing_params)
ρv_forcing = Forcing(rayleigh_damping_ρv, discrete_form=true, parameters=forcing_params)
ρw_forcing = Forcing(rayleigh_damping_ρw, discrete_form=true, parameters=forcing_params)

# Exponentially decaying radiative forcing
@inline radiative_profile(z, Q₀, ℓ) = Q₀ * exp(-z/ℓ)

radiative_params = (Q₀=100.0, ℓ=100.0)

@inline radiative_forcing_ρe(i, j, k, grid, clock, fields, p) =
    @inbounds (radiative_profile(grid.zC[k], p.Q₀, p.ℓ) - radiative_profile(grid.zC[k+1], p.Q₀, p.ℓ)) / grid.Δz

ρe_forcing = Forcing(radiative_forcing_ρe, discrete_form=true, parameters=radiative_params)

model = CompressibleModel(
              architecture = GPU(),
                      grid = grid,
                 advection = WENO5(),
                     gases = gas,
    thermodynamic_variable = Energy(),
                   closure = IsotropicDiffusivity(ν=1e-2, κ=1e-2),
                   forcing = (ρu=ρu_forcing, ρv=ρv_forcing, ρw=ρw_forcing, ρe=ρe_forcing)
)

g  = model.gravity
R, cₚ, cᵥ = gas.ρ.R, gas.ρ.cₚ, gas.ρ.cᵥ
uᵣ, Tᵣ, ρᵣ, sᵣ = gas.ρ.u₀, gas.ρ.T₀, gas.ρ.ρ₀, gas.ρ.s₀
ρe₀(x, y, z) = ρ₀(x, y, z) * (uᵣ + cᵥ * (T₀(x, y, z) - Tᵣ) + g*z) + randn()

# Set initial state
set!(model.tracers.ρ, ρ₀)
set!(model.tracers.ρe, ρe₀)
update_total_density!(model)

function print_progress(simulation)
    model, Δt = simulation.model, simulation.Δt
    tvar = model.thermodynamic_variable

    progress = 100 * model.clock.time / simulation.stop_time
    message = @sprintf("[%05.2f%%] iteration = %d, time = %s, CFL = %.4e, acoustic CFL = %.4e",
                       progress, model.clock.iteration, prettytime(model.clock.time), cfl(model, Δt),
                       acoustic_cfl(model, Δt))

    @info message

    return nothing
end

simulation = Simulation(model, Δt=0.02second, stop_time=1hour, iteration_interval=20,
                        progress=print_progress)

fields = Dict(
    "ρ"  => model.total_density,
    "ρu" => model.momenta.ρu,
    "ρv" => model.momenta.ρv,
    "ρw" => model.momenta.ρw,
    "ρe" => model.tracers.ρe
)

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, fields, filepath="dry_convection.nc", time_interval=10seconds)

# Save base state to NetCDF.
ds = simulation.output_writers[:fields].dataset
ds_ρ = defVar(ds, "ρ₀", Float32, ("xC", "yC", "zC"))
ds_ρe = defVar(ds, "ρe₀", Float32, ("xC", "yC", "zC"))

x, y, z = nodes((Cell, Cell, Cell), grid, reshape=true)
ds_ρ[:, :, :] = ρ₀.(x, y, z)
ds_ρe[:, :, :] = ρe₀.(x, y, z)

run!(simulation)


text:
function read_answers( filename )

    answers = Vector{String}[]
    open(filename) do f
        batch = String[]
        for line in readlines(f)
            if isempty(line)
                push!(answers, batch)
                batch = String[]
            else
                push!(batch, line)
            end
        end
        push!(answers, batch)
    end
    return answers

end

function counts_1( answers )

    return sum(length.(unique.(join.(answers))))

end

function counts_2( answers )

    res = 0
    for group in answers, c in 'a':'z'
        res += all([ c in ans for ans in group])
    end
    return res

end

answers = read_answers( "input06.txt" )
println(counts_1( answers ))
println(counts_2( answers ))

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
using Documenter, RandomMatrix, Plots, LinearAlgebra

ENV["GKSwstype"] = "100"
gr()

const PAGES = Any["Home" =>"index.md",
      "Random Matrix Models" =>[
            "I.I.D. Matrices"=>"IID Matrices.md",
            "Correlated Matrices"=>"Correlated Matrices.md",
            "Unitary Matrices"=>"Unitary.md",
            "Stochastic Matrices"=>"Stochastic.md",
            "Patterned Matrices "=>"Patterned.md",
      ],
      "Distributions"=>[
            "Scarlar-Valued Distributions"=> "ScalarValued.md",
            "Matrix-Valued Distributions"=>"MatrixValued.md"
            ],
      "Matrix and Spectral Transformations"=>[
            "Spectral Spacing"=>"Spectral Spacing.md",
            "Matrix Transformations"=>"matrixtransformation.md"
            ],

      "Under Construction" =>[
            "Under Construction1" =>"Others.md",
              "Under Construction2" =>"Gallery.md"
      ]

]

makedocs(modules = [RandomMatrix],
      sitename="RandomMatrix.jl",
      authors = "Wei Yang",
      pages = PAGES
      )

deploydocs(repo = "github.com/weiyang2048/RandomMatrix.jl.git",
devbranch = "main",
push_preview=true )


text:
module OSXNotifier

export notify

@osx_only begin
	using Homebrew

	if Homebrew.installed("terminal-notifier")
		terminalnotifier = joinpath(Homebrew.prefix("terminal-notifier"), "bin", "terminal-notifier")
		if !isfile(terminalnotifier)
			error("terminal-notifier binary does not exist: $terminalnotifier")
		end
	else
		error("OSXNotifier not properly installed. Please run Pkg.build(\"OSXNotifier\")")
	end

	function notify(message=""; title="Julia", subtitle="", open="", group="", sound=false, sender="org.julialang")
		if group != "" && sound != "" && sound != false
			run(`$terminalnotifier -sender $sender -message $message -title $title -subtitle $subtitle -open $open -group $group -sound $sound`)
		elseif sound != "" && sound != false
			run(`$terminalnotifier -sender $sender -message $message -title $title -subtitle $subtitle -open $open -sound $sound`)
		elseif group != ""
			run(`$terminalnotifier -sender $sender -message $message -title $title -subtitle $subtitle -open $open -group $group`)
		else
			run(`$terminalnotifier -sender $sender -message $message -title $title -subtitle $subtitle -open $open`)
		end
	end

	function remove(group="ALL"; sender="org.julialang")
		run(`$terminalnotifier -remove $group -sender $sender`)
	end

	#function list(group="ALL")
	#	lines = readall(`$terminalnotifier -list $group`)
	#	lines = split(lines, "\n")
	#	lines = map(line -> split(line, "\t"), lines)
	#	lines = map(line -> map(element -> convert(ASCIIString, element), line), lines)
	#	keys = lines[1]
	#	lines = lines[2:end-1]
	#	map(line -> Dict(zip(keys, line)), lines)
	#end
end

end # module


text:
# Not yet implemented
# Placeholder to shade samples of regression lines


text:
# test padding
@testset "padding dimension" begin
    @testset "padding 1D" begin
        # no meta dims
        input = ones(Float32, 10)
        padd = pad(input, (5,))
        padd = pad(input, 5)
        @test size(padd) == (20,)
        @test padd[1] == 0.0
        @test padd[6] == 1.0
        @test eltype(padd) == eltype(input)

        # one meta dim, different type
        input = ones(Float64, 10, 2)
        padd = pad(input, 5)
        @test size(padd) == (20, 2)
        @test minimum(padd[1,:] .== 0.0)
        @test minimum(padd[6,:] .== 1.0)
        @test eltype(padd) == eltype(input)

        # two meta dims, different type
        input = ones(ComplexF64, 10, 5, 3)
        padd = pad(input, 5)
        @test size(padd) == (20, 5, 3)
        @test minimum(padd[1,:,:] .== 0.0)
        @test minimum(padd[6,:,:] .== 1.0)
        @test eltype(padd) == eltype(input)

        # test the gradient is just an indicator with the right offset
        i = 6;
        grad = gradient(x -> pad(x, 5)[i,1,1], randn(5, 5, 2))[1]
        @test size(grad) == (5, 5, 2)
        @test grad[1,1,1] == 1
    end

    @testset "padding 2D" begin
        # no meta dims
        input = ones(Float32, 10, 10)
        padd = pad(input, (5, 5))
        @test size(padd) == (20, 20)
        @test padd[1,1] == 0.0
        @test padd[6,6] == 1.0
        @test eltype(padd) == eltype(input)

        # one meta dim, different type
        input = ones(Float64, 10, 10, 2)
        padd = pad(input, (5, 5))
        @test size(padd) == (20, 20, 2)
        @test minimum(padd[1,1,:] .== 0.0)
        @test minimum(padd[6,6,:] .== 1.0)
        @test eltype(padd) == eltype(input)

        # two meta dims, different type
        input = ones(ComplexF64, 10, 10, 5, 3)
        padd = pad(input, (5, 5))
        @test size(padd) == (20, 20, 5, 3)
        @test minimum(padd[1,1,:,:] .== 0.0)
        @test minimum(padd[6,6,:,:] .== 1.0)
        @test eltype(padd) == eltype(input)

        # test the gradient is just an indicator with the right offset
        i = 6; j = 7
        grad = gradient(x -> pad(x, (5, 6))[i,j,1,1], randn(5, 5, 5, 2))[1]
        @test size(grad) == (5, 5, 5, 2)
        @test grad[1,1,1,1] == 1
    end
end

@testset "boundaries" begin
    x = randn(10, 10, 10)
    bc = Pad((10, 10))
    xbc, usedInds = FourierFilterFlux.applyBC(x, bc, 2)
    @test size(xbc) == (10 + 2 * 10, 10 + 2 * 10, 10)
    @test usedInds == (11:20, 11:20)

    bc = Pad(10)
    xbc, usedInds = FourierFilterFlux.applyBC(x, bc, 1)
    @test size(xbc) == (10 + 2 * 10, 10, 10)
    @test usedInds == (11:20,)

    bc = FourierFilterFlux.Periodic()
    xbc, usedInds = FourierFilterFlux.applyBC(x, bc, 1)
    @test size(xbc) == (10, 10, 10)
    @test usedInds == (1:10,)

    bc = FourierFilterFlux.Sym();
    xbc, usedInds = FourierFilterFlux.applyBC(x, bc, 2);
    @test size(xbc) == (20, 20, 10)
    @test usedInds == (1:10, 1:10)

    eSz, nBound = FourierFilterFlux.effectiveSize((10, 10), Pad((-1, -1)));
    @test eSz == (16, 16)
    @test nBound == Pad(3, 3)
end


text:
#-----------------------------------------------------------------------
module cereal      # cereal module
#-----------------------------------------------------------------------

using LinearAlgebra

const RealVec{T<:Real} = Array{T,1}      # Defining vector datatype
const RealMtx{T<:Real} = Array{T,2}      # Defining matrix datatype

#-----------------------------------------------------------------------
function η( V1::RealVec , V2::RealVec )  # Computes Minkowski product
    nv1 = length(V1)
    nv2 = length(V2)

    met = -V1[1]*V2[1]

    if nv1==nv2
        for i=2:nv1
            met += V1[i]*V2[i]
        end
        return met
    else
        print("Vectors are of a different dimension.")
    end
end  # End η

#-----------------------------------------------------------------------
function Frame( X::RealMtx , erm::Bool=true )
    # Constructs spatial frame
    tpfl=typeof(X[1,1])       # Extracting floating point datatype
    s = 1       # Integer for constructing frame
    b = true    # Boolean variable - stops routine if not spacelike

    E = zeros(tpfl,4,3)

    for i=1:4
        for j=i+1:4
            V = X[:,i] - X[:,j]       # Computing difference between pts
            b = η( V , V ) > zero(tpfl)
            if b
                if s<=3
                    E[:,s]=V          # Frame vector
                end
                s += 1
            end
        end
    end

    if b
        return E
    else
        if erm
            println("Points not spacelike separated")
        end
        return zeros(tpfl,4,3)  # Returns zeros in case of error
    end
end  # End Frame

#-----------------------------------------------------------------------
function NormVecF( E::RealMtx , erm::Bool=true )
    # Constructs timelike normal to frame
    tpfl=typeof(E[1,1])
        NVr = [ E[2,3]*E[3,2]*E[4,1] - E[2,2]*E[3,3]*E[4,1] -
                E[2,3]*E[3,1]*E[4,2] + E[2,1]*E[3,3]*E[4,2] +
                E[2,2]*E[3,1]*E[4,3] - E[2,1]*E[3,2]*E[4,3]   ;
                #
                E[1,3]*E[3,2]*E[4,1] - E[1,2]*E[3,3]*E[4,1] -
                E[1,3]*E[3,1]*E[4,2] + E[1,1]*E[3,3]*E[4,2] +
                E[1,2]*E[3,1]*E[4,3] - E[1,1]*E[3,2]*E[4,3]   ;
                #
                E[1,1]*E[2,2]*E[4,3] - E[1,3]*E[2,2]*E[4,1] +
                E[1,2]*E[2,3]*E[4,1] + E[1,3]*E[2,1]*E[4,2] -
                E[1,1]*E[2,3]*E[4,2] - E[1,2]*E[2,1]*E[4,3]   ;
                #
                E[1,3]*E[2,2]*E[3,1] - E[1,2]*E[2,3]*E[3,1] -
                E[1,3]*E[2,1]*E[3,2] + E[1,1]*E[2,3]*E[3,2] +
                E[1,2]*E[2,1]*E[3,3] - E[1,1]*E[2,2]*E[3,3]   ]
        norm = η( NVr , NVr )
        if norm < zero(tpfl)
            if NVr[1] < zero(tpfl)
                NV = -NVr / sqrt(abs(norm))
            elseif NVr[1] > zero(tpfl)
                NV = NVr / sqrt(abs(norm))
            else
                if erm
                println("Time component of normal vector too small.")
                end
                NV = zeros(tpfl,4)  # Returns zeros in case of error
            end
        else
            if erm
                println("Normal vector not timelike.")
            end
            NV = zeros(tpfl,4)      # Returns zeros in case of error
        end
    return NV
end  # End NormVecF

#-----------------------------------------------------------------------
function LTM( NV::RealVec )
    # Constructs Lorentz transformation matrix
    tpfl=typeof(NV[1])
        l = one(tpfl)           # l has a value of one
        γ = NV[1]               # γ is the Lorentz factor
        δ = (γ - l)             # A useful quantity
        NVs = NV[2]^2 + NV[3]^2 + NV[4]^2      # Square of spatial part
        if NVs > 0              # Normalization factor check
            nf = 1/sqrt(NVs)    # Normalization factor
            β = sqrt(NVs/(γ*γ)) # β is v/c
        else
            nf = zero(tpfl)     # Set normalization factor to zero
            β = zero(tpfl)      # β is v/c
        end
        vx = nf*NV[2]     # x-component of unit vector
        vy = nf*NV[3]     # y-component of unit vector
        vz = nf*NV[4]     # z-component of unit vector
    return [  γ        -γ*β*vx       -γ*β*vy       -γ*β*vz      ;
             -γ*β*vx   l + δ*(vx^2)  δ*vx*vy       δ*vx*vz      ;
             -γ*β*vy   δ*vy*vx       l + δ*(vy^2)  δ*vy*vz      ;
             -γ*β*vz   δ*vz*vx       δ*vz*vy       l + δ*(vz^2) ]
end  # End LTM

#-----------------------------------------------------------------------
function IPfinder( Y::RealMtx )
    # Finds intersection of light cones in adapted frame
    tpfl=typeof(Y[1,1])

    x = zeros(tpfl,3,4)     # Create containers
    v = zeros(tpfl,3,3)

    x[1,:] = Y[2,:]    # Spatial points for corners of tetrahedron
    x[2,:] = Y[3,:]
    x[3,:] = Y[4,:]

    v[:,1] = x[:,2] - x[:,1]    # "Frame" vectors centered on x[:,1]
    v[:,2] = x[:,3] - x[:,1]
    v[:,3] = x[:,4] - x[:,1]

    B = [ x[:,2]⋅x[:,2] - x[:,1]⋅x[:,1] ; x[:,3]⋅x[:,3] - x[:,1]⋅x[:,1];
          x[:,4]⋅x[:,4] - x[:,1]⋅x[:,1] ] / 2    # Constructing B vector

    A = transpose([ v[:,1]  v[:,2]  v[:,3] ])    # Constructing A matrix

    xc = inv(A)*B       # Compute circumcenter
    rc = sqrt( (xc-x[:,1])⋅(xc-x[:,1]) )    # Distance to circumcenter
    tc = rc + Y[1,1]    # Compute time coordinate
    return [ tc ; xc[1] ; xc[2] ; xc[3] ]
end  # End IPfinder

#-----------------------------------------------------------------------
function locator( X::RealMtx , erm::Bool=true , erc::Bool=false )
    # Special relativistic locator code
    tpfl=typeof(X[1,1])

    Y = zeros(tpfl,4,4)         # Create containers
    Z = zeros(tpfl,4)

    ERR = false

    E = Frame(X,erm)           # Constructing spatial frame
    NV = NormVecF(E,erm)       # Constructing normal vector to frame
    if dot(NV,NV)>0 && ERR==false
        Λ = LTM(NV)             # Lorentz transformation matrix

        for i=1:4
            Y[:,i] = Λ * X[:,i] # Lorentz transformation
        end

        XP = IPfinder(Y)        # Intersection point in adapted frame
        ΛR = inv(Λ)             # Lorentz transformation matrix

        for i=1:4
            Z = ΛR * XP         # Lorentz transformation
        end
        if erc
            return (Z,true)
        else
            return Z
        end
    else
        if erm
            println("Errors encountered, returning zero vector.")
        end
        if erc
            return (zeros(tpfl,4),false)
        else
            return zeros(tpfl,4)    # Returns zero if errors encountered
        end
    end
end  # End locator

#-----------------------------------------------------------------------
end     # End scope of module cereal
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
module cerealtest    # Test module for cereal
#-----------------------------------------------------------------------

using LinearAlgebra
import ..cereal         # Importing functions from cereal

const RealVec{T<:Real} = Array{T,1}      # Defining vector datatype
const RealMtx{T<:Real} = Array{T,2}      # Defining matrix datatype

#-----------------------------------------------------------------------
function vgenerator( tpfl::DataType )    # Generates random 3-vector
    θ = pi*rand(tpfl)
    ϕ = 2*pi*rand(tpfl)
    V = [   sin(θ)*cos(ϕ)   ;
            sin(θ)*sin(ϕ)   ;
            cos(θ)          ]
    return V
end  # End vgenerator

#-----------------------------------------------------------------------
function nvgenerator( tpfl::DataType )   # Generates random 4-velocity
    v = vgenerator(tpfl)
    v = rand(tpfl)*rand(tpfl)*v
    V = [ one(tpfl) ; v[1] ; v[2] ; v[3] ]
    fl = one(tpfl) / sqrt(abs(cereal.η(V,V)))

    NV = fl*V
    return NV
end  # End nvgenerator

#-----------------------------------------------------------------------
function TVchk( Y::RealMtx , δV::Real=1e-14 )
    # Check tetrahedron volume is > |δV|    tpfl=typeof(Y[1,1])
    tpfl=typeof(Y[1,1])

    x = zeros(tpfl,3,4)     # Create containers
    v = zeros(tpfl,3,3)

    x[1,:] = Y[2,:]         # Spatial points for corners of tetrahedron
    x[2,:] = Y[3,:]
    x[3,:] = Y[4,:]

    v[:,1] = x[:,2] - x[:,1]    # "Frame" vectors centered on x[:,1]
    v[:,2] = x[:,3] - x[:,1]
    v[:,3] = x[:,4] - x[:,1]

    Vol = abs(det([ v[:,1]  v[:,2]  v[:,3] ]))/6    # Volume

    if Vol > abs(δV)
        return true
    else
        return false
    end
end  # End TVchk

#-----------------------------------------------------------------------
function pointgenerator( tpfl::DataType , δV::Real=1e-14 )
    # Generates four random points W all with same time coordinate
    b = false

    W = zeros(tpfl,4,4)

    while b == false
        V1 = (one(tpfl)+rand(tpfl))*vgenerator(tpfl)/(2*one(tpfl))
        V2 = (one(tpfl)+rand(tpfl))*vgenerator(tpfl)/(2*one(tpfl))
        V3 = (one(tpfl)+rand(tpfl))*vgenerator(tpfl)/(2*one(tpfl))
        V4 = (one(tpfl)+rand(tpfl))*vgenerator(tpfl)/(2*one(tpfl))

        Vt = rand(tpfl)

        W = [   Vt     Vt     Vt     Vt     ;
                V1[1]  V2[1]  V3[1]  V4[1]  ;
                V1[2]  V2[2]  V3[2]  V4[2]  ;
                V1[3]  V2[3]  V3[3]  V4[3]  ]

        b = TVchk(W,δV)
    end

    return W
end  # End pointgenerator

#-----------------------------------------------------------------------
function LT( W::RealMtx, NV::RealVec )
    # Lorentz transforms W to frame where vector NV is in t direction
    tpfl=typeof(W[1,1])
    l = one(tpfl)           # l has a value of one
    Z = zeros(tpfl,4,4)

    Λ = cereal.LTM(NV)
    β = rand(tpfl)/(tpfl(2))
    γ = l/sqrt( l - β^2 )
    δ = ( γ - l )

    vx = β*NV[1]            # components of unit vector
    vy = β*NV[2]
    vz = β*NV[3]

    Λ = [  γ        -γ*β*vx        -γ*β*vy       -γ*β*vz       ;
           -γ*β*vx   l + δ*(vx^2)  δ*vx*vy       δ*vx*vz       ;
           -γ*β*vy   δ*vy*vx       l + δ*(vy^2)  δ*vy*vz       ;
           -γ*β*vz   δ*vz*vx       δ*vz*vy       l + δ*(vz^2)  ]

    for i=1:4
        Z[:,i] = Λ * W[:,i]         # Lorentz transformation
    end
    return Z
end  # End LT

#-----------------------------------------------------------------------
function epgen( tpfl::DataType=Float64 , δV::Real=1e-14 )
    # Randomly generates emission pts
    W  = pointgenerator(tpfl,δV)
    NV = nvgenerator(tpfl)
    return LT(W,NV)
end  # End epgen

#-----------------------------------------------------------------------
function nsqr( V1::RealVec )        # Squared norm ratio
    return abs(cereal.η(V1,V1)/(V1⋅V1))
end  # End nsqr

#-----------------------------------------------------------------------
function single( q::Real , P::RealVec , X::RealMtx )
    # Checks if separation vectors are null
    tpfl=typeof(q)
    Y = copy(P)

    VT = [ nsqr( Y - X[:,1] ) ; nsqr( Y - X[:,2] ) ;
           nsqr( Y - X[:,3] ) ; nsqr( Y - X[:,4] ) ]

    av = ( VT[1] + VT[2] + VT[3] + VT[4] ) / (4*one(tpfl))

    if av < q
        return true
    else
        return ( X , av )
    end
end  # End single

#-----------------------------------------------------------------------
function full( iters::Number , q::Real , δV::Real=1e-14
               , erm::Bool=true , counter::Bool=true )
    # Main test function--q determines floating point datatype
    tpfl=typeof(q)
    lb=false
    mi = zero(Int64)

    for i=1:Int64(iters)
        B = true
        if counter
            print("\r$i")
        end
        while B
            X = epgen(tpfl,δV)
            P = cereal.locator(X,erm,true)
            b = single(q,P[1],X)
            if b != true
                print("\n", b[2],"\n")
                lb = true
                mi += 1
            end
            if P[2]
                B = false
            end
        end
    end
    if lb
        print("\rTest ended with ",mi," failed cases out of "
              ,iters," \n")
    else
        print("\rTest ended with zero failed cases out of ",iters," \n")
    end
end  # End full

#-----------------------------------------------------------------------
end     # End scope of module cerealtest
#-----------------------------------------------------------------------


text:
import Dolo

path = Dolo.pkg_path

@testset "testing discretizing functions" begin

    rho_z = 0.9
    rho_b = 0.9
    sig_z = 0.001
    sig_b = 0.001

    proc = Dolo.ProductProcess(
        Dolo.VAR1(rho_z, ones(1,1)*sig_z^2),
        Dolo.VAR1(rho_b, ones(1,1)*sig_b^2)
    )

    Dolo.discretize(proc)
    exo = Dolo.discretize(Dolo.GDP, proc)

    @test true

end


text:
using Libfreenect2
using Base.Test

import Libfreenect2: width, height, bytes_per_pixel, exposure, gain, gamma

@testset "Frame" begin
    frame = FramePtr(512, 424, 4)
    @test width(frame) == 512
    @test height(frame) == 424
    @test bytes_per_pixel(frame) == 4
    @test exposure(frame) == 0.0
    @test gain(frame) == 0.0
    @test gamma(frame) == 0.0
end

@testset "libfreenect2 basics" begin
    f = Freenect2()
    num_devices = enumerateDevices(f)
    @test num_devices > 0

    serial = getDefaultDeviceSerialNumber(f)
    @test !isempty(serial)

    device = openDevice(f, serial)

    firmware_version = getFirmwareVersion(device)
    @test !isempty(firmware_version)

    listener = SyncMultiFrameListenerPtr()
    setIrAndDepthFrameListener(device, listener)
    setColorFrameListener(device, listener)

    start(device)

    # NOTE: must be called after start(device)
    registration = Registration(getIrCameraParams(device),
        getColorCameraParams(device))
    # these are of the same memory layout to color stream
    undistorted = FramePtr(512, 424, 4, key=Libfreenect2.FRAME_COLOR)
    registered = FramePtr(512, 424, 4, key=Libfreenect2.FRAME_COLOR)

    frames = waitForNewFrame(listener)

    stop(device)

    color = frames[FrameType.COLOR]
    ir = frames[FrameType.IR]
    depth = frames[FrameType.DEPTH]

    apply(registration, color, depth, undistorted, registered)

    for frame in [color, ir, depth]
        @test bytes_per_pixel(frame) == 4
    end

    @test width(color) == 1920
    @test height(color) == 1080
    @test width(ir) == 512
    @test height(ir) == 424
    @test width(depth) == 512
    @test height(depth) == 424

    @testset "Convert frame to Julia array" begin
        colorarr = convert(Array{UInt8, 3}, color)
        @test size(colorarr) == (4, 1920, 1080)
        irarr = convert(Array{Float32,2}, ir)
        @test size(irarr) == (512, 424)
        deptharr = convert(Array{Float32,2}, depth)
        @test size(deptharr) == size(irarr)
    end

    release(listener, frames)
    close(device)
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libcsfml-graphics", "csfml-graphics-2"], :libcsfml_graphics),
    LibraryProduct(prefix, ["libsfml-window", "sfml-window-2"], :libsfml_window),
    LibraryProduct(prefix, ["libsfml-audio", "sfml-audio-2"], :libsfml_audio),
    LibraryProduct(prefix, ["libsfml-network", "sfml-network-2"], :libsfml_network),
    LibraryProduct(prefix, ["libsfml-system", "sfml-system-2"], :libsfml_system),
    LibraryProduct(prefix, ["libsfml-graphics", "sfml-graphics-2"], :libsfml_graphics),
    LibraryProduct(prefix, ["libcsfml-system", "csfml-system-2"], :libcsfml_system),
    LibraryProduct(prefix, ["libcsfml-network", "csfml-network-2"], :libcsfml_network),
    LibraryProduct(prefix, ["libcsfml-window", "csfml-window-2"], :libcsfml_window),
    LibraryProduct(prefix, ["libcsfml-audio", "csfml-audio-2"], :libcsfml_audio),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/Gnimuc/SFMLBuilder/releases/download/v2.5-4"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Windows(:i686) => ("$bin_prefix/SFML.v2.5.1.i686-w64-mingw32-gcc7.tar.gz", "39f8acac8c4c9e63cc6d649cfda24f5de9e68a7871abcf87c99f0791d9fc54cb"),
    MacOS(:x86_64) => ("$bin_prefix/SFML.v2.5.1.x86_64-apple-darwin14.tar.gz", "33a64920293273ca616ddf795460fe9fdf0a403cab67f8d231deaa3d30c4f99c"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/SFML.v2.5.1.x86_64-linux-gnu-gcc7.tar.gz", "b6ef09c057426c945644df018dd7614e2a9da9920f30ab82102fcb39b339cc6e"),
    Windows(:x86_64) => ("$bin_prefix/SFML.v2.5.1.x86_64-w64-mingw32-gcc7.tar.gz", "eceb442ca35057feacab099eaf614029033d4da958d8acf579de94c83dd38c25"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose, ignore_platform=true)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
module SpikeNetOpt

    using PyCall
    using OrderedCollections
    using LinearAlgebra
    using SpikeSynchrony
    using SpikingNeuralNetworks
    SNN = SpikingNeuralNetworks
    SNN.@load_units
    using Evolutionary#, Test, Random
    import DataStructures
    using JLD
    #using Reexport
    using UnicodePlots
    using Statistics
    using JLD
    using Distributed
    using SharedArrays
    using Plots
    using Evolutionary
    using Distributions
    using LightGraphs
    using Metaheuristics
    using ClearStacktrace
    using SharedArrays
    #using WaspNet
    #using SpikingNN
    using SignalAnalysis
    #using Reexport
    #using Spikes
    # Export Function names and struct names found in file utils
    #=
    export get_trains
    export get_vm
    export sim_net_darsnack

    function get_trains end
    function get_vm end
    function checkmodel end
    function make_net_from_graph_structure end
    function make_net end
    function spike_train_difference end
    function get_spikes end
    function get_trains end
    function sim_net_darsnack end
    =#
    include("utils.jl")
    #include("current_search.jl")
    include("sdo_network.jl")

    #include("sciunit.jl")

end


text:
function line2segments(points)
    npoints = length(points)
    indices = RPR.rpr_int[]
    count = 0
    for i in 1:npoints
        push!(indices, i -  1)
        count += 1
        if count == 4 && !(i == npoints)
            push!(indices, i - 1)
            count = 1
        end
    end
    missing_points = 4 - (length(indices) % 4)
    if missing_points != 4
        append!(indices, fill(last(indices), missing_points))
    end
    return indices
end

function to_rpr_object(context, matsys, scene, plot::Makie.Lines)
    points = decompose(Point3f, to_value(plot[1]))
    isempty(points) && return nothing
    npoints = length(points)
    indices = line2segments(points)
    radius = [plot.linewidth[] / 1000]
    curve = RPR.Curve(context, points, indices, radius, [Vec2f(0.0)], [length(indices) ÷ 4])
    material = RPR.MaterialNode(matsys, RPR.RPR_MATERIAL_NODE_DIFFUSE)
    set!(material, RPR.RPR_MATERIAL_INPUT_COLOR, to_color(plot.color[]))
    set!(curve, material)
    return curve
end

function to_rpr_object(context, matsys, scene, plot::Makie.LineSegments)
    points = decompose(Point3f, to_value(plot[1]))
    segments = TupleView{2,2}(RPR.rpr_int(0):RPR.rpr_int(length(points) - 1))
    indices = RPR.rpr_int[]

    for (a, b) in segments
        push!(indices, a, a, b, b)
    end

    nsegments = length(indices) ÷ 4

    radius = if plot.linewidth[] isa AbstractVector
        lw = plot.linewidth[]
        perseg = if length(lw) == length(points)
            lw[1:2:end]
        elseif length(lw) == nsegments
            lw
        else
            error("Length $(nsegments) doesn't match $(length(lw))")
        end
        Float32.(lw ./ 1000)
    else
        fill(Float32(plot.linewidth[] / 1000), nsegments)
    end

    curve = RPR.Curve(context, points, indices, radius, Vec2f.(0.0, LinRange(0, 1, nsegments)),
                      fill(1, nsegments))
    material = RPR.DiffuseMaterial(matsys)
    color = to_color(plot.color[])

    function set_color!(colorvec)
        tex = RPR.ImageTextureMaterial(matsys)
        ncols = length(colorvec)
        img = RPR.Image(context, reshape(colorvec, (ncols, 1)))
        tex.data = img
        material.color = tex
        return
    end

    if color isa AbstractVector{<:Colorant}
        set_color!(copy(color))
    elseif color isa AbstractVector{<:Number}
        sampler = Makie.sampler(to_colormap(plot.colormap[]), color; scaling=Makie.Scaling(identity, plot.colorrange[]))
        set_color!(collect(sampler))
    else
        material.color = to_color(color)
    end
    set!(curve, material.node)
    return curve
end


text:
# Buffers/strings

abstract type AbstractBuffer end
abstract type AbstractSizedBuffer end
abstract type AbstractUnsizedBuffer end

const BufPtr = Core.LLVMPtr{UInt8,0}

@inline @generated function create_buffer(::Val{N}) where N
    Context() do ctx
        T_i8 = LLVM.Int8Type(ctx)
        T_pi8 = LLVM.PointerType(T_i8)
        T_i64 = LLVM.Int64Type(ctx)
        T_buf = LLVM.ArrayType(T_i8, N)
        llvm_f, _ = create_function(T_pi8)
        mod = LLVM.parent(llvm_f)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            # Allocate stack buffer
            malloc_ft = LLVM.FunctionType(T_pi8, [T_i64])
            malloc_f = LLVM.Function(mod, "malloc", malloc_ft)
            buf_len = ConstantInt(T_i64, N)
            buf_ptr = call!(builder, malloc_f, [buf_len])

            # Zero-out buffer
            zero_value = ConstantInt(T_i8, 0)
            _memset!(builder, ctx, mod, buf_ptr, zero_value, buf_len, ConstantInt(LLVM.Int1Type(ctx), 0))

            ret!(builder, buf_ptr)
        end
        call_function(llvm_f, BufPtr)
    end
end

"""
    SizedBuffer <: AbstractSizedBuffer

Represents a buffer/string with a known size.
"""
struct SizedBuffer <: AbstractSizedBuffer
    ptr::BufPtr
    length::UInt32
end
SizedBuffer(buf::AbstractSizedBuffer) = SizedBuffer(pointer(buf), length(buf))
SizedBuffer(ptr::Ptr, length) = SizedBuffer(reinterpret(BufPtr, ptr), length)
Base.pointer(str::SizedBuffer) = str.ptr
Base.length(str::SizedBuffer) = str.length
@inline create_buffer(N::Int) = SizedBuffer(create_buffer(Val(N)), N)
macro create_string(str::String)
    N = length(str)+1
    :(SizedBuffer(create_string($(Val(Symbol(str)))), $N))
end

@inline @generated function create_string(::Val{str}) where str
    Context() do ctx
        _str = String(str)
        T_i8 = LLVM.Int8Type(ctx)
        T_pi8 = LLVM.PointerType(T_i8)
        T_i64 = LLVM.Int64Type(ctx)
        T_i1 = LLVM.Int1Type(ctx)
        T_buf = LLVM.ArrayType(T_i8, length(_str)+1)
        llvm_f, _ = create_function(T_pi8)
        mod = LLVM.parent(llvm_f)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            # Allocate string
            str_ptr = globalstring_ptr!(builder, _str)
            gv = operands(str_ptr)[1]
            #set_used!(mod, gv)

            # Allocate stack buffer
            malloc_ft = LLVM.FunctionType(T_pi8, [T_i64])
            malloc_f = LLVM.Function(mod, "malloc", malloc_ft)
            buf_ptr = call!(builder, malloc_f, [ConstantInt(T_i64, length(_str)+1)])

            # Copy string into stack buffer
            memcpy_ft = LLVM.FunctionType(LLVM.VoidType(ctx), [T_pi8, T_pi8, T_i64, T_i1])
            memcpy_f = LLVM.Function(mod, "llvm.memcpy.p0i8.p0i8.i64", memcpy_ft)
            call!(builder, memcpy_f, [buf_ptr, str_ptr, ConstantInt(T_i64, length(_str)+1), ConstantInt(T_i1, 0)])

            # Return stack buffer
            ret!(builder, buf_ptr)
        end
        call_function(llvm_f, BufPtr)
    end
end

"""
    UnsizedBuffer <: AbstractUnsizedBuffer

Represents a buffer/string with an unknown size.
"""
struct UnsizedBuffer <: AbstractUnsizedBuffer
    ptr::BufPtr
end
UnsizedBuffer(buf::AbstractBuffer) = UnsizedBuffer(pointer(buf))
UnsizedBuffer(ptr::Ptr) = UnsizedBuffer(reinterpret(BufPtr, ptr))
Base.pointer(str::UnsizedBuffer) = str.ptr
Base.length(str::UnsizedBuffer) = missing


text:
using Random
using Distributions
using JLD2, FileIO
using DelimitedFiles

using RCall
@rimport base
@rimport FlowSOM
@rimport flowCore

# Path to scratch
RESULTS_DIR = "/scratchdata/alui2/cytof/results/"

# CB data path
CB_PATH = "$(RESULTS_DIR)/cb/best/reduced_data/reduced_cb.jld2"

"""
read a single object from a jld2 file
"""
function loadSingleObj(objPath)
  data = load(objPath)
  return data[collect(keys(data))[1]]
end

"""
replace missing values (NaN) in yi with x
"""
function replaceMissing(yi, x)
  out = deepcopy(yi)
  out[isnan.(out)] .= x
  return out
end

function flowsanitize(Y)
  @rput Y
  return R"""
  colnames(Y) = 1:NCOL(Y)
  flowCore::flowFrame(Y)
  """
end

### Main ###
function main()
  # Read data
  y = loadSingleObj(CB_PATH)
  I = length(y)
  N = [size(yi, 1) for yi in y]
  Y = replaceMissing(vcat(y...), -6)
  J = size(Y, 2)
  ff_Y = flowsanitize(Y)

  println("Running FlowSOM ...")
  @time fsom = FlowSOM.FlowSOM(ff_Y,
                               colsToUse=1:J,  # columns to use
                               maxMeta=20,  # Meta clustering option
                               seed=42)  # Seed for reproducible results:

  idx_upper = cumsum(N)
  idx_lower = [1; idx_upper[1:end-1] .+ 1]
  idx = [idx_lower idx_upper]

  fsmeta = fsom[:metaclustering]
  fsclus = fsmeta[Int.(convert(Matrix, fsom[:FlowSOM][:map][:mapping])[:, 1])]

  outpath = "results/cb/flowsom"
  mkpath(outpath)
  for i in 1:I
    # print output
    clus_i = fsclus[idx[i, 1]:idx[i, 2]]

    writedlm("$(outpath)/clustering_$(i).txt", clus_i)
  end
end

main()


text:
function cheb2jac(c_cheb::AbstractVector{T},α::T,β::T,plan::ChebyshevJacobiPlan{BACKWARD,T}) where T<:AbstractFloat
    M,N,nM₀,αN,K = getconstants(plan)
    i₁,i₂,j₁,j₂ = getindices(plan)
    p₁,p₂,rp,c₁,c₂,um,vm,cfs,θ,tempcos,tempsin,tempcosβsinα,tempmindices,cnαβ,cnmαβ,w,anαβ,c_cheb2,pr = getplan(plan)

    # Initialize c_cheb2
    @inbounds for i=1:N+1 c_cheb2[i] = c_cheb[i] end
    @inbounds for i=N+2:2N+1 c_cheb2[i] = 0 end

    # Initialize values at Clenshaw-Curtis points
    v_jac = zero(c_cheb2)

    # Perform transposed DCT-I
    applyTN!(c_cheb2,p₁)

    # Scale values by Clenshaw-Curtis weights
    @inbounds for i=1:2N+1 c_cheb2[i] *= w[i] end

    k=0
    for κ = 1:K
        k = κ
        copyto!(cnmαβ,cnαβ)

        j₁[k] < nM₀ && (k-=1; break)

        @inbounds for i=i₁[k]:i₁[k+1]-1
            forward_recurrence!(pr,j₂[k],θ[i],tempcos[i],tempsin[i],rp)
            cc2i = c_cheb2[i]
            for j=1:j₂[k] v_jac[j] += cc2i*pr[j] end
        end

        @inbounds for i=i₂[k+1]+1:i₂[k]
            forward_recurrence!(pr,j₂[k],θ[i],tempcos[i],tempsin[i],rp)
            cc2i = c_cheb2[i]
            for j=1:j₂[k] v_jac[j] += cc2i*pr[j] end
        end

        for m=0:M-1
            # Compute u_m(θ) and v_m(θ)
            compute_umvm!(um,vm,cfs,α,β,tempcos,tempsin,tempcosβsinα,m,θ,i₁[k+1]:i₂[k+1])

            # Multiply point-wise by u_m(θ) and v_m(θ) for valid indices
            init_c₁c₂!(c₁,c₂.parent,um,vm,c_cheb2,i₁[k+1],i₂[k+1])

            # Apply planned DCT-I and DST-I in-place
            applyTN!(c₁,p₁);applyUN!(c₂,p₂)

            # Compute diagonal 2N-scaling multiplied by local coefficients and zero out excess
            @inbounds for j=j₁[k]:j₂[k] v_jac[j] += cnmαβ[j]*(c₁[j]+c₂.parent[j]) end

            # Update C_{n,m}^{α,β} by recurrence in m
            @inbounds for j=1:2N+1 cnmαβ[j] /= 2(2j+α+β+m) end
        end
    end

    # Finish off recurrence
    @inbounds for i=i₁[k+1]:i₂[k+1]
        forward_recurrence!(pr,j₂[k+1],θ[i],tempcos[i],tempsin[i],rp)
        cc2i = c_cheb2[i]
        for j=1:j₂[k+1] v_jac[j] += cc2i*pr[j] end
    end

    # Scale by orthonormality constants
    ret = zero(c_cheb)
    @inbounds for i=1:N+1 ret[i] = v_jac[i]/anαβ[i] end

    ret
end