text:
using Colors
# http://www.glvisualize.com/examples/meshes/
baselen = 0.1f0
dirlen = 2f0
# create an array of differently colored boxes in the direction of the 3 axes
rectangles = [
(HyperRectangle{3,Float32}(Vec3f0(-baselen/2), Vec3f0(baselen)), RGBA(0f0,0f0,0f0,1f0)),
(HyperRectangle{3,Float32}(Vec3f0(baselen/2, -baselen/2, -baselen/2), Vec3f0(dirlen, baselen, baselen)), RGBA(1f0,0f0,0f0,1f0)),
(HyperRectangle{3,Float32}(Vec3f0(-baselen/2, baselen/2, -baselen/2), Vec3f0(baselen, dirlen, baselen)), RGBA(0f0,1f0,0f0,1f0)),
(HyperRectangle{3,Float32}(Vec3f0(-baselen/2, -baselen/2, baselen/2), Vec3f0(baselen, baselen, dirlen)), RGBA(0f0,0f0,1f0,1f0))
]
# convert to an array of normal meshes
# note, that the constructor is a bit weird. GLNormalMesh takes a tuple of
# a geometry and a color. This means, the geometry will be converted to a GLNormalMesh
# and the color will be added afterwards, so the resulting type is a GLNormalColorMesh
meshes = map(GLNormalMesh, rectangles)
# merge them into one big mesh
# the resulting type is a GLNormalAttributeMesh, since we merged meshes with different
# attributes (colors). An array of the colors will be created and each vertex in the
# mesh will be asigned to one of the colors found there.
colored_mesh = merge(meshes)
view(visualize(colored_mesh), window)


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

ENV["GKSwstype"] = "nul"
using Plots
using Measures

using Printf
using Statistics

using Oceananigans
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary

#####
##### The Bickley jet
#####

Ψ(y) = - tanh(y)
U(y) = sech(y)^2

# A sinusoidal tracer
C(y, L) = sin(2π * y / L)

# Slightly off-center vortical perturbations
ψ̃(x, y, ℓ, k) = exp(-(y + ℓ/10)^2 / 2ℓ^2) * cos(k * x) * cos(k * y)

# Vortical velocity fields (ũ, ṽ) = (-∂_y, +∂_x) ψ̃
ũ(x, y, ℓ, k) = + ψ̃(x, y, ℓ, k) * (k * tan(k * y) + y / ℓ^2)
ṽ(x, y, ℓ, k) = - ψ̃(x, y, ℓ, k) * k * tan(k * x)

"""
    run_bickley_jet(output_time_interval = 2, stop_time = 200, arch = CPU(), Nh = 64, ν = 0,
                    momentum_advection = VectorInvariant())

Run the Bickley jet validation experiment until `stop_time` using `momentum_advection`
scheme or formulation, with horizontal resolution `Nh`, viscosity `ν`, on `arch`itecture.
"""
function run_bickley_jet(; output_time_interval = 2, stop_time = 200, arch = CPU(), Nh = 64, ν = 0, advection = WENO5())

    # Regular model
    grid = RectilinearGrid(size=(Nh, Nh), halo=(3, 3),
                                  x = (-2π, 2π), y=(-2π, 2π),
                                  topology = (Periodic, Bounded, Flat))

    regular_model = NonhydrostaticModel(architecture = arch,
                                        advection = advection,
                                        timestepper = :RungeKutta3,
                                        grid = grid,
                                        tracers = :c,
                                        closure = IsotropicDiffusivity(ν=ν, κ=ν),
                                        coriolis = nothing,
                                        buoyancy = nothing)

    # Non-regular model
    solid(x, y, z) = y > 2π

    expanded_grid = RectilinearGrid(size=(Nh, Int(5Nh/4)), halo=(3, 3),
                                           x = (-2π, 2π), y=(-2π, 3π),
                                           topology = (Periodic, Bounded, Flat))

    immersed_grid = ImmersedBoundaryGrid(expanded_grid, GridFittedBoundary(solid))

    immersed_model = NonhydrostaticModel(architecture = arch,
                                         advection = advection,
                                         timestepper = :RungeKutta3,
                                         grid = immersed_grid,
                                         tracers = (:c, :mass),
                                         closure = IsotropicDiffusivity(ν=ν, κ=ν),
                                         coriolis = nothing,
                                         buoyancy = nothing)

    # ** Initial conditions **
    #
    # u, v: Large-scale jet + vortical perturbations
    #    c: Sinusoid

    # Parameters
    ϵ = 0.1 # perturbation magnitude
    ℓ = 0.5 # Gaussian width
    k = 0.5 # Sinusoidal wavenumber

    # Total initial conditions
    uᵢ(x, y, z) = U(y) + ϵ * ũ(x, y, ℓ, k)
    vᵢ(x, y, z) = ϵ * ṽ(x, y, ℓ, k)
    cᵢ(x, y, z) = C(y, grid.Ly)

    set!(regular_model, u=uᵢ, v=vᵢ, c=cᵢ)
    set!(immersed_model, u=uᵢ, v=vᵢ, c=cᵢ, mass=1)

    wall_clock = [time_ns()]

    function progress(sim)
        @info(@sprintf("Iter: %d, time: %.1f, Δt: %.3f, wall time: %s, max|u|: %.2f",
                       sim.model.clock.iteration,
                       sim.model.clock.time,
                       sim.Δt.Δt,
                       prettytime(1e-9 * (time_ns() - wall_clock[1])),
                       maximum(abs, sim.model.velocities.u.data.parent)))

        wall_clock[1] = time_ns()

        return nothing
    end

    models = (immersed_model, regular_model)
    @show experiment_name = "bickley_jet_Nh_$(Nh)_$(typeof(regular_model.advection).name.wrapper)"

    for m in models
        wizard = TimeStepWizard(cfl=0.1, Δt=0.1 * grid.Δx, max_change=1.1, max_Δt=10.0)

        simulation = Simulation(m, Δt=wizard, stop_time=stop_time, iteration_interval=10, progress=progress)

        # Output: primitive fields + computations
        u, v, w, c = merge(m.velocities, m.tracers)
        ζ = ComputedField(∂x(v) - ∂y(u))
        outputs = merge(m.velocities, m.tracers, (ζ=ζ,))

        output_name = m.grid isa ImmersedBoundaryGrid ?
                            "immersed_" * experiment_name :
                            "regular_" * experiment_name

        @show output_name

        simulation.output_writers[:fields] =
            JLD2OutputWriter(m, outputs,
                             schedule = TimeInterval(output_time_interval),
                             prefix = output_name,
                             field_slicer = nothing,
                             force = true)

        @info "Running a simulation of an unstable Bickley jet with $(Nh)² degrees of freedom..."

        start_time = time_ns()

        run!(simulation)
    end

    return experiment_name
end

"""
    visualize_bickley_jet(experiment_name)

Visualize the Bickley jet data associated with `experiment_name`.
"""
function visualize_immersed_bickley_jet(experiment_name)

    @info "Making a fun movie about an unstable Bickley jet..."

    regular_filepath = "regular_" * experiment_name * ".jld2"
    immersed_filepath = "immersed_" * experiment_name * ".jld2"

    regular_u_timeseries = FieldTimeSeries(regular_filepath,  "u")
    immersed_u_timeseries = FieldTimeSeries(immersed_filepath, "u")

    regular_v_timeseries = FieldTimeSeries(regular_filepath,  "v")
    immersed_v_timeseries = FieldTimeSeries(immersed_filepath, "v")

    regular_ζ_timeseries = FieldTimeSeries(regular_filepath,  "ζ")
    immersed_ζ_timeseries = FieldTimeSeries(immersed_filepath, "ζ")

    regular_c_timeseries = FieldTimeSeries(regular_filepath,  "c")
    immersed_c_timeseries = FieldTimeSeries(immersed_filepath, "c")

    regular_grid = regular_c_timeseries.grid
    immersed_grid = immersed_c_timeseries.grid

    xu, yu, zu = nodes(regular_u_timeseries)
    xv, yv, zv = nodes(regular_v_timeseries)
    xζ, yζ, zζ = nodes(regular_ζ_timeseries)
    xc, yc, zc = nodes(regular_c_timeseries)

    anim = @animate for (i, t) in enumerate(regular_c_timeseries.times)

        @info "    Plotting frame $i of $(length(regular_c_timeseries.times))..."

        Nx, Ny, Nz = size(regular_grid)

        regular_u = regular_u_timeseries[i]
        regular_v = regular_v_timeseries[i]
        regular_ζ = regular_ζ_timeseries[i]
        regular_c = regular_c_timeseries[i]

        immersed_u = immersed_u_timeseries[i]
        immersed_v = immersed_v_timeseries[i]
        immersed_ζ = immersed_ζ_timeseries[i]
        immersed_c = immersed_c_timeseries[i]

        regular_ui = interior(regular_u)[:, :, 1]
        regular_vi = interior(regular_v)[:, :, 1]
        regular_ζi = interior(regular_ζ)[:, :, 1]
        regular_ci = interior(regular_c)[:, :, 1]

        Nx, Ny = size(regular_ui)

        immersed_ui = interior(immersed_u)[1:Nx, 1:Ny, 1]
        immersed_vi = interior(immersed_v)[1:Nx, 1:Ny+1, 1]
        immersed_ζi = interior(immersed_ζ)[1:Nx, 1:Ny+1, 1]
        immersed_ci = interior(immersed_c)[1:Nx, 1:Ny, 1]

        difference_ui = immersed_ui .- regular_ui
        difference_vi = immersed_vi .- regular_vi
        difference_ζi = immersed_ζi .- regular_ζi
        difference_ci = immersed_ci .- regular_ci

        kwargs = Dict(:aspectratio => 1,
                      :linewidth => 0,
                      :colorbar => :none,
                      :ticks => nothing,
                      :clims => (-1, 1),
                      :xlims => (-regular_grid.Lx/2, regular_grid.Lx/2),
                      :ylims => (-regular_grid.Ly/2, regular_grid.Ly/2))

        regular_u_plot = heatmap(xu, yu, clamp.(regular_ui, -1, 1)'; color = :balance, kwargs...)
        regular_v_plot = heatmap(xv, yv, clamp.(regular_vi, -1, 1)'; color = :balance, kwargs...)
        regular_ζ_plot = heatmap(xζ, yζ, clamp.(regular_ζi, -1, 1)'; color = :balance, kwargs...)
        regular_c_plot = heatmap(xc, yc, clamp.(regular_ci, -1, 1)'; color = :thermal, kwargs...)

        immersed_u_plot = heatmap(xu, yu, clamp.(immersed_ui, -1, 1)'; color = :balance, kwargs...)
        immersed_v_plot = heatmap(xv, yv, clamp.(immersed_vi, -1, 1)'; color = :balance, kwargs...)
        immersed_ζ_plot = heatmap(xζ, yζ, clamp.(immersed_ζi, -1, 1)'; color = :balance, kwargs...)
        immersed_c_plot = heatmap(xc, yc, clamp.(immersed_ci, -1, 1)'; color = :thermal, kwargs...)

        δlim = 0.1
        kwargs[:clims] = (-δlim, δlim)

        difference_u_plot = heatmap(xu, yu, clamp.(difference_ui, -δlim, δlim)'; color = :balance, kwargs...)
        difference_v_plot = heatmap(xv, yv, clamp.(difference_vi, -δlim, δlim)'; color = :balance, kwargs...)
        difference_ζ_plot = heatmap(xζ, yζ, clamp.(difference_ζi, -δlim, δlim)'; color = :balance, kwargs...)
        difference_c_plot = heatmap(xc, yc, clamp.(difference_ci, -δlim, δlim)'; color = :thermal, kwargs...)

        r_u_title = @sprintf("regular u at t = %.1f", t)
        r_v_title = @sprintf("regular v at t = %.1f", t)
        r_ζ_title = @sprintf("regular ζ at t = %.1f", t)
        r_c_title = @sprintf("regular c at t = %.1f", t)

        i_u_title = @sprintf("immersed u at t = %.1f", t)
        i_v_title = @sprintf("immersed v at t = %.1f", t)
        i_ζ_title = @sprintf("immersed ζ at t = %.1f", t)
        i_c_title = @sprintf("immersed c at t = %.1f", t)

        d_u_title = @sprintf("Δu at t = %.1f", t)
        d_v_title = @sprintf("Δv at t = %.1f", t)
        d_ζ_title = @sprintf("Δζ at t = %.1f", t)
        d_c_title = @sprintf("Δc at t = %.1f", t)

        plot(regular_u_plot, regular_v_plot, regular_ζ_plot, regular_c_plot,
             immersed_u_plot, immersed_v_plot, immersed_ζ_plot, immersed_c_plot,
             difference_u_plot, difference_v_plot, difference_ζ_plot, difference_c_plot,
             title = [r_u_title r_v_title r_ζ_title r_c_title i_u_title i_v_title i_ζ_title i_c_title d_u_title d_v_title d_ζ_title d_c_title],
             layout = (3, 4), size = (2000, 1000))
    end

    mp4(anim, "differences_" * experiment_name * ".mp4", fps = 8)
end

"""
    analyze_bickley_jet(experiment_name)

Analyze the Bickley jet data associated with `experiment_name`.
"""

function analyze_immersed_bickley_jet(experiment_name)

    @info "    Analyzing IBM Results for Velocity and Tracer Concentration... "

    regular_filepath = "regular_" * experiment_name * ".jld2"
    immersed_filepath = "immersed_" * experiment_name * ".jld2"

    regular_v_timeseries = FieldTimeSeries(regular_filepath,  "v")
    immersed_v_timeseries = FieldTimeSeries(immersed_filepath, "v")

    immersed_m_timeseries = FieldTimeSeries(immersed_filepath, "mass")
    regular_c_timeseries = FieldTimeSeries(regular_filepath,  "c")

    xv, yv, zv = nodes(regular_v_timeseries)
    xm, ym, zm = nodes(immersed_m_timeseries)

    # Finding the rms surface normal velocity (should be zero)
    function rms_normal(r_v_series, im_v_series)

        time_amt = length(r_v_series.times);
        last_y_fluid = size(r_v_series)[2];

        r_norm_rms = sqrt.(sum(r_v_series[:,last_y_fluid,1,:].^2, dims = 2)./time_amt)
        im_norm_rms = sqrt.(sum(im_v_series[:,last_y_fluid,1,:].^2, dims = 2)./time_amt)

        return r_norm_rms, im_norm_rms
    end

    regular_norm_rms, immersed_norm_rms = rms_normal(regular_v_timeseries, immersed_v_timeseries);

    # largest value for plotting
    max_norm = round(maximum(immersed_norm_rms), sigdigits = 2, RoundUp);

    @info "    Plotting the surface normal velocity..."

    norm_plot = plot(regular_v_timeseries.grid.xᶜᵃᵃ, regular_norm_rms, label = "regular", yformatter = :scientific,
    color = :red, lw = 3, xlabel = "x", ylabel = "Vⁿ", legend = :bottomright)
    plot!(regular_v_timeseries.grid.xᶜᵃᵃ, immersed_norm_rms, label = "immersed", color = :blue, lw = 3)
    plot!(yticks = 0:max_norm/5:max_norm, guidefontsize = 14, titlefont=14,legendfont = 10,
        tickfont = 8)

    # Finding the area integrated concentration in time

    function area_int_concentration(r_m_series, im_m_series, xm, ym)

        last_y = size(r_m_series)[2];

        im_C = (sum(im_m_series[1:last_y,1:last_y,1,:], dims = 1:2)[1,1,:])*(ym[2]-ym[1])*(xm[2]-xm[1]);

        return im_C
    end

    immersed_concent = area_int_concentration(regular_c_timeseries, immersed_m_timeseries, xm, ym);

    #taking the percent leakage between IBM and nonIBM (which does not change from initial)
    percent_leakage_m = (abs.(immersed_concent[1] .- immersed_concent) ./ immersed_concent[1]) * 100;

    @info "    Plotting the tracer concentration leakage..."

    concent_diff = plot(immersed_m_timeseries.times, percent_leakage_m, yformatter = :scientific,
        color = :blue, xlabel = "t", ylabel = "% Concentration Leakage", legend = false, lw = 3,
        guidefontsize = 14, titlefont=14,legendfont = 10, tickfont = 8)

    results_plot = plot(norm_plot, concent_diff,
             title = ["Normal Velocity" "% Change in Area Integrated Concentration"],
             layout = (1, 2), size = (1400, 500), left_margin=10mm, bottom_margin=10mm)

    Plots.savefig(results_plot, "Analysis_" * experiment_name * ".png")

end

advection = WENO5()
experiment_name = run_bickley_jet(advection=advection, Nh=64, stop_time=200)
visualize_immersed_bickley_jet(experiment_name)
analyze_immersed_bickley_jet(experiment_name)

ext: jl


text:
using TypeDBClient: Label, scoped_name, RelationType, Remote
using TypeDBClient.typedb.protocol: _Type
using Test

@testset "Label" begin

    let x = Label("a")
        @test x.name == "a"
        @test x.scope === nothing
        @test scoped_name(x) == "a"
    end

    let y = Label("s", "a")
        @test y.name == "a"
        @test y.scope === "s"
        @test scoped_name(y) == "s:a"
    end

end

@testset "RelationType" begin
    type_proto = _Type(label = "a", root = false)
    relation_type = RelationType(type_proto)
end


text:
using LowRank, Test

begin
  @time @testset "10 pts, 2d" begin include("uniform2d.jl") end
  @time @testset "30 pts, 2d" begin include("N_30_d_2.jl") end
  @time @testset "50 pts, 5d" begin include("N_50_d_5.jl") end
  @time @testset "non-zero error" begin include("non_zero_error.jl") end
end


text:
#script that estimates the dir bin gas network model on emid data and evaluates GAS
# forecasting performances

using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,DynNets
using PyCall; pygui(:qt); using PyPlot
#estimate and save for half the datase (after LTRO) or whole data?

## Load dataj
using JLD2
 halfPeriod = false
 fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
 loadFilePartialName = "Weekly_eMid_Data_from_"
 halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
 @load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks,
                YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

#

#density plot
matY_T = YeMidWeekly_T[:,:,3:end]
 matA_T = matY_T.>0
    N = length(matA_T[1,:,1])
    T = size(matA_T)[3]
    Ttrain = 100#round(Int, T/2) #70 106 #
    threshVar = 0.00#5
    degsIO_T = [sumSq(matA_T,2);sumSq(matA_T,1)]
 allFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(degsIO_T); identPost = false,identIter= true )
 sizeLab = 30
# density plot
# density_T = [sum(matA_T[:,:,t]) for t=1:T]./(N*(N-1))
# sizeLab = 30
#  close()
#  #subplot(2,1,1) ;
#  plot(density_T)
#  ylabel("Network Density",size = sizeLab)
#  #xlabel("Time",size = sizeLab)
#  grid()
#  plt[:axvline](x=Ttrain,color = "r",linestyle = "--")
#  S_T = [sum(matY_T[:,:,t]) for t=1:T]

#Estimate gas model on train sample
allFitConstTrain,~,~ =  StaticNets.estimate( StaticNets.ErgmDirBin1(meanSq(degsIO_T[:,1:Ttrain],2)) )
 modGasDirBin1_eMidTrain = DynNets.SdErgmDirBin1(degsIO_T[:,1:Ttrain],"FISHER-DIAG")
 estTargDirBin1_eMidTrain,~ = DynNets.estimateTarg(modGasDirBin1_eMidTrain;SSest = allFitSS )
 gasParEstOnTrain = estTargDirBin1_eMidTrain
 modAllObs =  DynNets.SdErgmDirBin1(degsIO_T,"FISHER-DIAG")
 GasforeFit,~ = DynNets.score_driven_filter_or_dgp(modAllObs,[gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])
 gasforeFit = Float64.(GasforeFit)

load_fold = "./data/multiStepForecast/"
 @load(load_fold* "aucForVariousNsampleAndNsteps_final.jld",
          aucGas, valsNsteps,valsNsample)

N_Nsteps,N_Nsample,Nforetype = size(aucGas)

#@load(save_fold* "aucForVariousNsampleAndNsteps_testMat.jld",
#       aucGas, valsNsteps,valsNsample)

sum(mean(expMat_T_means[:,:,2,1,:],3).>0.5)
close()
  tmp_means = [sum(expMat_T_means[:,:,2,1,i].>0.5) for i=1:N_Nsample]
  tmp_Y = [sum(expMat_T_Y[:,:,2,1,i].>0.5) for i=1:N_Nsample]
 subplot(2,1,1)
 indSteps = 1
 title("N steps =  $(valsNsteps[indSteps])")
 plot(valsNsample,tmp_means,".-")
 #plot(valsNsample,aucGas[indSteps,:,:],".-")
 subplot(2,1,2)
 indSam = 3
 title("N sample =  $(valsNsample[indSam])")
 plot(valsNsample,tmp_Y,".-")
 grid()
 #plot(tmp_Y)
 #plot(valsNsteps, aucGas[:,indSam,1:2],".-")

mean(tmp)
expMat_T[:,:,:,1] == expMat_T[:,:,:,4]

close()#
  remMat =squeeze(prod(.!matA_T,3),3)#
  sizeLeg = 15
  maxNsteps = 50
  legTex = []
  aucStat = zeros(maxNsteps)
    for n=1:maxNsteps
          shiftVal = n
          tmpTm1 = StaticNets.nowCastEvalFitNet(   allFitSS[:,Ttrain+1:end],matA_T[:,:,Ttrain+1:end];mat2rem = remMat,shift = shiftVal)
          aucStat[n] = tmpTm1[3]
          legTex = [legTex; "t-$(shiftVal) Par  AUC = $(round(tmpTm1[3],3))"]
    end

  legend(legTex,fontsize = sizeLeg)
  title("T train = $(Ttrain) of $(T)  ",size = sizeLab)

  xlabel("False Positive Rate",size = sizeLab)
  ylabel("False Negative Rate",size = sizeLab)

close()

    legTex = "ERGM"
    Tfore = 16
    plot(1:Tfore,aucStat[1:Tfore,:],"r")
    for indSam=2
    plot(valsNsteps[1:Tfore],aucGas[1:Tfore,indSam,1])
     legTex = [legTex; "SDERGM " ]
    end
    legend(legTex,fontsize = sizeLeg)
    title("Multi Steps Ahead Forecast   ",size = sizeLab)
    xlabel("Number Of Steps Ahead",size = sizeLab)
    ylabel("Area Under The Curve",size = sizeLab)
    grid()

#
# function forecastEvalGasNetDirBin1(obsNet_T::BitArray{3}, gasParEstOnTrain::Array{Array{Float64,1},1},Ttrain::Int,Nsteps::Int;Nsample::Int=100, thVarConst = 0.0005,mat2rem=falses(obsNet_T[:,:,1]) )
#     # this function evaluates the forecasting performances of the GAS model over the train sample.
#     # to do so uses gas parameters previously estimated and observations at time t of the test sample
#     # to forecast tv par at t+1, then using these estimates and t+1 observations goes to t+2 and so on
#
#     @show N2 = length(gasParEstOnTrain[1]);N = round(Int,N2/2)
#     T = length(obsNet_T[1,1,:])
#     @show Ttest = T-Ttrain
#     testObsNet_T = obsNet_T[:,:,Ttrain+1:end]
#     degsIO_T = [sumSq(obsNet_T,2); sumSq(obsNet_T,1)]
#
#
#
#     #disregard predictions for degrees that are constant in the training sample
#     # gas shoud manage but AR fitness no and I dont want to advantage gas
#     # to do so i need to count the number of links that are non constant (nnc)
#     Nc,isConIn,isConOut = StaticNets.defineConstDegs(degsIO_T[:,1:Ttrain];thVarConst = thVarConst )
#     noDiagIndnnc = putZeroDiag(((.!isConIn).*(.!isConOut')).* (.!mat2rem)    )
#
#     @show Nlinksnnc =sum(noDiagIndnnc)
#     # forecast fitnesses using Gas parameters and observations
#     foreFit,~ = score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])
#
#     TRoc = Ttest-1
#     #storage variables
#     foreVals = 100 * ones(Nsteps*Nlinksnnc*(TRoc))
#     realVals = trues( Nsteps*Nlinksnnc*(TRoc))
#
#     lastInd=1
#
#     for t=(1+Nsteps):TRoc
#         adjMat = testObsNet_T[:,:,t]   #[testObsNet_T[50:end,:,t];testObsNet_T[1:49,:,t]] #a shuffle to test that the code does not not work
#         adjMat_tm1 = testObsNet_T[:,:,t-1]
#         indZeroR = sum(adjMat_tm1,2).==0
#         indZeroC = sum(adjMat_tm1,1).==0
#         indZeroMat = indZeroR.*indZeroC
#     #    println(sum(indZeroMat)/(length(indZeroC)^2))
#
#         expMat = StaticNets.expMatrix(StaticNets.fooErgmDirBin1,foreFit[:,Ttrain+1:end][:,t])
#
#          tmpAllMat,~ = multiSteps_predict_score_driven_par(modGasDirBin1_eMidTrain,N,foreFit[:,Ttrain+1:end][:,t-Nsteps],
#             gasParEstOnTrain[1],gasParEstOnTrain[2],gasParEstOnTrain[3],Nsample,Nsteps)
#
#          expMat=tmpAllMat[:,:,Nsteps+1]
#
#
#         #expMat[indZeroMat] = 0
#
#         foreVals[lastInd:lastInd+Nlinksnnc-1] = expMat[noDiagIndnnc]
#         realVals[lastInd:lastInd+Nlinksnnc-1] = adjMat[noDiagIndnnc]
#         lastInd += Nlinksnnc
#     end
#     return realVals,foreVals,foreFit
#  end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
using Octavian, VectorizationBase, ProgressMeter
using Octavian: StaticFloat64
function matmul_pack_ab!(C, A, B, ::Val{W₁}, ::Val{W₂}, ::Val{R₁}, ::Val{R₂}) where {W₁, W₂, R₁, R₂}
    M, N = size(C); K = size(B,1)
    zc, za, zb = Octavian.zstridedpointer.((C,A,B))
    nspawn = min(Threads.nthreads(), VectorizationBase.num_cores())
    @elapsed(
        Octavian.matmul_pack_A_and_B!(
            zc, za, zb, StaticInt{1}(), StaticInt{0}(), M, K, N, nspawn,
            StaticFloat64{W₁}(), StaticFloat64{W₂}(), StaticFloat64{R₁}(), StaticFloat64{R₂}()
        )
    )
end

function bench_size(Cs, As, Bs, ::Val{W₁}, ::Val{W₂}, ::Val{R₁}, ::Val{R₂}) where {W₁, W₂, R₁, R₂}
    if length(first(Cs)) < length(last(Cs))
        matmul_pack_ab!(first(Cs), first(As), first(Bs), Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    else
        matmul_pack_ab!(last(Cs), last(As), last(Bs), Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    end
    gflop = 0.0
    for (C,A,B) ∈ zip(Cs,As,Bs)
        M, K, N = Octavian.matmul_sizes(C, A, B)
        # sleep(0.5)
        t = matmul_pack_ab!(C, A, B, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
        gf = 2e-9M*K*N / t
        gflop += gf
    end
    gflop / length(As)
end
matrix_sizes(s::Int) = (s,s,s)
matrix_sizes(MKN::NTuple{3,Int}) = MKN
size_range(l, u, len) = round.(Int, exp.(range(log(l), stop = log(u), length = len)))
function matrix_range(l, u, len, ::Type{T} = Float64) where {T}
    matrix_range(size_range(l, u, len), T)
end
function matrix_range(S, ::Type{T} = Float64) where {T}
    Alen = 0; Blen = 0; Clen = 0;
    for s ∈ S
        M, K, N = matrix_sizes(s)
        Alen = max(Alen, M*K)
        Blen = max(Blen, K*N)
        Clen = max(Clen, M*N)
    end
    Abuf = rand(T, Alen)
    Bbuf = rand(T, Blen)
    Cbuf = rand(T, Clen)
    As = Vector{Base.ReshapedArray{T, 2, SubArray{T, 1, Vector{T}, Tuple{Base.OneTo{Int}}, true}, Tuple{}}}(undef, length(S))
    Bs = similar(As); Cs = similar(As);
    for (i,s) ∈ enumerate(S)
        M, K, N = matrix_sizes(s)
        As[i] = reshape(view(Abuf, Base.OneTo(M * K)), (M, K))
        Bs[i] = reshape(view(Bbuf, Base.OneTo(K * N)), (K, N))
        Cs[i] = reshape(view(Cbuf, Base.OneTo(M * N)), (M, N))
    end
    Cs, As, Bs
end

T = Float64
min_size = round(Int, sqrt(0.65 * Octavian.VectorizationBase.cache_size(Val(3)) / sizeof(T)))
max_size = round(Int, sqrt( 32  * Octavian.VectorizationBase.cache_size(Val(3)) / sizeof(T)))

SR = size_range(max_size, min_size, 100);
const CsConst, AsConst, BsConst = matrix_range(SR, T);

function matmul_objective(params)
    print("Params: ", params, "; ")
    W₁, W₂, R₁, R₂ = params
    # print("(W₁ = $(round(W₁, sigdigits=4)); W₂ = $(round(W₂, sigdigits=4)); R₁ = $(round(R₁, sigdigits=4)); R₂ = $(round(R₂, sigdigits=4))); ")
    gflop = bench_size(CsConst, AsConst, BsConst, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    println(gflop)
    - gflop
end

using Optim
hours = 60.0*60.0; days = 24hours;
init = Float64[Octavian.W₁Default(), Octavian.W₂Default(), Octavian.R₁Default(), Octavian.R₂Default()]

opt = Optim.optimize(
    matmul_objective, init, ParticleSwarm(lower = [0.001, 0.01, 0.3, 0.4], upper = [0.2, 2.0, 0.9, 0.99]),
    Optim.Options(iterations = 10^6, time_limit = 8hours)
);


text:
#!/usr/bin/env julia

using Cumulants
using SymmetricTensors
using CumulantsUpdates
using JLD2
using FileIO
using ArgParse
using Distributed
import CumulantsUpdates: cnorms

function cumnorms(c::Vector{SymmetricTensor{Float64}})
  map(norm, c[1:2])..., cnorms(c)...
end

function comptime(dm::DataMoments{Float64}, Xup::Matrix{Float64})
  t = time_ns()
  c = cumulantsupdate!(dm, Xup)
  cumnorms(c)
  Float64(time_ns()-t)/1.0e9
end

function precomp(m::Int)
  X = randn(15, 10)
  dm = DataMoments(X[1:10,:], m, 4)
  cumulantsupdate!(dm, X[10:15,:])
end

function savect(u::Vector{Int}, nvec::Vector{Int}, m::Int, p::Int, b::Int)
  comptimes = zeros(length(nvec), length(u))
  precomp(m)
  i = 1
  for n in nvec
    X = randn(maximum(u)+10, n)
    println("n = ", n)
    dm = DataMoments(X, m, b)
    for k in 1:length(u)
      Xup = randn(u[k], n)
      a = comptime(dm, Xup)
      comptimes[i, k] = u[k]/a
      println("u = ", u[k], "f =", div(u[k],a), "Hz")
    end
    i += 1
  end
  filename = replace("res/$(m)_$(u)_$(nvec)_$(p)_$(b)_maxf.jld2", "["=>"")
  filename = replace(filename, "]"=>"")
  filename = replace(filename, " "=>"")
  compt = Dict{String, Any}("cumulants"=> comptimes)
  push!(compt, "t" => u)
  push!(compt, "n" => nvec)
  push!(compt, "m" => m)
  push!(compt, "x" => "n")
  save(filename, compt)
end

function main(args)
  s = ArgParseSettings("description")
  @add_arg_table s begin
      "--order", "-d"
        help = "d, the order of cumulant, ndims of cumulant's tensor"
        default = 4
        arg_type = Int
      "--nvar", "-n"
        help = "n number of marginal variables"
        nargs = '*'
        default = [20, 24, 28]
        help = "n, numbers of marginal variables"
        arg_type = Int
      "--tup", "-u"
        help = "u, numbers of data updates"
        nargs = '*'
        default = [10000, 20000]
        arg_type = Int
      "--nprocs", "-p"
        help = "number of processes"
        default = 3
        arg_type = Int
      "--blocksize", "-b"
         help = "the size of blocks of the block structure"
         default = 4
         arg_type = Int
    end
  parsed_args = parse_args(s)
  m = parsed_args["order"]
  n = parsed_args["nvar"]
  u = parsed_args["tup"]
  p = parsed_args["nprocs"]
  b = parsed_args["blocksize"]
  if p > 1
    addprocs(p)
    eval(Expr(:toplevel, :(@everywhere using CumulantsUpdates)))
  end
  println("number of workers = ", nworkers())
  savect(u, n, m, p, b)
end

main(ARGS)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 232c24de-0e30-11eb-1544-53cb0b45ec9b
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		"Colors",
		"Plots",
		"PlutoUI",
		])
	using Colors
	using Plots
	using PlutoUI
end

# ╔═╡ beed7b48-0e2b-11eb-2805-edf0be651fc3
md"""
### Lecture 22: Solving Partial Differential Equations (PDEs) Numerically
**Part I: One-dimensional advection and diffusion**

#### 1) Advection

Consider the following differential equation for temperature $T(x,t)$,

$\frac{\partial T}{\partial t} = - U \frac{\partial T}{\partial x}$

which describes the rate of change of temperature due to a uniform flow of speed $U$ in the $\mathbf{x}$ direction.

**Notation**:
Since this differential equation includes derivatives in both time $t$ **and** space $x$, we replace the *ordinary* deriative symbol $d$ with a *partial* derivative symbol $\partial$, and call the equation a **partial differential equation (PDE)**.

#### 1.1) Discretizing in time: forward finite difference
**A.k.a. *Euler's Method*$\,$**

As in our solution of the "zero-dimenisonal" energy balance model ODE, we begin by discretizing time using a *forward finite difference* scheme:

$\frac{T_{n+1}-T_{n}}{\Delta t} = - U \frac{\partial T}{\partial x},$

where $\Delta t \equiv t_{n+1} - t_{n}$ is a uniform discrete **timestep**.
"""

# ╔═╡ e74322be-12d6-11eb-3f27-f17821669699
md"
#### 1.2) Discretizing in space: centered finite difference

We now need to discretize the spatial derivative. Unlike for differencing in time, we no longer have a prefered direction for the discretization and so we use a *centered finite difference*:

$\frac{T_{n+1,\, i}-T_{n,\, i}}{\Delta t} = - U \frac{T_{n,\, i+1} - T_{n,\, i-1}}{2 \Delta x},$

where $\Delta x \equiv x_{i+1} - x_{i} = x_{i} - x_{i-1}$ is a uniform discrete **grid spacing**.
"

# ╔═╡ c9a1719a-12d6-11eb-1c8c-13e1f75bc852
md"**Note:** You should think of the *centered* finite difference average of the *right* (forward) and *left* (backward) finite differences:

$\frac{1}{2}\left(\frac{T_{n,\, i+1} - T_{n,\, i}}{\Delta x} + \frac{T_{n,\, i} - T_{n,\, i-1}}{\Delta x} \right) = \frac{T_{n,\, i+1} - T_{n,\, i-1}}{2 \Delta x}$

"

# ╔═╡ c3f99436-0f49-11eb-33b6-e78c150d4ff5
md"""
#### 1.3) Algorithmic procedure: time-stepping

Just as before, we re-order the equation to solve for the temperature of the $i$th grid cell at the next timestep $t_{n+1}$,

$T_{n+1,\, i} = T_{n,\, i} + \Delta t \left(-U \frac{T_{n,\, i+1} - T_{n,\, i-1}}{2 \Delta x} \right),$

which allows us to recursively advance in time, for each $i$, given some initial condition $T_{0,\, i}$ for all $i$.

#### 1.4) Numerical grid and boundary conditions
Consider solving this on a grid $i \in [1,\, N_{i}]$. How do we handle the edge case $i=1$, where $T_{n+1,\, 1}$ depends on $T_{n+1,\, 0}$, which is undefined?

For both this point $i=1$ and the other extreme $i=N_{i}$, we need a **boundary condition**. There are a number of ways of doing this, and we will explore some alternative kinds of boundary conditions later, but the simplest is a *periodic boundary condition*, which wraps the grid around by setting $T_{n,\, 0} = T_{n,\, N_{i}}$ and $T_{n,\, 1} = T_{n,\, N_{i}+1}$
"""

# ╔═╡ 0967dc46-0f4d-11eb-13b2-2b2042595b6e
md"#### 2) Numerical implementation"

# ╔═╡ 19273c8a-0f4d-11eb-3caf-678f073c6b1a
md"##### 2.1) Setting up the model parameters

To keep things simple, let us consider our temperature equation as a model of temperature variations in a one-dimensional ocean current of length $L=1$ m and with a speed $U = 1$ m/s.

The choice of the **discretization resolution** $N_{i}$ is up to the modeller, but here we make it $N_{i} = 10$ because we will be able to easily pick out each of the individual grid cells in plots below.
"

# ╔═╡ e700c752-0e2d-11eb-2c61-3959bdba269a
begin
	nx = 10
	Lx = 1.
	Δx = Lx/nx
	Δt = 0.001
	U = 1.

	x = Δx/2.:Δx:Lx
end;

# ╔═╡ 5b1b6802-0e2e-11eb-1fa6-1ddd5478cc54
begin
	# Initial conditions
	T = sin.(2π*x);
	t = [0.]
end;

# ╔═╡ 4463cf8c-0ef8-11eb-2c6e-45d982d5687a
function advect(T)
	return U*(circshift(T, (1)) .- circshift(T, (-1)))/(2Δx)
end

# ╔═╡ 67cb0ed2-0efa-11eb-2645-495dba94ea64
function timestep!(t, T)
	T .+= Δt*(advect(T))
	t .+= Δt
end

# ╔═╡ ba1ea938-0f44-11eb-0fac-cd1afe507d57
timestepButton = @bind go Button("Timestep")

# ╔═╡ 18a3f9dc-0e6d-11eb-0b31-0b5296c2e83b
begin
	⏩ = nothing
	go
	nT = 50
	for i = 1:nT
		timestep!(t, T)
	end
end;

# ╔═╡ f2c7638a-0e34-11eb-2210-9f9b0c3519fe
function temperature_heatmap(T)
	p = plot(xticks=x, yticks=nothing, size=(700,90))
	plot!(p, x, [0.], reshape(T, (size(T)...,1))', st=:heatmap, clims=(-1., 1.))
end;

# ╔═╡ b19bfb58-0f4e-11eb-218f-8d930b7afff6
md"#### 2) Diffusion

$\frac{\partial T}{\partial t} = \kappa \frac{\partial^{2} T}{\partial x^{2}}$

Again, we discretize the equation by considering the derivatives one at a time:

$\frac{T_{n+1,\, i} - T_{n,\, i}}{\Delta t} = \kappa \frac{\partial^{2} T}{\partial x^{2}}$

$\frac{T_{n+1,\, i} - T_{n,\, i}}{\Delta t} = \kappa \left( \frac{\frac{\partial T}{\partial x}|_{n,\, i+0.5} - \frac{\partial T}{\partial x}|_{n,\, i-0.5}}{\Delta x} \right)$

$\frac{T_{n+1,\, i} - T_{n,\, i}}{\Delta t} = \kappa \left( \frac{\frac{T_{n,\, i+1} - T_{n,\, i}}{\Delta x} - \frac{T_{n,\, i} - T_{n,\, i-1}}{\Delta x}}{\Delta x}\right)$

$\frac{T_{n+1,\, i} - T_{n,\, i}}{\Delta t} = \kappa \left( \frac{T_{n,\, i+1} - 2 T_{n,\, i} + T_{n,\, i-1}}{(\Delta x)^{2}}\right)$
"

# ╔═╡ 37bed6e8-12da-11eb-1c36-c7f92fe732e4
κ = 0.05

# ╔═╡ 3e1d44b6-12da-11eb-16ac-2524d0e0d900
function diffuse(T)
	return κ*(circshift(T, (1)) .- 2*T .+ circshift(T, (-1)))/(Δx^2)
end

# ╔═╡ 43d1272a-2113-11eb-1d03-c948d11a8423
# function timestep!(t, T)
# 	T .+= Δt*(diffuse(T))
# 	t .+= Δt
# end

# ╔═╡ 74a36fd4-2113-11eb-2d3e-3be5e4318dea
timestepButton

# ╔═╡ d2c101c8-0f4e-11eb-0c4f-0972be924ce1
md"#### 3) Advection-Diffusion"

# ╔═╡ 7c490fd2-2113-11eb-3239-f72cc77c0b53
# function timestep!(t, T)
# 	T .+= Δt*(advect(T) .+ diffuse(T))
# 	t .+= Δt
# end

# ╔═╡ 894b97f4-2113-11eb-3b86-2986a6cd4fd1
timestepButton

# ╔═╡ b00f7bec-12da-11eb-3e20-7769cc2bd0ec
md"
#### 4) Finite differences as kernels acting on an array

##### 4.1) Analogies between image processing and PDEs

##### 4.1.1) Blurring (Gaussian kernel) vs. Laplacian operator

##### 4.1.2) Edge detection (Sobel kernel) vs. Gradient operator
"

# ╔═╡ 70ae3138-0f44-11eb-1a8d-d3b5a39c1b42
md"""#### Pluto Book-keeping"""

# ╔═╡ 845e102a-0f44-11eb-3935-2d9a0efecedc
as_svg(x) = PlutoUI.Show(MIME"image/svg+xml"(), repr(MIME"image/svg+xml"(), x))

# ╔═╡ 00cc530a-0e35-11eb-133b-ef8e30ea7b12
begin
	⏩
	p1 = plot(x, T, label="Temperature", ylim=[-1.1, 1.1], xlim=[0., 1.], marker=:c)
	annotate!(p1, [(0.05, 0.9, string("t = ", round(t[1], digits=2)))])
	p2 = temperature_heatmap(T)
	p = plot(p1, p2, layout=grid(2, 1, heights=[0.7 , 0.3]), size=(680,250))
	as_svg(p)
end

# ╔═╡ 4c5b62fa-2113-11eb-14ff-2bd776efe96c
as_svg(p)

# ╔═╡ 869c4be8-2113-11eb-0703-f905fe93142f
as_svg(p)

# ╔═╡ Cell order:
# ╟─beed7b48-0e2b-11eb-2805-edf0be651fc3
# ╟─e74322be-12d6-11eb-3f27-f17821669699
# ╟─c9a1719a-12d6-11eb-1c8c-13e1f75bc852
# ╟─c3f99436-0f49-11eb-33b6-e78c150d4ff5
# ╟─0967dc46-0f4d-11eb-13b2-2b2042595b6e
# ╟─19273c8a-0f4d-11eb-3caf-678f073c6b1a
# ╠═e700c752-0e2d-11eb-2c61-3959bdba269a
# ╠═5b1b6802-0e2e-11eb-1fa6-1ddd5478cc54
# ╠═4463cf8c-0ef8-11eb-2c6e-45d982d5687a
# ╠═67cb0ed2-0efa-11eb-2645-495dba94ea64
# ╠═ba1ea938-0f44-11eb-0fac-cd1afe507d57
# ╠═00cc530a-0e35-11eb-133b-ef8e30ea7b12
# ╠═18a3f9dc-0e6d-11eb-0b31-0b5296c2e83b
# ╠═f2c7638a-0e34-11eb-2210-9f9b0c3519fe
# ╟─b19bfb58-0f4e-11eb-218f-8d930b7afff6
# ╠═37bed6e8-12da-11eb-1c36-c7f92fe732e4
# ╠═3e1d44b6-12da-11eb-16ac-2524d0e0d900
# ╠═43d1272a-2113-11eb-1d03-c948d11a8423
# ╠═74a36fd4-2113-11eb-2d3e-3be5e4318dea
# ╠═4c5b62fa-2113-11eb-14ff-2bd776efe96c
# ╟─d2c101c8-0f4e-11eb-0c4f-0972be924ce1
# ╠═7c490fd2-2113-11eb-3239-f72cc77c0b53
# ╠═894b97f4-2113-11eb-3b86-2986a6cd4fd1
# ╠═869c4be8-2113-11eb-0703-f905fe93142f
# ╟─b00f7bec-12da-11eb-3e20-7769cc2bd0ec
# ╟─70ae3138-0f44-11eb-1a8d-d3b5a39c1b42
# ╠═232c24de-0e30-11eb-1544-53cb0b45ec9b
# ╠═845e102a-0f44-11eb-3935-2d9a0efecedc

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
using Conda, BinDeps, Compat, VersionParsing
using Compat.Test

Conda.update()

env = :test_conda_jl
@test Conda.exists("curl", env)
Conda.add("curl", env)

curlvers = Conda.version("curl",env)
@test curlvers >= v"5.0"
@test Conda.exists("curl==$curlvers", env)

exe = Compat.Sys.iswindows() ? ".exe" : ""

curl_path = joinpath(Conda.bin_dir(env), "curl" * exe)
@test isfile(curl_path)

@test "curl" in Conda.search("cu*", env)

Conda.rm("curl", env)
@test !isfile(curl_path)

pythonpath = joinpath(Conda.PYTHONDIR, "python" * exe)
@test isfile(pythonpath)
pyversion = read(`$pythonpath -c "import sys; print(sys.version)"`, String)
@test pyversion[1:1] == Conda.MINICONDA_VERSION

Conda.add_channel("foo", env)
@test Conda.channels(env) == ["foo", "defaults"]
# Testing that calling the function twice do not fail
Conda.add_channel("foo", env)

Conda.rm_channel("foo", env)
@test Conda.channels(env) == ["defaults"]


text:
using LinearAlgebraicRepresentation, ViewerGL, SparseArrays
Lar = LinearAlgebraicRepresentation; GL = ViewerGL
using Base.union

# 3D Boolean example generation
#-------------------------------------------------------------------------------
n,m,p = 1,1,1
V,(VV,EV,FV,CV) = Lar.cuboidGrid([n,m,p],true)
cube = V,FV,EV

V,FV = Lar.sphere()()
EV = Lar.simplexFacets(FV)
sphere = V,FV,EV

# three cubes in "assembly"
assembly = Lar.Struct([ cube,
    Lar.t(.3,.4,.25), Lar.r(pi/5,0,0), Lar.r(0,0,pi/12), cube,
    Lar.t(-.2,.4,-.2), Lar.r(0,pi/5,0), Lar.r(0,pi/12,0), cube ])

V,FV,EV = Lar.struct2lar(assembly)
GL.VIEW([ GL.GLGrid(V,FV), GL.GLFrame ]);

W, (copEV, copFE, copCF), boolmatrix = Lar.bool3d(assembly)
Matrix(boolmatrix)
#three-chains = [ for k = 1:3]

A = boolmatrix[:,2]
B = boolmatrix[:,3]
C = boolmatrix[:,4]
AorB = A .| B
AandB = A .& B
AxorB = AorB .& (.! AandB)
AorBorC = A .| B .| C
AorBorC = .|(A,B,C)
AandBandC = A .& B .& C
AandBandC = .&(A,B,C)
AminusBminusC = .&(A, .!(B .| C)) # A - B - C

unione = .|(A,B,C)
unione = Matrix(copCF)' * Int.(unione) # coord vector of Faces
intersection = Matrix(copCF)' * Int.(AandBandC) # coord vector of Faces
difference = Matrix(copCF)' * Int.(AminusBminusC) # coord vector of Faces

V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF) # whole assembly
Fs = unione
V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF, Fs) # part of assembly

# EV = Lar.cop2lar(copEV)
# EVor = [ev for (k,ev) in enumerate(EV) if abs(union[k])==1 ]
# EVand = [ev for (k,ev) in enumerate(EV) if abs(intersection[k])==1 ]
# EVxor = [ev for (k,ev) in enumerate(EV) if abs(xor[k])==1 ]

GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.,1.,1.,99,1));
meshes = GL.GLExplode(V,CVs[2:end],1.5,1.5,1.5,99,1);
GL.VIEW( push!( meshes, GL.GLFrame) );

# GL.VIEW(GL.GLExplode(V,[EVor],1.,1.,1.,99,1));

GL.VIEW([ GL.GLGrid(V,EV, GL.COLORS[1],0.5), GL.GLFrame ]);


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

using Markdown
using InteractiveUtils

# ╔═╡ b7433dc1-846b-4173-ad51-f7598cd263a1
using FFTW, Plots, SampledSignals

# ╔═╡ 37f59728-aeae-4b15-ba8c-a90368dfa8f0
begin
    fs = 48000
	freqs = Array(LinRange(0, 2*pi*200, fs))
end

# ╔═╡ cfc3b26d-d2da-4128-9ef2-a9494c4c25e5
function plotfreqs(tone, title)
	freqs = fft(tone) |> fftshift
	domain = fftfreq(length(tone), fs) |> fftshift
	plot(domain, 10*log10.(abs.(freqs)/fs), title=title, xlimit=(0, 2000), ylimit=(-50, 0))
end

# ╔═╡ b70c60fd-b29a-46f3-8daf-e6d7277156ca
function plotfreqs!(tone, title)
	freqs = fft(tone) |> fftshift
	domain = fftfreq(length(tone), fs) |> fftshift
	plot!(domain, 10*log10.(abs.(freqs)/fs), title=title, xlimit=(0, 2000), ylimit=(-50, 0))
end

# ╔═╡ baafd0a2-b449-11eb-0654-172c64b81fd6
md"# The Harmonic Series
Periodic sounds have a base frequency that is the reciprocal of the period,
and also contain components (overtones) of integer multiples of the base frequency.
The amplitudes of the overtones determine the timbre of the sound.
"

# ╔═╡ 7a4f0f06-6b37-4a3c-81e4-8004462da33a
function note(t)
	t = t/pi%2-1
	t = 3*sqrt(3)/2*t*(t*t-1)
end

# ╔═╡ ef073029-a63e-4625-aecf-0cde96af5686
tone = SampleBuf(note.(freqs), fs)

# ╔═╡ 9997ac3d-3cae-4bb5-8290-6ca53f431533
plotfreqs(tone, "Frequencies of a single 200 Hz note")

# ╔═╡ dcd90ce9-dfc0-42ed-990a-55710c32b34e
md"## Intervals
If we play two unrelated notes at the same time we get two harmonic series"

# ╔═╡ 7f363c48-a194-4494-bb59-d5e00db70bd2
begin
   ratio=1+rand()
   plotfreqs(note.(freqs), "First note")
   plotfreqs!(note.(ratio*freqs), "Second note")
end

# ╔═╡ 9e45b6c3-a674-4b1f-b93e-11b4a99a22da
dissonant = SampleBuf((note.(freqs) + note.(ratio*freqs))/2, fs)

# ╔═╡ 07312f0d-59d9-45ee-b60a-c6ecfc2bc0df
# The frequencies of the combined note is the sum of each notes frequencies
plotfreqs(dissonant, "Two notes frequencies")

# ╔═╡ bfd2d6b6-2863-41ee-9b6a-c99ed19f8d7d
md"""## Octave
If we use frequencies in a $1:2$ ratio the harmonic series line up.
Every frequency is an overtone of the lower note.
This single harmonic series has more pronounced even overtones.
Because there is a single harmonic series,
these two notes combined are equivalent to a single note
with a different timbre.
With real world sounds different timbres, the sound source locations,
and innacuracies of analog tuning make this equivalence inexact.
"""

# ╔═╡ 3813605c-e5c8-4bd3-9813-f553d22d47b2
begin
   plotfreqs(note.(freqs), "")
   plotfreqs!(note.(2*freqs), "Octave")
end

# ╔═╡ 6a831a92-a222-48e8-abfb-72bc9f3e0fb9
octave = SampleBuf(note.(freqs)/2 + note.(2*freqs)/2, fs)

# ╔═╡ 0b9c89bf-291e-4b93-8632-f2d5117d25a4
md"If we take a tone with only odd harmonics we can construct an octave with no harmonics in common but that still sounds like a single note because of the shared harmonic series. The notes have to be scaled so that the curve of the harmonics line up to sound cohesive."

# ╔═╡ a4882ae0-32ed-4db4-ba0c-165081128f2d
function odd(t)
	t = t/pi%2-1
	t = 4*t*(abs(t)-1)
end

# ╔═╡ e3faa985-c01a-471c-ae10-0aaaee8a6b43
odd_octave = SampleBuf(odd.(freqs) + odd.(2*freqs)/8 + odd.(4*freqs)/64, fs)

# ╔═╡ 6049ed44-b19f-41a9-b017-860e80d62017
begin
    plotfreqs(odd.(freqs), "")
    plotfreqs!(odd.(2*freqs)/8, "")
	plotfreqs!(odd.(4*freqs)/64, "Odd harmonics combined")
end

# ╔═╡ ccb62252-fcf1-4fc8-b486-1df69ea0c1e3
md"Any $1:n$ ratio gives a similar lining up like the octave.
Unison $1:1$ being a special case that gives the exact original sound but louder.
However, since for natural timbres the higher harmonics sound less loud,
higher intervals harmonics blend less and stand out more from their root note.
If the higher note is made softer they start to blend again.
For example a $1:3$ ratio, or tritive sounds like this."

# ╔═╡ 14af9192-4b8a-49c3-9775-b4dfef343780
tritive = SampleBuf(note.(freqs)/2 + note.(3*freqs)/2, fs)

# ╔═╡ 02a98eb5-4021-4ea0-a0f8-b43c811d87b2
blend_tritive = SampleBuf(note.(freqs)/2 + note.(3*freqs)/8, fs)

# ╔═╡ 31b44424-b3f9-4ea0-8d99-28738f1fc857
md"## Fifth
Though only integer ratios will give harmonics that line up exactly
any rational ratio will have overlap.
A $2:3$ ratio, the simplest non-integer interval, gives a fifth.
The fifth is an octave below a tritave
and so the fifth and the tritave have a similar quality.
Going up a fifth $(3/2)$ is the inverse of going down a fifth $(2/3)$
a characteristic shared by all intervals."

# ╔═╡ 82582e9b-fcc0-49ae-9f3b-558e70d8f0b8
fifth = SampleBuf(note.(freqs)/2 + note.(3/2*freqs)/2, fs)

# ╔═╡ db7c4951-972b-43db-a647-6e61c8ebdf03
begin
   plotfreqs(note.(freqs), "")
   plotfreqs!(note.(3/2*freqs), "Fifth")
end

# ╔═╡ 53773643-1b58-4eee-ac11-bb30b7535fec
md"Every second harmonic of the higher note lines up
with every third harmonic of the lower note.
Unlike playing notes in unison or octaves,
the two notes have distinct pitches,
but unlike a random interval it's consonant.
"

# ╔═╡ 151b3663-524e-465f-9426-3f85689e9e27
md"## The Fundamental
Because of how the overtones line up,
every overtone is part of the harmonic series
of a note that isn't being played.
for the fifth this is one octave below the lower note.
This corresponds to the period of the combined sounds,
which repeats whenever the period of both notes line up.

For an interval $n:m$ that is rational and in reduced terms,
the fundamental is as $1$ in the ratio: $1:n:m$.
So the fundamental is at the inteval 1/n from the first note
and 1/m from the second.
The ratio $1:n$ is a measure of consonance and is equivalent to the porportion of harmonics in the second note that are also in the first.
Alternately the ratio $1:m$ is is another measure and is equivalent to the porportion
of harmonics of the first note present in the second.
Normally the lowest note of a collection
is chosen as the root, the reference point for calculating intervals,
but for chord inversions using a different note is often more convenient.
This measure extends naturally to collections of more than 2 notes."

# ╔═╡ 24e5caba-ce1c-4436-8ca2-647b5438658a
begin
    plotfreqs(note.(freqs), "")
    plotfreqs!(note.(3/2*freqs), "")
    plotfreqs!(note.(freqs/2), "Fifth with fundamental")
end

# ╔═╡ 8adb4c0a-a29e-4cee-bb49-ee285f90ff70
md"""## The Common Overtone
Like the fundamental of a collection of notes
is the highest note that has all of the notes as harmonics,
the dual to the fundamental is the common overtone:
The lowest note that is a harmonic of all of the notes.
The harmonic series of this note contains the frequencies
common to all the harmonic series of the other notes.

For an interval $n:m$ that is rational and in reduced terms,
the ratio with the dual fundamental is $n:m:nm$.
This is the least common multiple of $n,m,$ dual to the fundamental
which is the greatest common denominator.
The ratios $n:nm = 1:m$, and $m:nm = 1:n$ so the
common overtone gives the same measures of consonance as the
fundamental but compared with the opposite notes.
For collections of more than two notes this difference leads to
collections with simple ratios to the fundamental being different
from collections with simple ratios to their common overtone."""

# ╔═╡ 8378910d-3a2c-441e-a3b3-716d819eced9
begin
    plotfreqs(note.(freqs), "")
    plotfreqs!(note.(3/2*freqs), "")
    plotfreqs!(note.(freqs*3), "Fifth with common overtone")
end

# ╔═╡ 732ab799-970e-4ae7-a81d-cda4c79873f4
md" # Fourth
The fourth $3:4$ is an octave equivalent to the fifth.
Taking a fifth and raising the lower note by an octave
(or equivalently lowering the higher note by an octave)
gives a fourth."

# ╔═╡ df7e7e47-fe34-4250-bb0b-670fee0b56de
fourth = SampleBuf(note.(freqs)/2 + note.(freqs*4/3)/2, fs)

# ╔═╡ 518f1a8b-977b-4008-b1b7-667f0475b42c
begin
    plotfreqs(note.(freqs), "")
    plotfreqs!(note.(4/3*freqs), "")
    plotfreqs!(note.(freqs/3), "Fourth with fundamental")
end

# ╔═╡ 849e0964-2c80-4fa6-9320-5fa0c595d79a
begin
    plotfreqs(note.(freqs), "")
    plotfreqs!(note.(4/3*freqs), "")
    plotfreqs!(note.(freqs*4), "Fourth with common overtone")
end

# ╔═╡ 06b01905-7249-4492-a5b2-b9b43f1908c1
md"Other harmonic intervals work by similar rules.
Each has an integer ratio, where the interval going up is the dual of the interval going down. When treating octaves as equivalent intervals can be normalized by shifting them by octaves into the range $[1, 2]$.
Each interval $m/n$ in $[1,2]$ has an inverse also in $[1,2]$: $2(n/m)$.
An interval and it's 'octave inverse' combine to give an octave
and they will have similar qualities."

# ╔═╡ Cell order:
# ╠═b7433dc1-846b-4173-ad51-f7598cd263a1
# ╠═37f59728-aeae-4b15-ba8c-a90368dfa8f0
# ╠═cfc3b26d-d2da-4128-9ef2-a9494c4c25e5
# ╠═b70c60fd-b29a-46f3-8daf-e6d7277156ca
# ╟─baafd0a2-b449-11eb-0654-172c64b81fd6
# ╠═7a4f0f06-6b37-4a3c-81e4-8004462da33a
# ╠═ef073029-a63e-4625-aecf-0cde96af5686
# ╟─9997ac3d-3cae-4bb5-8290-6ca53f431533
# ╟─dcd90ce9-dfc0-42ed-990a-55710c32b34e
# ╟─7f363c48-a194-4494-bb59-d5e00db70bd2
# ╠═9e45b6c3-a674-4b1f-b93e-11b4a99a22da
# ╟─07312f0d-59d9-45ee-b60a-c6ecfc2bc0df
# ╟─bfd2d6b6-2863-41ee-9b6a-c99ed19f8d7d
# ╟─3813605c-e5c8-4bd3-9813-f553d22d47b2
# ╠═6a831a92-a222-48e8-abfb-72bc9f3e0fb9
# ╟─0b9c89bf-291e-4b93-8632-f2d5117d25a4
# ╠═a4882ae0-32ed-4db4-ba0c-165081128f2d
# ╠═e3faa985-c01a-471c-ae10-0aaaee8a6b43
# ╠═6049ed44-b19f-41a9-b017-860e80d62017
# ╟─ccb62252-fcf1-4fc8-b486-1df69ea0c1e3
# ╠═14af9192-4b8a-49c3-9775-b4dfef343780
# ╠═02a98eb5-4021-4ea0-a0f8-b43c811d87b2
# ╟─31b44424-b3f9-4ea0-8d99-28738f1fc857
# ╠═82582e9b-fcc0-49ae-9f3b-558e70d8f0b8
# ╟─db7c4951-972b-43db-a647-6e61c8ebdf03
# ╟─53773643-1b58-4eee-ac11-bb30b7535fec
# ╟─151b3663-524e-465f-9426-3f85689e9e27
# ╟─24e5caba-ce1c-4436-8ca2-647b5438658a
# ╟─8adb4c0a-a29e-4cee-bb49-ee285f90ff70
# ╟─8378910d-3a2c-441e-a3b3-716d819eced9
# ╟─732ab799-970e-4ae7-a81d-cda4c79873f4
# ╠═df7e7e47-fe34-4250-bb0b-670fee0b56de
# ╟─518f1a8b-977b-4008-b1b7-667f0475b42c
# ╠═849e0964-2c80-4fa6-9320-5fa0c595d79a
# ╟─06b01905-7249-4492-a5b2-b9b43f1908c1

max_stars_repo_name: jnha/music

max_issues_repo_name: jnha/music

max_forks_repo_name: jnha/music


text:
include("../imports.jl")

tbbackend() = TensorBoardBackend(mktempdir())

@testset ExtendedTestSet "`LogMetrics`" begin
    cb = LogMetrics(tbbackend())
    learner = testlearner(Metrics(accuracy), Recorder(), cb)
    @test_nowarn fit!(learner, 1)
end

@testset ExtendedTestSet "`LogHyperParams`" begin
    cb = LogHyperParams(tbbackend())
    learner = testlearner(Recorder(), Scheduler(), cb)
    @test_nowarn fit!(learner, 1)
end

@testset ExtendedTestSet "`LogHistograms`" begin
    cb = LogHistograms(tbbackend(), freq = 5)
    learner = testlearner(Recorder(), Scheduler(), cb)
    @test_nowarn fit!(learner, 5)
end

@testset ExtendedTestSet "`LogVisualization`" begin
    cb = LogVisualization(tbbackend(), freq = 5) do batch
        return rand(RGB, 50, 50)
    end
    learner = testlearner(Recorder(), Scheduler(), cb)
    @test_nowarn fit!(learner, 5)
end


text:
module Ramnet

include("Utils.jl")

include("Encoders/Encoders.jl")

include("Partitioners/Partitioners.jl")

using Ramnet.Partitioners
export partition, RandomPartitioner, random_partitioning

include("Models/Models.jl")
using Ramnet.Models
export train!, predict

export Discriminator, BitDiscriminator, RegressionDiscriminator, GeneralizedRegressionDiscriminator

export MultiDiscriminatorClassifier

end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GDB_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GDB")
JLLWrappers.@generate_main_file("GDB", UUID("a8b33d9f-0f8b-5197-9986-f85cbaa50c6c"))
end  # module GDB_jll


text:
#Посчитать число всех горизонтальных прямолинейных перегородок (вертикальных - нет)
include("library.jl")
function count_gorizont_obstacles(r)
    back_path = BackPath(r, (Sud, West))
    print(counterPartition!(r))
    BackPath(r, (Sud, West))
    back!(r, back_path)
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
"""
    smssvd(X, d, stdThresholds=10 .^ range(-2,stop=0,length=100); nbrIter=10, maxSignalDim=typemax(Int)) -> (U,Σ,V,ps,signalDimensions,selectedVariables)

Computes the SubMatrix Selection Singular Value Decomposition (SMSSVD) of a matrix X.

**Inputs**

* `X`: PxN data matrix (P variables and N samples).
* `d`: Number of dimension in the result.
* `stdThresholds`: Standard deviation filtering thresholds at which the Projection Score will be evaluated. Expressed as fractions of the maximum standard deviation of a variable in X. Should be a vector of increasing values between 0 and 1.

**Outputs**

* `U`: Pxd matrix of left singular vectors (variable representation).
* `S`: d-vector of singular values.
* `V`: Nxd matrix of right singular vectors (sample representation).
* `ps`: d x length(stdThresholds) matrix with all Projection Scores. Useful for plotting.
* `signalDimensions`: Vector with the number of dimensions in each signal detected by SMSSVD.
* `selectedVariables`: For each signal, a bitmask showing the variables selected by Projection Score.
"""
function smssvd(X, d::Integer, stdThresholds=10 .^ range(-2,stop=0,length=100); nbrIter=10, maxSignalDim=typemax(Int))
    σMax = maximum(std(X,dims=2)) # Always base the variable filtering on the original σ's

    U = zeros(size(X,1),d)
    Σ = zeros(d)
    V = zeros(size(X,2),d)

    ps = zeros(d, length(stdThresholds))

    signalDimensions = Int[]
    selectedVariables = Vector{BitArray}()

    k = 1
    while k<=d
        # Filter variables and choose dimension by optimizing over Projection Score
        σ = dropdims(std(X,dims=2),dims=2) / σMax # always keep the same scale

        dmax = min(d-k+1, maxSignalDim)
        PS = projectionscorefiltered(X, 1:dmax, stdThresholds, nbrIter=nbrIter, σ=σ)
        dims,σInd = Tuple(CartesianIndices(size(PS))[argmax(PS)]) # Use projection score both to choose dimension and threshold for this signal
        σThresh = stdThresholds[σInd]

        r = k:k+dims-1
        ps[r,:] = PS[1:dims,:]
        push!(signalDimensions, dims)

        push!(selectedVariables,σ.>=σThresh)
        Y = X[selectedVariables[end],:]

        # Find the subspace Π that we are interested in.
        K = Symmetric(Y'Y)
        M = size(K,1)
        _,Π = eigen(K, M-dims+1:M) # only get the largest eigenvalues and vectors

        # Project X onto the subspace v and compute SVD. For dims=1, this is identical to uσ:=Xv.
        F = svd(X*Π) # solve for smaller matrix expressed in the basis of the subspace Π
        UΠ,ΣΠ,VΠ = F.U,F.S,F.V

        U[:,r] = UΠ
        Σ[r]   = ΣΠ
        V[:,r] = Π*VΠ # expand to original basis

        # make matrix orthogonal to previous component
        X = X - UΠ*(UΠ'X) # multiplication order avoiding big matrices

        k += dims
    end

    U,Σ,V,ps,signalDimensions,selectedVariables
end

# convenience method useful when calling smssvd from R
function smssvd(X, d::Vector{T}, stdThresholds=10 .^ range(-2,stop=0,length=100); kwargs...) where T<:Integer
    @assert length(d)==1
    smssvd(X, d[1], stdThresholds; kwargs...)
end


text:
Base.@kwdef struct Test
    a::Float64 = 1E12
    b::Float64 = pi
end


text:
# This module contains the necessary methods to use Fermi Arrays within TensorOperations

import TensorOperations: contract!
import TensorOperations: IndexTuple
import TensorOperations: IndexError
import TensorOperations: similarstructure_from_indices
import TensorOperations: scalar
import TBLIS

function tblis_set_num_threads(N)
    TBLIS.set_num_threads(N)
end

function set_num_threads(N = Threads.nthreads())
    tblis_set_num_threads(N)
    LinearAlgebra.BLAS.set_num_threads(N)
end

function oind2eins(oindA::NTuple{NAo}, cindA::NTuple{NAc},
          oindB::NTuple{NBo}, cindB::NTuple{NBc},
          tindC::NTuple{NCt}) where {NAo, NAc, NBo, NBc, NCt}

    # This function converts the numerical scheme used in TensorOperations
    # to represent contractions into an Einsum scheme.

    # Check contraction consistency.
    NAo + NBo == NCt || throw(IndexError("number of outer index not consistent."))
    NAc == NBc || throw(IndexError("number of contracted index not consistent."))

    cPadding = 'a' - 'A'
    einA = zeros(Int8, NAo+NAc)
    einB = zeros(Int8, NBo+NBc)

    # Outer indices.
    for i = 1:NAo
        einA[oindA[i]] = i
    end
    for i = 1:NBo
        einB[oindB[i]] = i + NAo
    end

    # Contracted indices.
    for i = 1:NAc
        einA[cindA[i]] = i + cPadding
        einB[cindB[i]] = i + cPadding
    end

    einA = string((einA .+'A')...)
    einB = string((einB .+'A')...)
    einC = string((tindC.+'A')...) # C has direct conversion relations.
    einA, einB, einC
end

function contract!(α,
          A::FermiMDArray, conjA::Symbol,
          B::FermiMDArray, conjB::Symbol,
          β,
          C::FermiMDArray,
          oindA::IndexTuple, cindA::IndexTuple,
          oindB::IndexTuple, cindB::IndexTuple,
          tindC::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    if Fermi.Options.get("tblis")
        # Check permutation consistency.
        # This check is copied from stridedarray.jl
        pA = (oindA...,cindA...)
        (length(pA) == ndims(A) && isperm(pA)) ||
            throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))
        pB = (oindB...,cindB...)
        (length(pB) == ndims(B) && isperm(pB)) ||
            throw(IndexError("invalid permutation of length $(ndims(B)): $pB"))
        (length(tindC) == ndims(C) && isperm(tindC)) ||
            throw(IndexError("invalid permutation of length $(ndims(C)): $tindC"))

        # Complex conjugate is not supported by TBLIS yet.
        # Here we do a intermediate copying-conversion.
        if conjA == :C || conjA == :A || conjB == :C || conjB == :A
            throw(Fermi.Options.FermiException("conj operations not suported for TBLIS"))
        end
        einA, einB, einC = oind2eins(oindA, cindA,
                                     oindB, cindB,
                                     tindC)

        CType = eltype(C.data)
        _A = TBLIS.TTensor{CType}(A.data,CType(α))
        _B = TBLIS.TTensor{CType}(B.data)
        _C = TBLIS.TTensor{CType}(C.data,CType(β))

        TBLIS.mul!(_C, _A, _B, einA, einB, einC)
        return C
    else
        contract!(α, A.data, conjA, B.data, conjB, β, C.data, oindA, cindA, oindB, cindB, tindC)
        return C
    end
end

function contract!(α,
          A::FermiMDArray, conjA::Symbol,
          B::FermiMDArray, conjB::Symbol,
          β,
          C::AbstractArray,
          oindA::IndexTuple, cindA::IndexTuple,
          oindB::IndexTuple, cindB::IndexTuple,
          tindC::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    contract!(α, A, conjA, B, conjB, β, FermiMDArray(C), oindA, cindA, oindB, cindB, tindC, syms)
end

function contract!(α,
          A::FermiMDArray, conjA::Symbol,
          B::AbstractArray, conjB::Symbol,
          β,
          C::AbstractArray,
          oindA::IndexTuple, cindA::IndexTuple,
          oindB::IndexTuple, cindB::IndexTuple,
          tindC::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    contract!(α, A, conjA, FermiMDArray(B), conjB, β, FermiMDArray(C), oindA, cindA, oindB, cindB, tindC, syms)
end

function contract!(α,
          A::AbstractArray, conjA::Symbol,
          B::FermiMDArray, conjB::Symbol,
          β,
          C::AbstractArray,
          oindA::IndexTuple, cindA::IndexTuple,
          oindB::IndexTuple, cindB::IndexTuple,
          tindC::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    contract!(α, FermiMDArray(A), conjA, B, conjB, β, FermiMDArray(C), oindA, cindA, oindB, cindB, tindC, syms)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing
using LightXML

@testset "create new Xdmf object" begin
    r = Xdmf()
    expected = "<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\" Version=\"2.1\"/>"
    @test string(r.xml) == expected
end

@testset "put and get to Xdmf, low level" begin
    io = Xdmf()
    # h5
    write(io.hdf, "/Xdmf/Domain/Geometry", [1 2 3])
    @test isapprox(read(io.hdf, "/Xdmf/Domain/Geometry"), [1 2 3])
    # xml
    obj = new_child(io.xml, "Domain")
    set_attribute(obj, "Name", "Test Domain")
    obj2 = find_element(io.xml, "Domain")
    @test attribute(obj2, "Name") == "Test Domain"
end

@testset "put and get to xdmf" begin
    xdmf = Xdmf()
    domain = new_child(xdmf, "Domain")
    grid = new_child(domain, "Grid")
    grid["CollectionType"] = "Temporal"
    grid["GridType"] = "Collection"

    frame1 = new_child(grid, "Grid")
    new_child(frame1, "Time", "Value" => 0.0)
    X1 = new_child(frame1, "Geometry", Dict("Type" => "XY"))

    frame2 = new_child(grid, "Grid", "Name" => "Frame 2")
    new_child(frame2, "Time", "Value" => 1.0)
    X2 = new_child(frame2, "Geometry", "Type" => "XY")

    add_child(grid, frame1)
    add_child(grid, frame2)

    dataitem = new_dataitem(xdmf, "/Domain/Grid/Grid/2/Geometry", [1.0, 2.0])
    add_child(X2, dataitem)

    println(xdmf.xml)
    @test has_child(xdmf.xml, "Domain")
    @test isa(get_child(xdmf.xml, "Domain"), XMLElement)
    @test !has_attribute(xdmf.xml, "Domain")
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Time/Value"), 0.0)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Time/Value"), 1.0)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[@Name=Frame 2]/Time/Value"), 1.0)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Geometry/DataItem"), [1.0, 2.0])
end

@testset "higher level xdmf" begin
    e1 = Element(Quad4, 1, [1, 2, 3, 4])
    e2 = Element(Quad4, 2, [5, 6, 7, 8])
    p1 = Problem(Elasticity, "Body 1", 2)
    p2 = Problem(Elasticity, "Body 2", 2)
    push!(p1, e1)
    push!(p2, e2)
    #update!(p1)

    e3 = Element(Seg2, 3, [1, 2])
    e4 = Element(Seg2, 4, [3, 4])
    e5 = Element(Seg2, 5, [5, 6])
    p3 = Problem(Dirichlet, "Fixed BC", 2, "displacement")
    p4 = Problem(Contact, "Contact between bodies 1 and 2", 2, "displacement")
    push!(p3, e3)
    push!(p4, e4)
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0],
        5 => [0.0, 2.0],
        6 => [1.0, 2.0],
        7 => [1.0, 3.0],
        8 => [0.0, 3.0])
    u = Dict{Int64, Vector{Float64}}(
        1 => [0.1, 0.1],
        2 => [0.1, 0.1],
        3 => [0.1, 0.1],
        4 => [0.1, 0.1],
        5 => [0.1, 0.1],
        6 => [0.1, 0.1],
        7 => [0.1, 0.1],
        8 => [0.1, 0.1])
    n = Dict{Int64, Vector{Float64}}(
        3 => [0.0, 1.0],
        4 => [0.0, 1.0])
    R = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 1.0],
        2 => [0.0, 1.0])
    update!(e4, "master elements", [e3])
end

#=

@testset "save results to disk, linear solver" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [0.0,0.0],
        2 => [1.0,0.0],
        3 => [1.0,1.0],
        4 => [0.0,1.0])
    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", X)
    update!(element, "temperature thermal conductivity", 6.0)
    update!(element, "temperature load", 0.0 => 12.0)
    update!(element, "temperature load", 1.0 => 24.0)
    problem = Problem(Heat, "one element heat problem", 1)
    problem.properties.formulation = "2D"
    push!(problem, element)
    boundary_element = Element(Seg2, [1, 2])
    update!(boundary_element, "geometry", X)
    update!(boundary_element, "temperature 1", 0.0)
    bc = Problem(Dirichlet, "fixed", 1, "temperature")
    push!(bc, boundary_element)
    xdmf = Xdmf()
    solver = Solver(Linear, problem, bc)
    solver.xdmf = xdmf

    solver.time = 0.0
    solver()
    empty!(problem.assembly)
    solver.time = 1.0
    solver()

    info(solver("temperature", 0.0))
    info(solver("temperature", 1.0))
    info(element("temperature load", [0.0, 0.0], 0.0))
    info(element("temperature load", [0.0, 0.0], 1.0))

    info("h5 file = $(h5file(xdmf))")
    E = read(xdmf.hdf, "/Topology/Quad4/Element IDs")
    C = read(xdmf.hdf, "/Topology/Quad4/Connectivity")
    N = read(xdmf.hdf, "/Node IDs")
    X = read(xdmf.hdf, "/Geometry")
    T1 = read(xdmf.hdf, "/Results/Time 0.0/Nodal Fields/Temperature")
    T2 = read(xdmf.hdf, "/Results/Time 1.0/Nodal Fields/Temperature")
    @test isapprox(E, [-1])
    @test isapprox(C, [0 1 2 3])
    @test isapprox(N, [1, 2, 3, 4])
    X_expected = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]'
    T1_expected = [0.0 0.0 1.0 1.0]
    T2_expected = [0.0 0.0 2.0 2.0]
    @test isapprox(X, X_expected)
    @test isapprox(T1, T1_expected)
    @test isapprox(T2, T2_expected)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Time/Value"), 0.0)
    @test read(xdmf, "/Domain/Grid/Grid/Geometry/Type") == "XY"
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Geometry/DataItem"), X_expected)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology/DataItem"), [0 1 2 3])
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology[@TopologyType=Polyline]/DataItem"), [0 1])
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[1]/Attribute[@Name=Temperature]/DataItem"), T1_expected)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Attribute[@Name=Temperature]/DataItem"), T2_expected)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Topology/DataItem"), [0 1 2 3])
end

@testset "save results to disk, nonlinear solver" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [0.0,0.0],
        2 => [1.0,0.0],
        3 => [1.0,1.0],
        4 => [0.0,1.0])
    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", X)
    update!(element, "temperature thermal conductivity", 6.0)
    update!(element, "temperature load", 0.0 => 12.0)
    update!(element, "temperature load", 1.0 => 24.0)
    problem = Problem(Heat, "one element heat problem", 1)
    problem.properties.formulation = "2D"
    push!(problem, element)
    boundary_element = Element(Seg2, [1, 2])
    update!(boundary_element, "geometry", X)
    update!(boundary_element, "temperature 1", 0.0)
    bc = Problem(Dirichlet, "fixed", 1, "temperature")
    push!(bc, boundary_element)
    solver = Solver(Nonlinear, problem, bc)
    solver.xdmf = Xdmf()

    solver.time = 0.0
    solver()
    solver.time = 1.0
    solver()

    xdmf = get(solver.xdmf)
    info("h5 file = $(h5file(xdmf))")
    E = read(xdmf.hdf, "/Topology/Quad4/Element IDs")
    C = read(xdmf.hdf, "/Topology/Quad4/Connectivity")
    N = read(xdmf.hdf, "/Node IDs")
    X = read(xdmf.hdf, "/Geometry")
    T11 = read(xdmf.hdf, "/Results/Time 0.0/Iteration 1/Nodal Fields/Temperature")
    T12 = read(xdmf.hdf, "/Results/Time 0.0/Iteration 2/Nodal Fields/Temperature")
    T21 = read(xdmf.hdf, "/Results/Time 1.0/Iteration 1/Nodal Fields/Temperature")
    T22 = read(xdmf.hdf, "/Results/Time 1.0/Iteration 2/Nodal Fields/Temperature")
    X_expected = [
        0.0 0.0
        1.0 0.0
        1.0 1.0
        0.0 1.0]
    T1_expected = [0.0 0.0 1.0 1.0]
    T2_expected = [0.0 0.0 2.0 2.0]
    @test isapprox(T12, T1_expected)
    @test isapprox(T22, T2_expected)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Time/Value"), 0.0)
    @test read(xdmf, "/Domain/Grid/Grid/Geometry/Type") == "XY"
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Geometry/DataItem"), X_expected')
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology/DataItem"), [0 1 2 3])
    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology[@TopologyType=Polyline]/DataItem"), [0 1])
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[1]/Attribute[@Name=Temperature]/DataItem"), T1_expected)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Attribute[@Name=Temperature]/DataItem"), T2_expected)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0)
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Topology/DataItem"), [0 1 2 3])
end

=#


text:
Base.init_depot_path()
Base.init_load_path()

@eval Module() begin
    Base.include(@__MODULE__, "CSV2Parquet.jl")
    for (pkgid, mod) in Base.loaded_modules
        if !(pkgid.name in ("Main", "Core", "Base"))
            eval(@__MODULE__, :(const $(Symbol(mod)) = $mod))
        end
    end
    for statement in readlines("app_precompile.jl")
        try
            Base.include_string(@__MODULE__, statement)
        catch
            # See julia issue #28808
            Core.println("failed to compile statement: ", statement)
        end
    end
end # module

empty!(LOAD_PATH)
empty!(DEPOT_PATH)

max_stars_repo_name: xiaodaigh/csv-to-parquet

max_issues_repo_name: xiaodaigh/csv-to-parquet

max_forks_repo_name: xiaodaigh/csv-to-parquet


text:
using Random, DataFrames, DataFramesMeta, CSV

# Create jitter data of horizontal axis in Figure 4 and Figure 7
function jitter_csv(;
                    # Nalgos = 11, Ncombinations = 30,
                    Nalgos = 9, Ncombinations = 42,
                    savepath = homedir(),
                    # filestr="gauss"
                    filestr="multi"
                    )
   df = DataFrame()
   from = 0.
   for i in 1:Nalgos
       df[Symbol("x"*string(i))] = rand(Ncombinations) .+ from
       from += 1.5
    end
    CSV.write(joinpath(savepath, "jitter_"* filestr *".csv"),
        df, delim = " ")
end


text:
module RSUtils

export supnorm,
       fixedpoint,
       analytic_histogram,
       sequential_compare,
       is_increasing,
       newton

using Distributions: Distributions

include("fixedpoint.jl")
include("optimization.jl")
include("statistics.jl")
include("predicates.jl")

end #module


text:
# This file is a part of RunStatistics.jl, licensed under the MIT License (MIT).

export squares_pvalue, squares_cdf

const log_factorial = Vector{Float64}(undef, 0)

function cachechi2(T_obs::Real, N::Int)

    @argcheck 0 < N

    res = zeros(N + 1)
    res[1] = NaN

    for i = firstindex(res) + 1 :lastindex(res)

        res[i] = logcdf(Chisq(i - 1), T_obs)
    end

    return res
end

"""
    squares_cdf(T_obs::Real, N::Integer)

Compute P(T < T_obs | N), the value of the cumulative distribution of the Squares statistic `T` at the
value `T_obs` observed in `N` independent trials with gaussian probability.

`T_obs` is the value of the test statistic for the observed data set; i.e., the largest `χ^2` of
any run of consecutive observed values above the expectation in a sequence of `N` independent trials
with Gaussian uncertainty.

`N` is the total number of data points.

The calculation implements equations (16) and (17) from

Frederik Beaujean and Allen Caldwell. *A Test Statistic for
Weighted Runs.* Journal of Statistical Planning and Inference 141,
no. 11 (November 2011): 3437–46. doi:10.1016/j.jspi.2011.04.022

https://arxiv.org/abs/1005.3233.

"""
function squares_cdf(T_obs::Real, N::Integer)

    T = promote_type(typeof(T_obs), typeof(N))

    if N >= length(log_factorial)

        sizehint!(log_factorial, N)

        if isempty(log_factorial)
            push!(log_factorial, 0.0)
        end

        for i = length(log_factorial):N
            push!(log_factorial, last(log_factorial) + log(i))
        end

    end

    log_cumulative = cachechi2(T_obs, N)

    logpow2N1 = (N <= (8 * sizeof(Int) - 1)) ? log((1 << N) - 1) : N * log(2)

    p = zero(T)

    for r = 1:N

        Mmax = min(r, N - r + 1)
        poch = 0.0

        for M = 1:Mmax

            poch += log(N - r + 2 - M)
            scale = poch - logpow2N1
            ppi = 0.0

            g = init_partition(r, M)
            n = g.c
            y = g.y

            done = false

            while !done

                h = g.h
                ppartition = 0.0

                for l = 2:(h + 1)
                    ppartition += n[l] * log_cumulative[y[l]+1] - log_factorial[n[l]+1]
                end

                ppi += exp(ppartition)

                done = next_partition!(g)
            end

            p += exp(scale + log(ppi))
        end
    end

    @assert p < 1
    return p
end

"""
    squares_pvalue(T_obs::Real, N::Integer)

Compute P(T >= `T_obs` | `N`), the p value for the Squares test statistic `T` being larger or equal to `T_obs`,
the value of the Squares statistic observed in `N` datapoints.

The Squares statistic `T` denotes the largest `χ^2` of any run of consecutive successes (above expectation) in a
sequence of `N` independent trials with Gaussian uncertainty.

Via `squares_cdf()` this function implements equations (16) and (17) from

Frederik Beaujean and Allen Caldwell. *A Test Statistic for
Weighted Runs.* Journal of Statistical Planning and Inference 141,
no. 11 (November 2011): 3437–46. doi:10.1016/j.jspi.2011.04.022

https://arxiv.org/abs/1005.3233.

"""
function squares_pvalue(T_obs::Real, N::Integer)

    return 1 - squares_cdf(T_obs, N)

end


text:
# This file includes mapping of interpolation interval to different subintervals.

using FractalTools
using Plots

# Construct interpolation data
xi, dx, xf = 0, 0.1, 1
x = collect(xi : dx : xf)
y = sin.(2π * x)
d = 0.01 * ones(length(x) - 1)

# Interpolate the data
interp = interpolate(x, y, d)

# Get transformations
ws = interp.ifs.ws

# Get an example transformation
pts = hcat.(eachcol([x'; y']))

# Plots
plt = plot(layout=(2,1))
for (i, w) in enumerate(ws)
    v = getindex.(w.(pts), 1)
    plot!(v, ones(length(v)) * i, label="", subplot=1)
end
for (i, w) in enumerate(ws)
    v = getindex.(w.(pts), 2)
    plot!(ones(length(v)) * (i - 1), v, label="", subplot=2)
end
display(plt)


text:
"""
```
hessizero{T<:AbstractFloat}(fcn::Function, x::Vector{T};
                            check_neg_diag::Bool=false,
                            verbose::Symbol=:none,
                            distr::Bool=true)
```

Compute Hessian of function `fcn` evaluated at `x`.

### Arguments
- `check_neg_diag`: Throw an error if any negative diagonal elements are detected.
- `verbose`: Print verbose output
- `distr`: Use available parallel workers to increase performance.
"""
function hessizero{T<:AbstractFloat}(fcn::Function,
                                    x::Vector{T};
                                    check_neg_diag::Bool=false,
                                    verbose::Symbol=:none,
                                    distr::Bool=true)
    n_para = length(x)
    hessian  = zeros(n_para, n_para)

    # Compute diagonal elements first
    if distr && nworkers() > 1
        diag_elements = @sync @parallel (hcat) for i = 1:n_para
            hess_diag_element(fcn, x, i; check_neg_diag=check_neg_diag, verbose=verbose)
        end
        for i = 1:n_para
            hessian[i, i] = diag_elements[i]
        end
    else
        for i=1:n_para
            hessian[i,i] = hess_diag_element(fcn, x, i; check_neg_diag=check_neg_diag,
                                             verbose=verbose)
        end
    end

    # Now compute off-diagonal elements
    # Make sure that correlations are between -1 and 1
    # invalid_corr indexes elements that are invalid
    invalid_corr = Dict{Tuple{Int,Int}, Float64}()

    # Build indices to iterate over
    n_off_diag_els = Int(n_para*(n_para-1)/2)
    off_diag_inds = Vector{Tuple{Int,Int}}(n_off_diag_els)
    k=1
    for i=1:(n_para-1), j=(i+1):n_para
        off_diag_inds[k] = (i,j)
        k = k+1
    end

    # Iterate over off diag elements
    if distr
        off_diag_out = @sync @parallel (hcat) for (i,j) in off_diag_inds
            σ_xσ_y = sqrt(abs(hessian[i, i]*hessian[j, j]))
            hess_offdiag_element(fcn, x, i, j, σ_xσ_y; verbose=verbose)
        end
        # Ensure off_diag_out is array
        off_diag_out = hcat(off_diag_out)
    else
        off_diag_out = Array{Tuple{T, T},1}(n_off_diag_els)
        for (k,(i,j)) in enumerate(off_diag_inds)
            σ_xσ_y = sqrt(abs(hessian[i, i]*hessian[j, j]))
            off_diag_out[k] = hess_offdiag_element(fcn, x, i, j, σ_xσ_y; verbose=verbose)
        end
    end

    # Fill in values
    for k=1:n_off_diag_els
        (i,j) = off_diag_inds[k]
        (value, ρ_xy) = off_diag_out[k]

        hessian[i,j] = value
        hessian[j,i] = value

        if ρ_xy < -1 || 1 < ρ_xy
            invalid_corr[(i, j)] = ρ_xy
        end
    end

    has_errors = false
    if !isempty(invalid_corr)
        println("Errors: $invalid_corr")
        has_errors = true
    end

    return hessian, has_errors
end

# Compute diag element
function hess_diag_element{T<:AbstractFloat}(fcn::Function,
                                              x::Vector{T},
                                              i::Int;
                                              ndx::Int=6,
                                              check_neg_diag::Bool=false,
                                              verbose::Symbol=:none)
    # Setup
    n_para = length(x)
    dxscale  = ones(n_para, 1)
    dx       = exp.(-(6:2:(6+(ndx-1)*2))')
    hessdiag = zeros(ndx, 1)

    # Computation
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Hessian element: ($i, $i)")
    end

    # Diagonal element computation
    for k = 3:4
        paradx    = copy(x)
        parady    = copy(x)
        paradx[i] = paradx[i] + dx[k]*dxscale[i]
        parady[i] = parady[i] - dx[k]*dxscale[i]

        fx  = fcn(x)
        fdx = fcn(paradx)
        fdy = fcn(parady)

        hessdiag[k]  = -(2fx - fdx - fdy) / (dx[k]*dxscale[i])^2
    end

    if VERBOSITY[verbose] >= VERBOSITY[:high]
        println("Values: $(hessdiag)")
    end

    value = (hessdiag[3]+hessdiag[4])/2

    if check_neg_diag && value < 0
        error("Negative diagonal in Hessian")
    end

    if VERBOSITY[verbose] >= VERBOSITY[:high]
        println("Value used: $value")
    end

    return value
end

# Compute off diag element
function hess_offdiag_element{T<:AbstractFloat}(fcn::Function,
                                                 x::Vector{T},
                                                 i::Int,
                                                 j::Int,
                                                 σ_xσ_y::T;
                                                 ndx::Int=6,
                                                 verbose::Symbol=:none)
    # Setup
    n_para = length(x)
    dxscale  = ones(n_para, 1)
    dx       = exp.(-(6:2:(6+(ndx-1)*2))')
    hessdiag = zeros(ndx, 1)

    # Computation
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Hessian element: ($i, $j)")
    end

    for k = 3:4
        paradx      = copy(x)
        parady      = copy(x)
        paradx[i]   = paradx[i] + dx[k]*dxscale[i]
        parady[j]   = parady[j] - dx[k]*dxscale[j]
        paradxdy    = copy(paradx)
        paradxdy[j] = paradxdy[j] - dx[k]*dxscale[j]

        fx    = fcn(x)
        fdx   = fcn(paradx)
        fdy   = fcn(parady)
        fdxdy = fcn(paradxdy)

        hessdiag[k]  = -(fx - fdx - fdy + fdxdy) / (dx[k]*dx[k]*dxscale[i]*dxscale[j])
    end

    if VERBOSITY[verbose] >= VERBOSITY[:high]
        println("Values: $(hessdiag)")
    end

    value = (hessdiag[3]+hessdiag[4])/2

    if value == 0 || σ_xσ_y == 0
        ρ_xy = 0
    else
        ρ_xy = value / σ_xσ_y
    end

    if ρ_xy < -1 || 1 < ρ_xy
        value = 0
    end

    if VERBOSITY[verbose] >= VERBOSITY[:high]
        println("Value used: $value")
        println("Correlation: $ρ_xy")
    end

    return value, ρ_xy
end


text:
#= Instructions:

- Pkg.add("PkgBenchmark.jl")
- using PkgBenchmark

- results = benchmarkpkg("IntervalArithmetic")
- show(results)

To compare the current version with a given release use `judge`:

- results = judge("IntervalArithmetic", tag)
- show(results)

To export the benchmark results to a markdown file, do:

- export_markdown("results.md", results)
=#

using BenchmarkTools, IntervalArithmetic

const SUITE = BenchmarkGroup() # parent BenchmarkGroup to contain our suite

S = SUITE["Constructors"] = BenchmarkGroup()
S["Interval"] = @benchmarkable Interval(1, 2)

S = SUITE["Intervals"] = BenchmarkGroup()
S["@interval"] = @benchmarkable @interval(0.1)
S["pi"] = @benchmarkable @interval(pi)
S["expression"] = @benchmarkable @interval(sin(0.1) + cos(0.2))

S = SUITE["Arithmetic"] = BenchmarkGroup()
begin
    a = Interval(1, 2)
    b = Interval(3, 4)

    for op in (+, -, *, /)
        S[string(op)] = @benchmarkable $(op)($a, $b)
    end
end

S = SUITE["Elementary functions"] = BenchmarkGroup()
begin
    for op in (exp, log, sin, tan)
        S[string(op)] = @benchmarkable $(op)($a)
    end
end

S = SUITE["Sum"] = BenchmarkGroup()
begin
    sum1(N) = sum(Interval(i, i+1) for i in 1:N)
    sum2(N) = (one = Interval(1.0); sum(one / (i^2) for i in 1:N))
    S["Sum1"] = @benchmarkable sum1(1000)
    S["Sum2"] = @benchmarkable sum2(1000)
end


text:
using ClimateMachine

ClimateMachine.init()

using ClimateMachine.CartesianFields: RectangularElement, assemble

@testset "Cartesian fields" begin

    Nx = 3
    Ny = 4
    Nz = 5

    data = rand(Nx, Ny, Nz)
    x = repeat(range(0.0, 1.0, length = Nx), 1, Ny, Nz)
    y = repeat(range(0.0, 1.1, length = Ny), Nx, 1, Nz)
    z = repeat(range(0.0, 1.2, length = Nz), Nx, Ny, 1)

    element = RectangularElement(data, x, y, z)

    @test size(element) == (Nx, Ny, Nz)
    @test maximum(element) == maximum(data)
    @test minimum(element) == minimum(data)
    @test maximum(abs, element) == maximum(abs, data)
    @test minimum(abs, element) == minimum(abs, data)
    @test element[1, 1, 1] == data[1, 1, 1]

    east_element = RectangularElement(2 .* data, x .+ x[end, 1, 1], y, z)
    north_element = RectangularElement(3 .* data, x, y .+ y[1, end, 1], z)
    top_element = RectangularElement(4 .* data, x, y, z .+ z[1, 1, end])

    west_east = assemble(Val(1), element, east_element)
    south_north = assemble(Val(2), element, north_element)
    bottom_top = assemble(Val(3), element, top_element)

    @test west_east.x[1, 1, 1] == x[1, 1, 1]
    @test west_east.x[end, 1, 1] == 2 * x[end, 1, 1]

    @test south_north.y[1, 1, 1] == y[1, 1, 1]
    @test south_north.y[1, end, 1] == 2 * y[1, end, 1]

    @test bottom_top.z[1, 1, 1] == z[1, 1, 1]
    @test bottom_top.z[1, 1, end] == 2 * z[1, 1, end]

    northeast_element =
        RectangularElement(5 .* data, x .+ x[end, 1, 1], y .+ y[1, end, 1], z)

    # Remember that matrix literals are transposed, so "northeast"
    # is the bottom right corner (for example).
    four_elements = [
        element north_element
        east_element northeast_element
    ]

    four_elements = reshape(four_elements, 2, 2, 1)

    four_way = assemble(four_elements)

    @test four_way.x[end, 1, 1] == west_east.x[end, 1, 1]
    @test four_way.y[1, end, 1] == south_north.y[1, end, 1]
    @test four_way.z[1, 1, end] == z[1, 1, end]
end


text:
function main(r::Robot)
    path = go_to_left_bottom_corner_and_get_path(r)
    for i ∈ (Ost, Nord, West, Sud)
        mark_direction(r, i)
    end
    go_back(r, path)
end

function go_to_left_bottom_corner_and_get_path(r::Robot) ::Array
    path=[]
    while !isborder(r,Sud) || isborder(r,West)
        if !isborder(r,Sud)
            move!(r,Sud)
            push!(path,Nord)
        end
        if !isborder(r,West)
            move!(r,West)
            push!(path,Ost)
        end
    end
    return path
end

function mark_direction(r,side::HorizonSide)
    while !isborder(r,side)
        putmarker!(r)
        move!(r,side)
    end
    putmarker!(r)
end

function go_back(r,path::Array)
    n=length(path)
    while n>0
        move!(r,path[n])
        n=n-1
    end
end

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
using Distributed
if nprocs() == 1
    addprocs(3)
end
# Number of logical CPU cores available in the system.
# println("CPU cores: ", Sys.CPU_CORES)
Sys.cpu_summary()
# Number of available processes.
# println("nprocs: ", nprocs())
# println("nworkers: ", nworkers())
println("nprocs: ", nprocs())
println("nworkers: ", nworkers())

n = 20000000
# Ordinary for loop.
let
    global sum1
    sum1 = 0
    @time for i = 1:n
         sum1 += Int(rand(Bool));
    end
end

# Parallel for loop
# @time @parallel (+) for i = 1:20000
@time sum2 = @distributed (+) for i = 1:n
    Int(rand(Bool));
end

# Predefined function.
@time sum3 =sum(rand(0:1, n))

# For a not so small amount of work:
# n = 200000000;
# n = 20000;
# @time @parallel (+) for i = 1:n
@time sum4 = @distributed (+) for i = 1:n
    Int(rand(Bool))
end

# println(sum1)
# println(sum2)
# println(sum3)
# println(sum4)

# DO NOT run sum(rand(0:1, n)) or the ordinary for loop with 'n': You may run
# out of memory.


text:
import Voronoi

using Random
using Printf
include("Draw.jl")

function randf(start, finish, n)
   v = rand(n)
   return map(x -> start + x*(finish-start), v)
end

function gradientDescent()
   WIDTH = 100.0
   HEIGHT = 100.0

   seed = 100
   println("Random seed: ", seed)
   Random.seed!(seed)

   command = nothing
   for n in 2:100
      points = convert(Array{Tuple{Real, Real}}, collect(zip(randf(1, WIDTH-1, n), randf(1, HEIGHT-1, n))))

      iterations = 1
      fortunes = 1

      time = @elapsed V, f, ∇f, ξ = Voronoi.Optimization.init(points, WIDTH, HEIGHT)

      μ = 0.5
      α = 0.01
      ϵ = 10^-1
      λ₀ = 0.01

      times = time

      while ξ > ϵ
         @printf("f = %.10g\t\t|∇f| = %.10g\t\ttime: %s seconds\n", f, ξ, time)
         Draw.init(WIDTH, HEIGHT)
         Draw.voronoiDiagram(V)

         d = ∇f

         time = @elapsed V, f, ∇f, ξ, forts = Voronoi.Optimization.lineSearch(points, d, α, λ₀, μ, ϵ, V, f, ∇f, ξ)

         times += time
         fortunes += forts
         iterations += 1

         if command != "a"
            println("Press Return for a step or enter \"a\" to animate until the end.")
            command = readline(stdin)
         else
            sleep(0.0000001) # make sure stuff is drawn
         end
      end

      println("n: ", n, "\t\titer: ", iterations, "\t\tfortunes: ", fortunes, "\t\ttime: ", times)
      #Draw.init(WIDTH, HEIGHT)
      #Draw.voronoiDiagram(V)
      #command = readline(stdin)
   end
end

gradientDescent()

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
module TestBlobsCompat

using Blobs
using Test

# basic sanity checks

struct Foo
    x::Int64
    y::Float32 # tests conversion and alignment
end

foo = Blobs.malloc_and_init(Foo)
@v foo.x = 1
@test @v(foo.x) == 1
@v foo.y = 2.5
@test (@v foo.y) == 2.5
@test (@v foo) == Foo(1,2.5)
# tests interior pointers
@test (@a foo) == foo
@test pointer(@a foo.y) == pointer(foo) + sizeof(Int64)

# tests ambiguous syntax
@test_throws LoadError eval(:(@v foo.y == 2.5))
@test_throws LoadError eval(:(@a foo.y == 2.5))

mpv = Blobs.malloc_and_init(BlobVector{Foo}, 3)
pv = @v mpv
pv[2] = Foo(2, 2.2)
@test pv[2] == Foo(2, 2.2)
@test length(pv) == 3
pv[1] = Foo(1, 1.1)
pv[3] = Foo(3, 3.3)
# tests iteration
@test collect(pv) == [Foo(1,1.1), Foo(2,2.2), Foo(3,3.3)]
# tests interior pointers
@test pointer(@a mpv[2]) - pointer(pv.data) == Blobs.self_size(Foo)

bbv = Blobs.malloc_and_init(BlobBitVector, 3)
pbv = @v bbv
pbv[2] = true
@test pbv[2] == true
@test pv[2] == Foo(2, 2.2)
pbv[1] = false
pbv[3] = false
# tests iteration
@test collect(pbv) == [false, true, false]
fill!(pbv, false)
@test collect(pbv) == [false, false, false]

# strings and unicode

s = "普通话/普通話"
mp = Blobs.malloc_and_init(BlobString, s)
p = @v mp
@test p == s
@test repr(p) == repr(s)
@test collect(p) == collect(s)
@test findfirst(p, "通") == findfirst(s, "通")

# test right-to-left

s = "سلام"
mp = Blobs.malloc_and_init(BlobString, s)
p = @v mp
@test p == s
@test repr(p) == repr(s)
@test collect(p) == collect(s)
@test findfirst(p, "ا") == findfirst(s, "ا")

# test string conversions

@test isbits(p)
@test String(p) isa String
@test String(p) == s

# test string slices

s = "aye bee sea"
mp = Blobs.malloc_and_init(BlobString, s)
p = @v mp
@test p[5:7] isa BlobString
@test p[5:7] == "bee"

# sketch of paged pmas

struct PackedMemoryArray{K,V}
     keys::BlobVector{K}
     values::BlobVector{V}
     mask::BlobBitVector
     count::Int
     #...other stuff
end

function Blobs.child_size(::Type{PackedMemoryArray{K,V}}, length::Int64) where {K,V}
    T = PackedMemoryArray{K,V}
    +(Blobs.child_size(fieldtype(T, :keys), length),
      Blobs.child_size(fieldtype(T, :values), length),
      Blobs.child_size(fieldtype(T, :mask), length))
  end

function Blobs.init(pma::Blob{PackedMemoryArray{K,V}}, free::Blob{Nothing}, length::Int64) where {K,V}
    free = Blobs.init(pma.keys, free, length)
    free = Blobs.init(pma.values, free, length)
    free = Blobs.init(pma.mask, free, length)
    fill!(pma.mask[], false)
    pma.count[] = 0
    free
end

pma = Blobs.malloc_and_init(PackedMemoryArray{Int64, Float32}, 3)
@test (@v pma.count) == 0
@test (@v pma.keys.length) == 3
# tests fill!
@test !any(@v pma.mask)
# tests pointer <-> offset conversion
@test unsafe_load(convert(Ptr{UInt64}, pointer(pma)), 1) == Blobs.self_size(PackedMemoryArray{Int64, Float32})
# tests nested interior pointers
pma2 = pma.mask[2]
@test pma2[] == false
pma2[] = true
@test pma2[] == true
@test pma.mask[2][] == true

# test api works ok with varying sizes

struct Quux
    x::BlobVector{Int}
    y::Float64
end

struct Bar
    a::Int
    b::BlobBitVector
    c::Bool
    d::BlobVector{Float64}
    e::Blob{Quux}
end

function Blobs.child_size(::Type{Quux}, x_len::Int64, y::Float64)
    T = Quux
    +(Blobs.child_size(fieldtype(T, :x), x_len))
end

function Blobs.child_size(::Type{Bar}, b_len::Int64, c::Bool, d_len::Int64, x_len::Int64, y::Float64)
    T = Bar
    +(Blobs.child_size(fieldtype(T, :b), b_len),
      Blobs.child_size(fieldtype(T, :d), d_len),
      Blobs.child_size(fieldtype(T, :e), x_len, y))
end

function Blobs.init(quux::Blob{Quux}, free::Blob{Nothing}, x_len::Int64, y::Float64)
    free = Blobs.init(quux.x, free, x_len)
    quux.y[] = y
    free
end

function Blobs.init(bar::Blob{Bar}, free::Blob{Nothing}, b_len::Int64, c::Bool, d_len::Int64, x_len::Int64, y::Float64)
    free = Blobs.init(bar.b, free, b_len)
    free = Blobs.init(bar.d, free, d_len)
    free = Blobs.init(bar.e, free, x_len, y)
    bar.c[] = c
    free
end

bar = Blobs.malloc_and_init(Bar, 10, false, 20, 15, 1.5)
quux = @v bar.e

@test length(@v bar.b) == 10
@test (@v bar.c) == false
@test length(@v bar.d) == 20
@test length(@v quux.x) == 15
@test (@v quux.y) == 1.5

end


text:
using DataStructures
using Base.Meta

import Base.ImmutableDict

const EClassId = Int64

# struct ENode{T, M}
mutable struct ENode{T}
    head::Any
    args::Vector{EClassId}
    # metadata::M
    # the nodes where it came from
    # proof_src::Vector{Tuple{Rule, ENode, Int}}
    # the enodes that were generated by this enode
    # proof_trg::Vector{Tuple{Rule, ENode, Int}}
    # age of the egraph when this enode was created
    # age::Int
    hash::Ref{UInt} # hash cache
end

# function ENode{T}(head, c_ids::AbstractVector{EClassId}, ps=[], pt=[], age=0) where {T}
function ENode{T}(head, c_ids::AbstractVector{EClassId}) where {T}
        # static_args = MVector{length(c_ids), Int64}(c_ids...)
    # m = metadata(e)
    # ENode{T}(head, c_ids, ps, pt, age, Ref{UInt}(0))
    ENode{T}(head, c_ids, Ref{UInt}(0))
end

ENode(a) = ENode{typeof(a)}(a, EClassId[])

ENode(a::ENode) =
    error("constructor of ENode called on enode. This should never happen")

function Base.:(==)(a::ENode, b::ENode)
    isequal(a.args, b.args) &&
    isequal(a.head, b.head)
end

# ===============================================================
# proof functions
# ===============================================================

# function setage!(n::ENode, age::Int)
#     n.age = age
# end

# function addproofsrc!(n::ENode{T}, rule::Rule, src::ENode, age::Int) where {T}
#     push!(n.proof_src, (rule, src, age))
# end

# function addprooftrg!(n::ENode{T}, rule::Rule, trg::ENode, age::Int) where {T}
#     push!(n.proof_trg, (rule, trg, age))
# end

# function hasproofdata(n::ENode)
#     !(isempty(n.proof_src) && isempty(n.proof_trg))
# end

# function mergeproof!(target::ENode, source::ENode)
#     append!(target.proof_src, source.proof_src)
#     append!(target.proof_trg, source.proof_trg)
# end

# This optimization comes from SymbolicUtils
# The hash of an enode is cached to avoid recomputing it.
# Shaves off a lot of time in accessing dictionaries with ENodes as keys.
function Base.hash(t::ENode{T}, salt::UInt) where {T}
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    # h′ = hash(t.args,  hash(t.metadata, hash(t.head, hash(T, salt))))
    h′ = hash(t.args,  hash(t.head, hash(T, salt)))
    t.hash[] = h′
    return h′
end

TermInterface.arity(n::ENode) = length(n.args)
# TermInterface.metadata(n::ENode) = n.metadata

termtype(x::ENode{T}) where T = T

function Base.show(io::IO, x::ENode{T}) where {T}
    print(io, "ENode{$T}(", x.head)
    n = arity(x)
    if n == 0
        print(io, ")")
        return
    else
        print(io, " ")
    end
    for i ∈ 1:n
        if i < n
            print(io, x.args[i], " ")
        else
            print(io, x.args[i], ")")
        end
    end
end


text:
log!(::Any, ::Any, ::Any, ::Any) = nothing
mutable struct BasicLogger{F}
    fbest::Vector{Float64}
    xbest::Vector{Vector{Float64}}
    fmedian::Vector{Float64}
    frange::Vector{Float64}
    times::Vector{Float64}
    callback::F
    verbosity::Int
end
function Base.show(io::IO, ::MIME"text/plain", s::BasicLogger)
    println(io, "BasicLogger")
end
function BasicLogger(x0; verbosity = 1,
                     callback = (o, y, fvals, perm) -> nothing)
    BasicLogger(Float64[], Vector{Float64}[], Float64[], Float64[], Float64[],
                callback, verbosity)
end
function log!(o::Optimizer{<:Any, <:BasicLogger}, y, fvals, perm)
    l = o.logger
    if length(l.fbest) == 0 || fvals[perm[1]] < l.fbest[end]
        push!(l.fbest, fvals[perm[1]])
        push!(l.xbest, o.p.mean + sigma(o.p) * y[:, perm[1]])
    end
    push!(l.fmedian, median(fvals))
    push!(l.frange, maximum(fvals) - minimum(fvals))
    l.callback(o, y, fvals, perm)
    l
end
best(l::BasicLogger) = l.fbest[end]


text:
############################################################################################
"""
$(TYPEDEF)

Utility struct to help assign boundaries to parameter.

# Fields
$(TYPEDFIELDS)
"""
struct Simplex <: AbstractConstraint
    len::Int64
    function Simplex(len::Integer)
        ArgCheck.@argcheck len > 0
        return new(len)
    end
end
Simplex(vec::AbstractVector) = Simplex(length(vec))

############################################################################################
function construct_transform(constraint::Simplex, val)
    transform = Bijectors.SimplexBijector{1, true}()
    return transform, Bijectors.inverse(transform)
end

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::Simplex,
    val::Vector{R},
) where {R<:Real}
    ArgCheck.@argcheck length(val) == constraint.len
    ArgCheck.@argcheck all(val[iter] > 0.0 for iter in eachindex(val))
    return true
end

############################################################################################
#= !NOTES:
    (1) Constrained's last element will sum up to 1.
    (2) Unconstrained's last element is irrelevant for Bijector{Simplex} in default method.
    (3) Consequently, we can flatten in length(x)-1 dimensions, and unflatten back to length(x) by summing up elements for length(x)'s element.
    This works for both constrained/unconstrained.
=#
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenStrict,
    constraint::C,
    x::Vector{R},
) where {
    T<:AbstractFloat,
    F<:FlattenTypes,
    R<:Real,
    C<:Union{Simplex, Distributions.Dirichlet, Bijectors.SimplexBijector}
}
    buffer_flat = zeros(T, length(x)-1)
    len_flat = length(x)
    len_unflat = len_flat-1
    function flatten_Simplex(x_vec::AbstractVector{R}) where {R<:Real}
        ArgCheck.@argcheck length(x_vec) == len_flat
        return fill_array!(buffer_flat, view(x_vec, 1:len_flat-1))
    end
    buffer_unflat = zeros(R, length(x))
    function unflatten_Simplex(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len_unflat
        return Simplex_from_flatten!(buffer_unflat, v)
    end
    return flatten_Simplex, unflatten_Simplex
end

function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenFlexible,
    constraint::C,
    x::Vector{R},
) where {
    T<:AbstractFloat,
    F<:FlattenTypes,
    R<:Real,
    C<:Union{Simplex, Distributions.Dirichlet, Bijectors.SimplexBijector}
}
    len_flat = length(x)
    len_unflat = len_flat-1
    function flatten_Simplex_AD(x_vec::AbstractVector{R}) where {R<:Real}
        ArgCheck.@argcheck length(x_vec) == len_flat
        buffer = zeros(R, len_flat-1)
        return fill_array!(buffer, view(x_vec, 1:len_flat-1))
    end
    function unflatten_Simplex_AD(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len_unflat
        buffer = zeros(eltype(v), length(x))
        return Simplex_from_flatten!(buffer, v)
    end
    return flatten_Simplex_AD, unflatten_Simplex_AD
end

############################################################################################
#=
!NOTE: For this constraint, we use a Bijector as Transformer, so we do not need to add a new functors
1. MyTransformer <: AbstractTransformer, MyInverseTransformer <: AbstractTransformer
2. define a function construct_transform(MyConstraint, val) -> MyTransformer, MyInverseTransformer
3. overload unconstrain and log_abs_det_jac on MyTransformer.
4. overload constrain on MyInverseTransformer.
=#

############################################################################################
#Export
export
    Simplex,
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/../" "${BASH_SOURCE[0]}" "$@" -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

module FaceDetection

import Base: size, getindex, LinearIndices
using Images: Images, coords_spatial

export to_integral_image, sum_region
export learn, get_feature_votes
export feature_types, HaarLikeObject, get_score, get_vote
export displaymatrix, notify_user, filtered_ls, load_image,
    ensemble_vote_all, reconstruct, get_random_image,
    generate_validation_image, get_faceness, determine_feature_size

include("HaarLikeFeature.jl")
include("Utils.jl") # Utils.jl exports HaarLikeFeature.jl functions
include("IntegralImage.jl")
include("AdaBoost.jl")

end # end module