text:
using RvSpectML
 using Statistics
 using Dates

make_plots = false
include("neid_solar_1_read.jl")
 order_list_timeseries = RvSpectML.make_order_list_timeseries(solar_data)
 order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
 lambda_range_with_data = (min = maximum(d->minimum(d.λ),solar_data), max = minimum(d->maximum(d.λ),solar_data) )
 RvSpectML.normalize_spectra!(order_list_timeseries,solar_data);

espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks","G2.espresso.mas")
 espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
 line_list_df = espresso_df |>
     @filter(lambda_range_with_data.min <= _.lambda ) |>
    @filter( _.lambda < lambda_range_with_data.max) |>
#    @filter( _.lambda < 6000.0 ) |>                       # Avoid tellurics at redder wavelengths
#    @filter( _.lambda >6157 || _.lambda < 6155  ) |>   # Avoid "line" w/ large variability
    DataFrame

# Setup to run CCF
mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=1.6)
 line_list = RvSpectML.CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
 ccf_plan = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list)
 v_grid = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan)

# Compute CCF's & measure RVs
tstart = now()
 @time ccfs = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
 println("# CCF runtime: ", now()-tstart)

make_plots = true
 rvs_ccf_gauss = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i],fit_type = "gaussian") for i in 1:length(order_list_timeseries) ]
 rvs_ccf_quad  = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i], fit_type = "quadratic") for i in 1:length(order_list_timeseries) ]
 rvs_ccf_cent  = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i], fit_type = "centroid") for i in 1:length(order_list_timeseries) ]
 rvs_ccf_best  = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i], fit_type = "bestfit") for i in 1:length(order_list_timeseries) ]
 println("RMS of RVs (ESPRESSO lines):  Gaussian: ",std(rvs_ccf_gauss), "   Quadratic: ", std(rvs_ccf_quad), "   Centroid: ", std(rvs_ccf_cent) )
 if make_plots
   using Plots
  plot(v_grid,ccfs,label=:none)
  xlabel!("v (m/s)")
  ylabel!("CCF")
 end

make_plots = true
 if make_plots
   t_idx = 1
   using Plots
   plot(v_grid,ccfs[:,t_idx],label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
 end

#make_plots = true
if make_plots
   nbin = 4
   plt_t = (order_list_timeseries.times .- minimum(order_list_timeseries.times) ) .* 24
   plt = scatter(plt_t,rvs_ccf_gauss.-mean(rvs_ccf_gauss),label="Gaussian")
   scatter!(plt,plt_t,rvs_ccf_quad.-mean(rvs_ccf_quad),label="Quadratic")
   #scatter!(plt,plt_t,rvs_ccf_best.-mean(rvs_ccf_best),label="Best fit")
   #scatter!(plt,plt_t,rvs_ccf_cent.mean(rvs_ccf_cent),label="Centroid")
   times_binned = RvSpectML.bin_times_consecutive(plt_t, nbin)
   rvs_ccf_gauss_binned = RvSpectML.bin_rvs_consecutive(rvs_ccf_gauss, nbin)
   rvs_ccf_quad_binned = RvSpectML.bin_rvs_consecutive(rvs_ccf_quad, nbin)
   rvs_ccf_cent_binned = RvSpectML.bin_rvs_consecutive(rvs_ccf_cent, nbin)
   println("# RMS of binned (n=",nbin,") RVs: ", std(rvs_ccf_gauss_binned), "  ", std(rvs_ccf_quad_binned), "  ", std(rvs_ccf_cent_binned))
   scatter!(times_binned,rvs_ccf_gauss_binned.-mean(rvs_ccf_gauss),label="Binned Gauss")
   scatter!(times_binned,rvs_ccf_quad_binned.-mean(rvs_ccf_quad),label="Binned Quad")
   #scatter!(times_binned,rvs_ccf_cent_binned.-mean(rvs_ccf_cent),label="Binned Cent")
   plot!(times_binned,rvs_ccf_gauss_binned.-mean(rvs_ccf_gauss),label=:none,color=3)
   plot!(times_binned,rvs_ccf_quad_binned.-mean(rvs_ccf_quad),label=:none,color=4)
   #plot!(times_binned,rvs_ccf_cent_binned.-mean(rvs_ccf_cent),label=:none,color=5)
   xlabel!(plt,"Time (hr)")
   ylabel!(plt,"RV (m/s)")
   display(plt)
 end

#=
# TODO: Need to update code below
# Try using Alex's clean lines
vald_filename = joinpath(ancilary_solar_data_path,"VALD_Fe1_DP_rejectTelluricSlope0.0_badLineFilterESPRESSO-strict-NEID-BIS_overlapcutoff6e-05_depthcutoff0.05_allowBlends0_wavesReiners_depthssolar_nbin1depth0.mas")
 vald_df = RvSpectML.read_mask_vald(vald_filename)
 ssn_out = RvSpectML.searchsortednearest(line_list_df.lambda,vald_df.lambda)
 line_list_clean = line_list_df[ssn_out,:]
 @assert maximum(line_list_clean.lambda .- vald_df.lambda) < 0.05
 clean_mask = hcat(line_list_clean.lambda./mask_entry_doppler_factor,line_list_clean.lambda.*mask_entry_doppler_factor, line_list_clean.weight)

 @time ccfs = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, line_list, mask_shape=mask_shape, plan=ccf_plan)
  rvs_ccf_gauss = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i],fit_type = "gaussian") for i in 1:length(order_list_timeseries) ]
  rvs_ccf_quad  = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i], fit_type = "quadratic") for i in 1:length(order_list_timeseries) ]
  rvs_ccf_cent  = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i], fit_type = "centroid") for i in 1:length(order_list_timeseries) ]
  rvs_ccf_best  = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i], fit_type = "bestfit") for i in 1:length(order_list_timeseries) ]
  println("RMS of RVs (Alex's lines):  Gaussian: ",std(rvs_ccf_gauss), "   Quadratic: ", std(rvs_ccf_quad), "   Centroid: ", std(rvs_ccf_cent) )
  if make_plots
    using Plots
   plot(v_grid,ccfs,label=:none)
   xlabel!("Δv (m/s)")
   ylabel!("CCF")
  end

=#


text:
module FakePackage

using InlineTest

greet() = "Hello World!"

const RUN = []

@testset "FakePackage inline begin-end" begin
    @test greet() == "Hello World!"
    @testset "inner" begin
        push!(RUN, 1)
        @test true
    end
end

@testset "FakePackage inline for" for i=1:1
    @test greet() == "Hello World!"
    push!(RUN, 2)
end

module Sub1

using InlineTest
using ..FakePackage: RUN

@testset "check that Sub1 is noticed by ReTest" begin
    push!(RUN, 3)
    @test true
end

end # Sub1

module Sub2
# Sub2 itself doesn't have any @testset, but has a submodule which has some

module SubSub

using InlineTest
using ...FakePackage: RUN

@testset "check that Sub2.SubSub is noticed by ReTest" begin
    push!(RUN, 4)
    @test true
end

end # SubSub
end # Sub2

end # FakePackage

ext: jl


text:
using Test
include("funnel.jl")

# analytical vs autodiff.
@test (l(ones(2)), ∇l(ones(2))) == dldx(ones(2))

max_stars_repo_name: cmerkatas/HMC

max_issues_repo_name: cmerkatas/HMC

max_forks_repo_name: cmerkatas/HMC


text: abstract type AbstractEstuary end


text:
using CellularAutomata

const size_space = 6
const gens = 5

glider = [[0, 0, 1, 0, 0] [0, 0, 0, 1, 0] [0, 1, 1, 1, 0]]
space = zeros(Bool, size_space, size_space)
space[1:size(glider, 1), 1:size(glider, 2)] = glider

ca = CellularAutomaton(Life((3, (2,3))), space, gens)

@test ca.evolution == cat([0 0 0 0 0 0; 0 0 1 0 0 0; 1 0 1 0 0 0; 0 1 1 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0],
[0 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 1 0 0; 0 1 1 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0],
[0 0 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 1 1 1 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0],
[0 0 0 0 0 0; 0 0 0 0 0 0; 0 1 0 1 0 0; 0 0 1 1 0 0; 0 0 1 0 0 0; 0 0 0 0 0 0],
[0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 1 0 0; 0 1 0 1 0 0; 0 0 1 1 0 0; 0 0 0 0 0 0], dims=3)


text:
# results

term_status_map = Dict(
    SCIP_STATUS_UNKNOWN => MOI.OPTIMIZE_NOT_CALLED,
    SCIP_STATUS_USERINTERRUPT => MOI.INTERRUPTED,
    SCIP_STATUS_NODELIMIT => MOI.NODE_LIMIT,
    SCIP_STATUS_TOTALNODELIMIT => MOI.NODE_LIMIT,
    SCIP_STATUS_STALLNODELIMIT => MOI.OTHER_LIMIT,
    SCIP_STATUS_TIMELIMIT => MOI.TIME_LIMIT,
    SCIP_STATUS_MEMLIMIT => MOI.MEMORY_LIMIT,
    SCIP_STATUS_GAPLIMIT => MOI.OPTIMAL,
    SCIP_STATUS_SOLLIMIT => MOI.SOLUTION_LIMIT,
    SCIP_STATUS_BESTSOLLIMIT => MOI.OTHER_LIMIT,
    SCIP_STATUS_RESTARTLIMIT => MOI.OTHER_LIMIT,
    SCIP_STATUS_OPTIMAL => MOI.OPTIMAL,
    SCIP_STATUS_INFEASIBLE => MOI.INFEASIBLE,
    SCIP_STATUS_UNBOUNDED => MOI.DUAL_INFEASIBLE,
    SCIP_STATUS_INFORUNBD => MOI.INFEASIBLE_OR_UNBOUNDED,
    SCIP_STATUS_TERMINATE => MOI.INTERRUPTED,
)

function MOI.get(o::Optimizer, ::MOI.TerminationStatus)
    return term_status_map[SCIPgetStatus(o)]
end

function MOI.get(o::Optimizer, attr::MOI.PrimalStatus)
    return if 1 <= attr.result_index <= MOI.get(o, MOI.ResultCount())
        MOI.FEASIBLE_POINT
    else
        MOI.NO_SOLUTION
    end
end

function MOI.get(::Optimizer, ::MOI.DualStatus)
    return MOI.NO_SOLUTION
end

function MOI.get(o::Optimizer, ::MOI.ResultCount)::Int
    status = SCIPgetStatus(o)
    if status in [SCIP_STATUS_UNBOUNDED, SCIP_STATUS_INFORUNBD]
        return 0
    end
    return SCIPgetNSols(o)
end

function MOI.get(o::Optimizer, ::MOI.RawStatusString)
    return string(SCIPgetStatus(o))
end

"Make sure that SCIP is currently in one of the allowed stages."
function assert_stage(o::Optimizer, stages)
    stage = SCIPgetStage(o)
    if !(stage in stages)
        error("SCIP is wrong stage ($stage, need $stages)!")
    end
end

"Make sure that the problem was solved (SCIP is in SOLVED stage)."
function assert_solved(o::Optimizer)
    # SCIP's stage is SOLVING when stopped by user limit!
    assert_stage(o, (SCIP_STAGE_PRESOLVING, SCIP_STAGE_SOLVING, SCIP_STAGE_PRESOLVED, SCIP_STAGE_SOLVED))

    # Check for invalid status (when stage is SOLVING).
    status = SCIPgetStatus(o)
    if status in (SCIP_STATUS_UNKNOWN,
                  SCIP_STATUS_USERINTERRUPT,
                  SCIP_STATUS_TERMINATE)
        error("SCIP's solving was interrupted, but not by a user-given limit!")
    end
end

"Make sure that: TRANSFORMED ≤ stage ≤ SOLVED."
assert_after_prob(o::Optimizer) = assert_stage(o, SCIP_Stage.(3:10))

function MOI.get(o::Optimizer, attr::MOI.ObjectiveValue)
    assert_solved(o)
    MOI.check_result_index_bounds(o, attr)
    sols = unsafe_wrap(Array{Ptr{SCIP_SOL}}, SCIPgetSols(o), SCIPgetNSols(o))
    return SCIPgetSolOrigObj(o, sols[attr.result_index])
end

function MOI.get(o::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    assert_solved(o)
    MOI.check_result_index_bounds(o, attr)
    sols = unsafe_wrap(Array{Ptr{SCIP_SOL}}, SCIPgetSols(o), SCIPgetNSols(o))
    return SCIPgetSolVal(o, sols[attr.result_index], var(o, vi))
end

function MOI.get(o::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{VI,<:BOUNDS})
    assert_solved(o)
    MOI.check_result_index_bounds(o, attr)
    sols = unsafe_wrap(Array{Ptr{SCIP_SOL}}, SCIPgetSols(o), SCIPgetNSols(o))
    return SCIPgetSolVal(o, sols[attr.result_index], var(o, VI(ci.value)))
end

function MOI.get(o::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{<:SAF,<:BOUNDS})
    assert_solved(o)
    MOI.check_result_index_bounds(o, attr)
    sols = unsafe_wrap(Array{Ptr{SCIP_SOL}}, SCIPgetSols(o), SCIPgetNSols(o))
    return SCIPgetActivityLinear(o, cons(o, ci), sols[attr.result_index])
end

function MOI.get(o::Optimizer, ::MOI.ObjectiveBound)
    assert_after_prob(o)
    return SCIPgetDualbound(o)
end

function MOI.get(o::Optimizer, ::MOI.RelativeGap)
    assert_stage(o, [SCIP_STAGE_PRESOLVING, SCIP_STAGE_SOLVING, SCIP_STAGE_SOLVED])
    return SCIPgetGap(o)
end

function MOI.get(o::Optimizer, ::MOI.SolveTimeSec)
    return SCIPgetSolvingTime(o)
end

function MOI.get(o::Optimizer, ::MOI.SimplexIterations)
    assert_stage(o, [SCIP_STAGE_PRESOLVING, SCIP_STAGE_SOLVING, SCIP_STAGE_SOLVED])
    return SCIPgetNLPIterations(o)
end

function MOI.get(o::Optimizer, ::MOI.NodeCount)
    return SCIPgetNNodes(o)
end


text:
# Abstract type for Local Spatial Autocorrelation
abstract type AbstractLocalSpatialAutocorrelation end

# Function to build the conditional randomization sample and calculate the local scores
function crand_local(permutations::Int, z::AbstractVector{T} where T<:Real, W::SpatialWeights, local_calc_function::Function, rng::AbstractRNG)::Matrix{Float64}
    # Build conditional permutations array
    ni = cardinalities(W)
    n = length(ni)
    maxni = maximum(ni)
    Cperms = zeros(Int, permutations, maxni)
    samplevec = 1:n-1
    for i in 1:permutations
        Cperms[i,:] = sample(rng, samplevec, maxni, replace = false)
    end

    # Calculate LISA for all the permutations
    lisaperms = zeros(n, permutations)

    Threads.@threads for i in 1:n
        bnoi = ones(Bool, n)
        bnoi[i] = false

        zi = z[i]
        znoi = z[bnoi]

        nni = ni[i]
        wi = weights(W, i)

        for p in 1:permutations
            zcrand = view(znoi, Cperms[p, 1:nni])
            lisaperms[i, p] = local_calc_function(zi, wi, zcrand)
        end
    end

    return lisaperms
end

"""
    issignificant(x, α, adjust = :none)
Return a vector of boolean values indicating if the local statistics are significant or not at the desired threshold ``α``.

p-values can be adjusted with the `adjust` parameter using the Bonferroni correction `:bonferroni` or controlling for the False Discovery Rate, `:fdr`
"""
function issignificant(x::AbstractLocalSpatialAutocorrelation, α::Float64; adjust::Symbol = :none)::Vector{Bool}
    p = pvalue(x)
    n = length(p)

    (adjust == :none || adjust == :bonferroni || adjust == :fdr) ||
        throw(ArgumentError("unknown p-value adjustment $(adjust)"))

    if adjust == :none
        return p .< α
    elseif adjust == :bonferroni
        return p .< (α / n)
    elseif adjust == :fdr
        psort = sort(p)
        pfdr = (1:n) .* α ./ n
        lower = psort .< pfdr
        threshold = findfirst(x -> x == false, lower)
        return p .< pfdr[threshold]
    end
end

"""
    assignments(x, α, adjust = :none)
Return a vector with the categories assigned by the local statistic with a significance threshold ``α``.

p-values can be adjusted with the `adjust` parameter using the Bonferroni correction `:bonferroni` or controlling for the False Discovery Rate, `:fdr`
"""
function assignments(x::AbstractLocalSpatialAutocorrelation, α::Float64; adjust::Symbol = :none)::Vector{Symbol}
    q = deepcopy(assignments(x))
    p = pvalue(x)

    q[.! issignificant(x, α, adjust = adjust)] .= :ns

    return q
end

function Base.show(io::IO, x::AbstractLocalSpatialAutocorrelation)

    println(io, testname(x) * " test of Spatial Autocorrelation")
    println(io, "--------------------------------------------")
    println(io, "")
    println(io, "Randomization test with ", size(scoreperms(x), 2), " permutations.")
    println(io, "Interesting locations at 0.05 significance level:")

    # Display clusters
    q = assignments(x)
    p = pvalue(x)

    labels = labelsorder(x)
    nl = length(labels)
    countcat = zeros(Int, nl)
    for i in 1:nl
        issig = issignificant(x, 0.05, adjust = :none)
        if labels[i] == :ns
            countcat[i] = count(.! issig)
        else
            countcat[i] = count(issig .& (q .== labels[i]))
        end
    end

    labelsstr = labelsnames(x)
    labelmaxlength = maximum(length.(labelsstr))
    for i in 1:nl
        if labels[i] != :ns
            println(io, " ", lpad(labelsstr[i], labelmaxlength), ": ", countcat[i])
        end
    end

end


text:
module Registrator

using Base64
using LibGit2
using UUIDs

# Remove all of a base64 string's whitespace before decoding it.
decodeb64(s::AbstractString) = String(base64decode(replace(s, r"\s" => "")))

include("slack.jl")
include("pull_request.jl")
include("Messaging.jl")
include("RegService.jl")
include("commentbot/CommentBot.jl")
include("webui/WebUI.jl")

end # module


text:
@testset "Testing Batchnorm" begin
    @testset "Testing Forward Training 2D" begin
        _src = randn(Float32, 100, 100)
        _scale_shift = randn(Float32, size(_src, 1), 2)

        src = OneDNN.Memory(_src)
        scale_shift = OneDNN.Memory(_scale_shift)

        # Run the OneDNN model
        nt = OneDNN.batchnorm_training(src, scale_shift; epsilon = 1f-5)

        # Some sanity checks.
        @test isapprox(OneDNN.materialize(nt.mean), Statistics.mean.(eachrow(_src)))
        @test isapprox(
            OneDNN.materialize(nt.variance),
            Statistics.var.(eachrow(_src); corrected = false),
        )

        # Create the Flux.jl model for comparison purposes.
        bn_flux = Flux.BatchNorm(
            size(_src, 1); affine = true, track_stats = false, ϵ = 1f-5
        )

        bn_flux.γ .= view(_scale_shift, :, 1)
        bn_flux.β .= view(_scale_shift, :, 2)

        expected = bn_flux(_src)
        @test isapprox(OneDNN.materialize(nt.dst), expected)
    end

    @testset "Testing Forward Training 4D" begin
        _src = randn(Float32, 10, 10, 10, 10)
        _scale_shift = randn(Float32, size(_src, 3), 2)

        src = OneDNN.Memory(_src)
        scale_shift = OneDNN.Memory(_scale_shift)

        # Run the OneDNN model
        nt = OneDNN.batchnorm_training(src, scale_shift; epsilon = 1f-5)

        # Some sanity checks.
        dims = (1, 2, 4)
        @test isapprox(
            OneDNN.materialize(nt.mean), view(Statistics.mean(_src; dims), 1, 1, :, 1)
        )
        @test isapprox(
            OneDNN.materialize(nt.variance),
            view(Statistics.var(_src; dims, corrected = false), 1, 1, :, 1),
        )

        # Create the Flux.jl model for comparison purposes.
        bn_flux = Flux.BatchNorm(
            size(_src, 3); affine = true, track_stats = false, ϵ = 1f-5
        )

        bn_flux.γ .= view(_scale_shift, :, 1)
        bn_flux.β .= view(_scale_shift, :, 2)

        expected = bn_flux(_src)
        @test isapprox(OneDNN.materialize(nt.dst), expected)
    end

    @testset "Testing Forward and Backward Pair" begin
        _scale_shift = randn(Float32, 100, 2)
        _src = randn(Float32, 100, 100)
        epsilon = 1f-5

        scale_shift = OneDNN.Memory(_scale_shift)
        src = OneDNN.Memory(_src)

        activations = [identity, Flux.relu, Flux.sigmoid, abs]
        for fn in activations
            @show fn
            bn_onednn = OneDNN.BatchNorm(scale_shift, fn; epsilon = epsilon)

            bn_flux = Flux.BatchNorm(
                size(_src, 1), fn; affine = true, track_stats = false, ϵ = 1f-5
            )
            bn_flux.γ .= view(_scale_shift, :, 1)
            bn_flux.β .= view(_scale_shift, :, 2)

            # Test results and pullbacks
            out_onednn, back_onednn = Zygote._pullback(bn_onednn, src)
            out_flux, back_flux = Zygote._pullback(bn_flux, _src)

            @test isapprox(OneDNN.materialize(out_onednn), out_flux)

            diff_onednn = back_onednn(out_onednn)
            diff_flux = back_flux(out_flux)
            @test length(diff_onednn) == length(diff_flux) == 2

            # Compare `diff_src`.
            # TODO: Figure out why there's such a difference between OneDNN and Flux ...
            @test isapprox(OneDNN.materialize(diff_onednn[2]), diff_flux[2]; rtol = 0.21)

            # scale and shift
            diff_scale_shift = OneDNN.materialize(diff_onednn[1].scale_shift)
            @test isapprox(diff_flux[1][].γ, view(diff_scale_shift, :, 1))
            @test isapprox(diff_flux[1][].β, view(diff_scale_shift, :, 2))
        end
    end
end


text:
## Min-cut algorithms using max-flow subroutine

# Min-cut from bfs on classical max-flow
function mincut(g::Network, request = :EdmondsKarp)
  if request == :EdmondsKarp
    val, flows = edmondsKarp(g)
  end
#  println(flows)
  cut = bfs(g,:MinCut,flows)
  return val, cut
end

# α_k of a cut
function partSum(g::Network,cutedges::Vector{Int},k::Int)
  Σ = Vector{Float64}()
  for i in cutedges
    push!(Σ,g.links[i].cap)
  end
  sort!(Σ)
  return sum(Σ[1:end-k])
end

function bestMinkCut(g::Network,cuts::Vector{Vector{Int}},k::Int)
  function aux(cutedges)
    partSum(g,cutedges,k)
  end
  return minimum(map(aux,cuts))
end


text:
module TestEnsembles

# uncomment two lines for testing distributed processing
# using Distributed
# addprocs(2)
# @everywhere using Pkg
# @everywhere Pkg.activate("working", shared=true)

# using Revise
using Test
using Random
using MLJ
using MLJBase
import MLJModels
using CategoricalArrays
using Distributions

@load KNNRegressor

## WRAPPED ENSEMBLES OF FITRESULTS

# target is :deterministic :multiclass false:
atom = MLJModels.DeterministicConstantClassifier()
L = ['a', 'b', 'j']
L2 = categorical(L)
ensemble = [L2[1], L2[3], L2[3], L2[2]]
n=length(ensemble)
weights = fill(1/n, n) # ignored by predict below
wens = MLJ.WrappedEnsemble(atom, ensemble)
X = MLJ.table(rand(3,5))
@test predict(wens, weights, X) == categorical(vcat(['j','j','j'],L))[1:3]

# target is :deterministic :continuous false:
atom = MLJModels.DeterministicConstantRegressor()
ensemble = Float64[4, 7, 4, 4]
weights = [0.1, 0.5, 0.2, 0.2]
wens = MLJ.WrappedEnsemble(atom, ensemble)
@test predict(wens, weights, X) ≈ [5.5, 5.5, 5.5]

# target is :probabilistic :multiclass false:
atom = ConstantClassifier()
L = categorical(['a', 'b', 'j'])
d1 = UnivariateFinite(L, [0.1, 0.2, 0.7])
d2 = UnivariateFinite(L, [0.2, 0.3, 0.5])
ensemble = [d2,  d1, d2, d2]
weights = [0.1, 0.5, 0.2, 0.2]
wens = MLJ.WrappedEnsemble(atom, ensemble)
X = MLJ.table(rand(2,5))
d = predict(wens, weights, X)[1]
@test pdf(d, 'a') ≈ 0.15
@test pdf(d, 'b') ≈ 0.25
@test pdf(d, 'j') ≈ 0.6

# target is :probabilistic :continuous false:
atom = ConstantRegressor()
d1 = Distributions.Normal(1, 2)
d2 = Distributions.Normal(3, 4)
ensemble = [d2,  d1, d2, d2]
weights = [0.1, 0.5, 0.2, 0.2]
wens = MLJ.WrappedEnsemble(atom, ensemble)
X = MLJ.table(rand(2,5))
d = predict(wens, weights, X)[1]

## ENSEMBLE MODEL

# target is :deterministic :multiclass false:
atom=MLJModels.DeterministicConstantClassifier()
X = MLJ.table(ones(5,3))
y = categorical(collect("asdfa"))
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
weights = rand(10)
weights = weights/sum(weights)
ensemble_model.weights = weights
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)
@test MLJBase.target_scitype(ensemble_model) == MLJBase.target_scitype(atom)

# target is :deterministic :continuous false:
atom = MLJModels.DeterministicConstantRegressor()
X = MLJ.table(ones(5,3))
y = Float64[1.0, 2.0, 1.0, 1.0, 1.0]
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
@test reduce(* , [x ≈ 1.0 || x ≈ 1.25 for x in fitresult.ensemble])
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
ensemble_model.bagging_fraction = 1.0
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
@test unique(fitresult.ensemble) ≈ [1.2]
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
weights = rand(10)
weights = weights/sum(weights)
ensemble_model.weights = weights
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)

# target is :deterministic :continuous false:
atom = MLJModels.DeterministicConstantRegressor()
Random.seed!(1234)
X = MLJ.table(randn(10,3))
y = randn(10)
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom, rng=1)
ensemble_model.out_of_bag_measure = [MLJ.rms,MLJ.rmsp]
ensemble_model.n = 2
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
# TODO: the following test fails in distributed version (because of multiple rng's ?)
@test report[:oob_estimates][1] ≈ 1.083490899041915
# @test MLJBase.output_is(ensemble_model) == MLJBase.output_is(atom)
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom,rng=Random.MersenneTwister(1))
ensemble_model.out_of_bag_measure = MLJ.rms
ensemble_model.n = 2
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)

# target is :probabilistic :multiclass false:
atom = ConstantClassifier()
X = MLJ.table(ones(5,3))
y = categorical(collect("asdfa"))
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.ProbabilisticEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
fitresult.ensemble
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
ensemble_model.bagging_fraction = 1.0
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
fitresult.ensemble
d = predict(ensemble_model, fitresult, MLJ.selectrows(X, test))[1]
@test pdf(d, 'a') ≈ 2/5
@test pdf(d, 's') ≈ 1/5
@test pdf(d, 'd') ≈ 1/5
@test pdf(d, 'f') ≈ 1/5
weights = rand(10)
weights = weights/sum(weights)
ensemble_model.weights = weights
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)
# @test MLJBase.output_is(ensemble_model) == MLJBase.output_is(atom)

# target is :probabilistic :continuous false:
atom = ConstantRegressor()
X = MLJ.table(ones(5,3))
y = Float64[1.0, 2.0, 2.0, 1.0, 1.0]
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.ProbabilisticEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
d1 = fit(Distributions.Normal, [1,1,2,2])
d2 = fit(Distributions.Normal, [1,1,1,2])
# @test reduce(* , [d.μ ≈ d1.μ || d.μ ≈ d2.μ for d in fitresult.ensemble])
# @test reduce(* , [d.σ ≈ d1.σ || d.σ ≈ d2.σ for d in fitresult.ensemble])
d=predict(ensemble_model, fitresult, MLJ.selectrows(X, test))[1]
for dc in d.components
    @test pdf(dc, 1.52) ≈ pdf(d1, 1.52) || pdf(dc, 1.52) ≈ pdf(d2, 1.52)
end
ensemble_model.bagging_fraction = 1.0
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
d = predict(ensemble_model, fitresult, MLJ.selectrows(X, test))[1]
d3 = fit(Distributions.Normal, y)
@test pdf(d, 1.52) ≈ pdf(d3, 1.52)
weights = rand(10)
weights = weights/sum(weights)
ensemble_model.weights = weights
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)
# @test MLJBase.output_is(ensemble_model) == MLJBase.output_is(atom)

# test generic constructor:
@test EnsembleModel(atom=ConstantRegressor()) isa Probabilistic
@test EnsembleModel(atom=MLJModels.DeterministicConstantRegressor()) isa Deterministic

## MACHINE TEST

N =100
X = (x1=rand(N), x2=rand(N), x3=rand(N))
y = 2X.x1  - X.x2 + 0.05*rand(N)

atom = KNNRegressor(K=7)
ensemble_model = EnsembleModel(atom=atom)
ensemble = machine(ensemble_model, X, y)
train, test = partition(eachindex(y), 0.7)
fit!(ensemble, rows=train); length(ensemble.fitresult.ensemble)
ensemble_model.n = 15
fit!(ensemble);
@test length(ensemble.fitresult.ensemble) == 15
ensemble_model.n = 10
fit!(ensemble);
@test length(ensemble.fitresult.ensemble) == 10
@test !isnan(predict(ensemble, MLJ.selectrows(X, test))[1])

# old Koala tests
# # check that providing fixed seed gives identical predictions each
# # time if trees are deterministic:
# tree.max_features = 0
# ensemble_model.rng_seed = 1234
# fit!(ensemble, train);
# p1 = predict(ensemble, X, test[1:10])
# fit!(ensemble, train);
# p2 = predict(ensemble, X, test[1:10])
# @test p1 == p2

# tree.max_features = 4
# fit!(ensemble);
# ensemble_model.weight_regularization = 0.5
# fit_weights!(ensemble);
# display(ensemble.report[falsermalized_weights])
# err(ensemble, test)
# fit!(ensemble, train, add=true);
# err(ensemble, test)
# u,v = weight_regularization_curve(ensemble, test, raw=false,
#                                   range=range(0, stop=1, length=21))
# UnicodePlots.lineplot(u,v)
# ensemble_model.weight_regularization = u[argmin(v)]
# fit_weights!(ensemble);
# err(ensemble, test)

end
true


text:
using IntervalLinearAlgebra
using Test

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


text:
using Gaugefields
using LinearAlgebra
using InteractiveUtils
using Random
using LatticeDiracOperators

function MDtest!(gauge_action,U,Dim,fermi_action,η,ξ)
    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients.
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 20
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0
    Random.seed!(123)

    numtrj = 2000
    for itrj = 1:numtrj
        @time accepted = MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,fermi_action,η,ξ)
        numaccepted += ifelse(accepted,1,0)

        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        println("$itrj plaq_t = $plaq_t")
        println("acceptance ratio ",numaccepted/itrj)
    end
end

function calc_action(gauge_action,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(gauge_action,U)/NC #evaluate_GaugeAction(gauge_action,U) = tr(evaluate_GaugeAction_untraced(gauge_action,U))
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end

function MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,fermi_action,η,ξ)
    Δτ = 1/MDsteps
    NC,_,NN... = size(U[1])

    gauss_distribution!(p)

    substitute_U!(Uold,U)
    gauss_sampling_in_action!(ξ,U,fermi_action)
    sample_pseudofermions!(η,U,fermi_action,ξ)
    Sfold = real(dot(ξ,ξ))
    println("Sfold = $Sfold")

    Sold = calc_action(gauge_action,U,p) + Sfold
    println("Sold = ",Sold)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,gauge_action)

        P_update!(U,p,1.0,Δτ,Dim,gauge_action)
        P_update_fermion!(U,p,1.0,Δτ,Dim,gauge_action,fermi_action,η)

        U_update!(U,p,0.5,Δτ,Dim,gauge_action)
    end
    Sfnew = evaluate_FermiAction(fermi_action,U,η)
    println("Sfnew = $Sfnew")
    Snew = calc_action(gauge_action,U,p) + Sfnew

    println("Sold = $Sold, Snew = $Snew")
    println("Snew - Sold = $(Snew-Sold)")

    accept = exp(Sold - Snew) >= rand()

    #ratio = min(1,exp(Snew-Sold))
    if accept != true #rand() > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,gauge_action)
    temps = get_temporary_gaugefields(gauge_action)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)

    end
end

function P_update!(U,p,ϵ,Δτ,Dim,gauge_action) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temps = get_temporary_gaugefields(gauge_action)
    dSdUμ = temps[end]
    factor =  -ϵ*Δτ/(NC)

    for μ=1:Dim
        calc_dSdUμ!(dSdUμ,gauge_action,μ,U)
        mul!(temps[1],U[μ],dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temps[1])
    end
end

function P_update_fermion!(U,p,ϵ,Δτ,Dim,gauge_action,fermi_action,η)  # p -> p +factor*U*dSdUμ
    #NC = U[1].NC
    temps = get_temporary_gaugefields(gauge_action)
    UdSfdUμ = temps[1:Dim]
    factor =  -ϵ*Δτ

    calc_UdSfdU!(UdSfdUμ,fermi_action,U,η)

    for μ=1:Dim
        Traceless_antihermitian_add!(p[μ],factor,UdSfdUμ[μ])
        #println(" p[μ] = ", p[μ][1,1,1,1,1])
    end
end

function test2()
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 0
    Dim = 4
    NC = 3
    nowing = true

    U = Initialize_4DGaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = 2.1/2
    push!(gauge_action,β,plaqloop)

    show(gauge_action)

    x = Initialize_pseudofermion_fields(U[1],"staggered",nowing = nowing)
    gauss_distribution_fermion!(x)
    params = Dict()
    params["Dirac_operator"] = "staggered"
    params["mass"] = 0.5
    params["eps_CG"] = 1.0e-19
    params["verbose_level"] = 2
    D = Dirac_operator(U,x,params)

    Nf = 4

    println("Nf = $Nf")
    parameters_action = Dict()
    parameters_action["Nf"] = Nf
    fermi_action = FermiAction(D,parameters_action)

    y = similar(x)

    MDtest!(gauge_action,U,Dim,fermi_action,x,y)

end

test2()


text:
# Test tetrahedron delta function integration
using TetrahedronIntegration
using StaticArrays
using Test
using TetrahedronIntegration: _evaluate_piecewise_polynomial

@testset "delta" begin
    e0 = 0.4
    v0 = [0.0, 0.0, 1.0]
    L = 0.3
    de = L * v0[3] / 2
    @test delta_parallelepiped(e0, e0, v0, L) ≈ 1.0 / L
    @test delta_parallelepiped(e0 + de - 1e-3, e0, v0, L) ≈ 1.0 / L
    @test delta_parallelepiped(e0 - de + 1e-3, e0, v0, L) ≈ 1.0 / L
    @test delta_parallelepiped(e0 + de + 1e-3, e0, v0, L) ≈ 0
    @test delta_parallelepiped(e0 - de - 1e-3, e0, v0, L) ≈ 0

    # Test small v0
    @test delta_parallelepiped(e0, e0, [0.0, 0.0, 0.0], L) ≈ 0
end

@testset "delta_polynomial" begin
    e1234 = SVector(1.5, -0.2, 0.7, 0.3)
    polys = delta_tetrahedron_polynomial(e1234)
    for e in -0.3:0.1:1.6
        @test delta_tetrahedron(e, e1234) ≈ _evaluate_piecewise_polynomial(e, polys) atol=10*eps(eltype(e))
    end

    e1234 = SVector(1.0, -0.5, 1.0, -0.5)
    polys = delta_tetrahedron_polynomial(e1234)
    for e in -1.0:0.5:1.0
        @test delta_tetrahedron(e, e1234) ≈ _evaluate_piecewise_polynomial(e, polys)
    end
end

max_stars_repo_name: jaemolihm/TetrahedronIntegration

max_issues_repo_name: jaemolihm/TetrahedronIntegration

max_forks_repo_name: jaemolihm/TetrahedronIntegration


text:
function Simple_AON_Louvain(H::hypergraph;startclusters="singletons",gamma=1.0,
	clusterpenalty=0,kmax=maximum(keys(H.E)),maxits::Int64=100,
	randflag::Bool = false,verbose=true,Z0 = collect(1:length(H.D)))
	"""
    Simplest version of AON Louvain, the doesn't require setting an initial
    intensity function Ω. Intensity function is generated from one of a few
	different choices:
		singletons: learn omega from all nodes being in a singleton cluster
		cliquelouvain: use clique expansion + graph louvain warm start with parameter gamma
		starlouvain: use star expansion + graph louvain warm start with parameter gamma
	"""

	n = length(H.D)
	if startclusters == "singletons"
		Z0 = collect(1:n)
	elseif startclusters == "cliquelouvain"
		Z0 = CliqueExpansionModularity(H,gamma)
	elseif startclusters == "starlouvain"
		Z0 = StarExpansionModularity(H,gamma)
	else
		error("There is no starting clustering called $startclusters. Choose from: [\"singletons\", \"cliquelouvain\",\"starlouvain\"]")
	end

	d = Float64.(H.D)
	n = length(d)
	kmax = maximum(keys(H.E))
	He2n, weights = hypergraph2incidence(H)
    e2n = incidence2elist(He2n);
    n2e = incidence2elist(SparseArrays.sparse(He2n'))
    m = length(e2n)
    elen = zeros(Int64,m)
    for e = 1:m
        elen[e] = length(e2n[e])
    end
	β, γ,omega = learn_omega_aon(e2n,weights,Z0,kmax,d,n)
	Zset = AON_Louvain(n2e,e2n,weights,d,elen,β,γ,kmax,randflag,maxits,verbose,Z0,clusterpenalty)
	Z = Zset[:,end]
	return Z
end

function AON_Louvain(H::hypergraph,β::Vector{Float64},γ::Vector{Float64};maxits::Int64=100,
	verbose::Bool=true,randflag::Bool = false,clusterpenalty=0.0,Z0 = collect(1:length(H.D)))
	"""
	Calling AON_Louvain without explicit intensity function Ω. Just requires
	cut and volume penalties β and γ. Assumes hyperedges are of size >= 2.
	"""
	if length(H.E[1]) > 0
		println("This code assumes a hypergraph where hyperedges have at least two nodes.")
	end
	if γ[1] != 0 || β[1] != 0
		println("Code assumes there are no hyperedges of size 1. Setting β[1] = γ[1] = 0.")
		β[1] = 0
		γ[1] = 0
	end
	d = Float64.(H.D)
	n = length(d)
	kmax = maximum(keys(H.E))
	@assert(kmax == length(β))
	He2n, weights = hypergraph2incidence(H)
    e2n = incidence2elist(He2n);
    n2e = incidence2elist(SparseArrays.sparse(He2n'))
    m = length(e2n)
    elen = zeros(Int64,m)
    for e = 1:m
        elen[e] = length(e2n[e])
    end
	Zset = AON_Louvain(n2e,e2n,weights,d,elen,β,γ,kmax,randflag,maxits,verbose,Z0,clusterpenalty)
	Z = Zset[:,end];
end

function alternate_hypergraph_storage(H::hypergraph)
	"""
	It will be quicker in some cases to store the hypergraph in an alternate format.
	"""
	d = Float64.(H.D)
	n = length(d)
	kmax = maximum(keys(H.E))
	He2n, weights = HyperModularity.hypergraph2incidence(H)
	e2n = incidence2elist(He2n);
	n2e = incidence2elist(SparseArrays.sparse(He2n'))
	m = length(e2n)
	elen = zeros(Int64,m)
    for e = 1:m
        elen[e] = length(e2n[e])
    end
	return d,n,kmax,e2n,n2e,weights,elen
end

function AON_Louvain(H::hypergraph,Ω::IntensityFunction;α,clusterpenalty=0,kmax=maximum(keys(H.E)),maxits::Int64=100,bigInt::Bool=true,verbose=true,scan_order="random",Z0 = collect(1:length(H.D)))
	"""
	Call AON_Louvain with a function call similar to the symmetric HMLL function.
	"""
	randflag = !(scan_order == "lexical")
    β, γ, e2n, n2e,w,d,elen = AON_Inputs(H,Ω.ω,α,kmax)
    Zset = AON_Louvain(n2e,e2n,w,d,elen,β,γ,kmax,randflag,maxits,verbose,Z0,clusterpenalty)
    Z = Zset[:,end]
    return Z
end

function AON_Louvain(node2edges::Vector{Vector{Int64}},
    edge2nodes::Vector{Vector{Int64}},w::Vector{Float64},
    d::Vector{Float64},elen::Vector{Int64},
    β::Vector{Float64},γ::Vector{Float64},
    kmax::Int64,randflag::Bool=false,maxits::Int64=100,verbose::Bool=true,
    Zwarm::Vector{Int64}=Vector{Int64}(),clusterpenalty = 0)
    """
    Supernode version of Hypergraph Louvain for planted partition model
        (all-or-nothing cut penalties).

        * node2edges[i] = list of edge labels that node i is in
        * edge2nodes[e] = list of node labels that edge e contains
        * w[e] = weight of hyperedge e
        * d[i] = degree or weight of node i
        * elen[e] = number of nodes in hyperedge e
                     note that this may be > length(edge2nodes[e])
                     if some of the nodes are supernodes
        * β = scaling values of for hyperedge cuts
        * γ = scaling values for volumes
        * kmax = maximum hyperedge size
        * Zwarm = warm start clustering (optional)
        * randflag = whether or not to permute node order
        * maxits = Maximum number of greedy passes over the node set
        * clusterpenalty = include term clusterpenalty*log(k) in objective
            so if clusterpenalty > 0, there is incentive to form fewer clusters
    """

    n = length(d)
    # Step 1: greedy moves until no more improvement
    Z, improved = ANHL_Step(node2edges,edge2nodes,w,d,elen,β,γ,kmax,
                            randflag,maxits,verbose,Zwarm,clusterpenalty)

    # Store all clusterings found
    if improved
        Zs = Z
        Z_old = copy(Z)
    else
        Zs = Z
    end

    # As long as something is still improving each time you call step 1, keep going
    while improved

        # Step 2: Collapse the clustering into supernodes
        @assert(minimum(Z_old) > 0)
        e2n = deepcopy(edge2nodes)
        Change_id_H!(e2n,Z_old)
        # Merge hyperedges that are the same and remove size 1 hyperedges
        e2n, wSuper, elenSuper =  Condense_H(e2n,w,elen)
        He2n = elist2incidence(e2n,maximum(Z_old))
        n2e = incidence2elist(He2n,true)
        dSuper = condense_d(d,Z_old)
        ########
        # Step 1: Go back to greedy local moves, this time on the reduced hypergraph
        Zsuper, improved = ANHL_Step(n2e,e2n,wSuper,dSuper,elenSuper,β,γ,
                                kmax,randflag,maxits,verbose,Vector{Int64}(),clusterpenalty)
        N = length(Zsuper)    # N = number of supernodes = number clusters from last round

        @assert(minimum(Zsuper)>0)
        @assert(length(Zsuper)==length(dSuper))
        # println("$N, $improved")
        # Undo the procedure that merged nodes into super nodes, so that you get a clustering vector of length n
        # Do this only if the last call to Step 1 led to at least one greedy move.
        if improved
            # Extract what that new clustering is
            Z_new = zeros(Int64,n)

            # For each supernode, place all original node IDs that make it up
            # into a cluster of the supernodes label
            for i = 1:N

                # Get the cluster that supernode i is in
                SuperI_Cluster = Zsuper[i]

                # Get individual node ID that are in supernode i.
                SuperI_nodes = findall(x->x==i,Z_old)
                Z_new[SuperI_nodes] .= SuperI_Cluster
            end
            @assert(minimum(Z_new)>0)
            Zs = [Zs Z_new]
            Z_old = copy(Z_new)
        end
    end

    return Zs
end

function ANHL_Step(node2edges::Vector{Vector{Int64}},edge2nodes::Vector{Vector{Int64}},
    w::Vector{Float64},d::Vector{Float64},elen::Vector{Int64},
    β::Vector{Float64},γ::Vector{Float64},kmax::Int64,
    randflag::Bool=false,maxits::Int64=100,verbose::Bool=true,Zwarm::Vector{Int64}=Vector{Int64}(),clusterpenalty=0)
    """
    Basic step Louvain algorithm: iterate through nodes and greedily move
    nodes to adjacent clusters. Does not form supernodes and does not recurse.

    Should work also for degenerate hyperedges, though some bugs are possible

    Features:
    * has weight for size of hyperedge rather than computing the weights
    * allows you to warm start with a clustering
    * allows you to randomize node order
    """
    if verbose println("One step of all-or-nothing HyperLouvain") end
    m = length(edge2nodes)
    n = length(node2edges)

    # Randomize node visit order
    if randflag
        node_order = Random.randperm(n)
        # NOTE: may be faster to re-arrange
        #   all the data structures so that we march through
        #   Z in order.
    else
        node_order = 1:n
    end

    if verbose println("") end

    # Warm start clustering
    if length(Zwarm) > 0
        Z = renumber(Zwarm)
    else
        # Default is to put all nodes in their own cluster
        Z = collect(1:n)
    end

	K = maximum(Z)   			# number of clusters
    ClusVol = zeros(K)			# volume of clusters
	ClusSize = zeros(K)			# size of each cluster
	for i = 1:n
		ClusVol[Z[i]] += d[i]
		ClusSize[Z[i]] += 1
	end

    # Store cut penalties for each edge
    cutpenalty = zeros(m)
    for e = 1:m
        k = elen[e]      				# size of the edge
        cutpenalty[e] = β[k]*w[e] 	# penalty for cutting it
    end

    # Store node neighbors of each node
    Neighbs = NeighborList(node2edges, edge2nodes)

    # Store all the edges for the node, minus that node id
    # This makes computing cut changes faster
    edgelists = Vector{Vector{Vector{Int64}}}()
    for i = 1:n
        iedges = Vector{Vector{Int64}}()    # each node has a list of hyperedges it's in, i.e. a list of lists of nodes.
        for e in node2edges[i]
            edge = edge2nodes[e]
            push!(iedges, setdiff(edge,i))  # don't save i itself
        end
        push!(edgelists,iedges)
    end

    improving = true
    iter = 0
    changemade = false

    # Initialize vols and cuts to track efficient local changes
    r = kmax
    toler = 1e-8

    mainstart = time()
    while improving && iter < maxits

        iter += 1
        if mod(iter,1) == 0
            if verbose println("Louvain Iteration $iter") end
        end
        improving = false

        # visit each node in turn
        for i = node_order

            # Cluster index for node i
            Ci_ind = Z[i]

            # Get the size of the cluster
			clus_size = ClusSize[Ci_ind]

            # Get the indices of i's neighbors--these define clusters we might move to
            Ni = Neighbs[i]

            # Get the neighboring clusters of i
            NC = unique(Z[Ni])

            # The default is to not move the node
            BestC = Z[i]
            BestImprove = 0

            # Set of edges that i is in
            Cv = node2edges[i]

            # Set of hyperedges that node i is in, not including i itself
            Cv_list = edgelists[i]

            # Volume of the set currently
            vS = ClusVol[Ci_ind]
            dv = d[i]

            for j = 1:length(NC)

				# Check how much it would improve to move i to to cluster j
                Cj_ind = NC[j]

                if Cj_ind == Ci_ind
                    continue
                end

                # Change in volume
                vJ = ClusVol[Cj_ind]
                Δvol = 0
                for k = 1:kmax
					# Better if this is smaller
                    Δvol += γ[k]*((vS-dv)^k + (vJ+dv)^k - vS^k - vJ^k)
                end

                # Change in cut
                Δcut = 0
                for eid = 1:length(Cv)
                    # change in cut for edge e when v moves from
                    # cluster Ci to cluster Cj
                    e = Cv[eid]
                    edge_noi = Cv_list[eid]

                    if elen[e] > 1
						we = cutpenalty[e]
                        Δcut += move_cut(i,Z,edge_noi,Ci_ind,Cj_ind,we)
                    end
                end

                # Change in objective due to change in cluster number
                Δclus = 0
                if clusterpenalty > 0 && clus_size == 1
                    Δclus = clusterpenalty*(log(K-1)-log(K))
                end
                change = Δcut + Δvol + Δclus # want this to be negative

                # Check if this is currently the best possible greedy move to make
                # The tolerance helps with making some behavior more stable:
                #   you only move if there's a numerically nonzero reason to move.
                if change < BestImprove - toler
                    BestImprove = change
                    BestC = Cj_ind
                    improving = true
                end
            end

            # Move i to the best new cluster, only if it strictly improves modularity
            if BestImprove < -toler

                # update clustering
                ci_old = Z[i]
                Z[i] = BestC

                # Remove i from its old cluster...
                ClusVol[ci_old] -= dv
				ClusSize[ci_old] -= 1

                # ...and add it to its new cluster
                ClusVol[BestC] += dv
				ClusSize[BestC] += 1
                changemade = true
                improving = true # we have a reason to keep iterating!

                if clus_size == 1
                    # we moved a singleton to another cluster
                    K -= 1
                end

            end
        end
    end
    mainloop = time()-mainstart
    if ~changemade
        improved = false
        if verbose println("No nodes moved clusters") end
    else
        improved = true
    end
    if verbose println("Main loop took $mainloop seconds") end
	Z = renumber(Z)
    return Z, improved

end

"""
Compute the change from moving node v from Ci to Cj
"""
function move_cut(v::Int64,Z::Vector{Int64},erest::Vector{Int64},Ci_ind::Int64,
    Cj_ind::Int64,we::Float64)

    # p = Z[erest]            # set of clusters in e, not counting e
    p1 = Z[erest[1]]
    na = notsame(erest,Z,p1)

    if na
        # if the edge is cut regardless of whether v is moved, return 0
        return 0
    else
        # p1 = p[1]   # this the the cluster all other nodes are in
        if Ci_ind == p1
            # then moving v cuts the edge
            return we
        elseif Cj_ind == p1
            # then moving v to Cj_ind "uncuts" the edge
            return -we
        else
            # moving doesn't matter
            return 0
        end
    end
end

"""
e is the set of nodes in the hyperedge. Z stores their clusters
"""
function notsame(e::Vector{Int64},Z::Vector{Int64},p1::Int64)
    for i = 2:length(e)
        if Z[e[i]] != p1
            return true
        end
    end
    return false
end

function Condense_H(edge2nodes::Vector{Vector{Int64}},w::Vector{Float64})
    """
    Condense Hypergraph.

        A hypergraph stored as an edgelist may contain some edges
        that are identical. This function finds them and condenses
        them instead into a weighted hypergraph edgelist
    """
    Hdict = Dict{Vector{Int64},Float64}()   # maps from set of nodes to a weight

    for e = 1:length(edge2nodes)
        edge = edge2nodes[e]
        sort!(edge)
        we = w[e]
        wt = get(Hdict,edge,0)
        Hdict[edge] = wt+we
    end

    newlist = Vector{Vector{Int64}}()
    for edge in keys(Hdict)
        push!(newlist,edge)
        push!(wnew,Hdict[edge])
    end

    return newlist, wnew
end

function Condense_H(edge2nodes::Vector{Vector{Int64}},w::Vector{Float64},label::Vector{Int64})
    """
    Condense Hypergraph.

        A hypergraph stored as an edgelist may contain some edges
        that are identical. This function finds them and condenses
        them instead into a weighted hypergraph edgelist.

        In this version, each hyperedge is additionally associated with a
        label, and we only merge nodes if they have the same node set
        AND the same label.

        The edge label is mainly meant to store original hyperedge size.
        Since some of the nodes may be supernodes, we need to keep track of
        the original hyperedge size.

        Also, get rid of size 1 hyperedges in the process
    """
    Hdict = Dict{Vector{Int64},Float64}()   # maps from set of nodes to a weight

    @assert(length(label)==length(edge2nodes))

    for e = 1:length(edge2nodes)
        edge = edge2nodes[e]
        lab = label[e]
        sort!(edge)
        ky = [lab; edge]    # the key is the label + edge
        we = w[e]
        wt = get(Hdict,ky,0)
        Hdict[ky] = wt+we
    end

    newlist = Vector{Vector{Int64}}()
    wnew = Vector{Float64}()
    labelsnew = Vector{Int64}() # indices of edges you keep
    for ky in keys(Hdict)
        edge = ky[2:end]
        lab = ky[1]
        if length(edge) > 1
            push!(newlist,edge)
            push!(wnew,Hdict[ky])
            push!(labelsnew,lab)
        end
    end

    return newlist, wnew, labelsnew
end

function Change_id_H!(edge2nodes::Vector{Vector{Int64}},Z::Vector{Int64},uniqueflag::Bool=true)
"""
Given a clustering vector Z, change node IDs in the hypergraph
so that node i becomes node Z[i]. This will potentially lead
to many degenerate hyperedges where the same node shows up multiple times.

This is useful as a super-node step in clustering heuristics like louvain.

If uniqueflag == true, don't store degenerate hyperedges. Instead, just save
one copy of the vector
"""

    for e = 1:length(edge2nodes)
        for i = 1:length(edge2nodes[e])
            edge2nodes[e][i] = Z[edge2nodes[e][i]]
        end
        if uniqueflag
            unique!(edge2nodes[e])
        end
    end

    return edge2nodes
end

function condense_d(d::Vector{Float64},Z::Vector{Int64})
    """
    Condense a degree vector to a weights vector by summing up volumes of
    clusters.
    """

    n = maximum(Z)
    dnew = zeros(n)
    for i = 1:n
        S = findall(x->x==i,Z)
        dnew[i] = sum(d[S])
    end
    return dnew
end

function Elist_to_Hypergraph(elist::Vector{Vector{Int64}}, maxsize::Int64=25)
    """
    Converts a binary hypergraph incidence matrix to type "hypergraph".
    Assumes edges are properly sorted already
    """
    E = Dict{Integer, Dict}()
    for edge in elist
        # sort!(edge)
        if length(edge) > maxsize; continue; end
        sz = length(edge)
        if !haskey(E, sz)
            E[sz] = Dict{}()
        end
        E[sz][edge] = 1
    end

    D = zeros(Int64, n)
    for (sz, edges) in E
        for (e, _) in edges
            D[e] .+= 1
        end
    end
    N = 1:n
    return hypergraph(N, E, D)
end

function Hmat_to_Hypergraph(H::SparseArrays.SparseMatrixCSC, maxsize::Int64=25)
    """
    Converts a binary hypergraph incidence matrix to type "hypergraph"
    """
    E = Dict{Integer, Dict}()
    elist = incidence2elist(H)
    for edge in elist
        sort!(edge)
        if length(edge) > maxsize; continue; end
        sz = length(edge)
        if !haskey(E, sz)
            E[sz] = Dict{}()
        end
        E[sz][edge] = 1
    end

    D = zeros(Int64, n)
    for (sz, edges) in E
        for (e, _) in edges
            D[e] .+= 1
        end
    end
    N = 1:n
    return hypergraph(N, E, D)
end

function AON_Inputs(H,ω,α,kmax)
    """
    Given a hyperedge H, ω function, and α which parameterizes ω,
    return all of the data structures you need to run the
    all-or-nothing louvain algorithm
    """
    β = zeros(kmax)
    γ = zeros(kmax)
    for k = 2:kmax
        β[k] = log(ω([1,k],α))-log(ω([0,k],α))
        γ[k] = ω([1,k],α)-ω([0,k],α)
    end

    He2n, edge_weights = hypergraph2incidence(H)
    e2n = incidence2elist(He2n);
    n2e = incidence2elist(SparseArrays.sparse(He2n'))
    m = length(e2n)
    edge_len = zeros(Int64,m)
    for e = 1:m
        edge_len[e] = length(e2n[e])
    end
    # deg = vec(sum(He2n,dims = 1)) # only works if hypergraph is unweighted
	deg = 1.0*H.D # coercion to float for compatibility with AON_louvain
    # @show H.D-deg
    # @assert(deg == H.D)

    return β, γ, e2n, n2e, edge_weights,deg,edge_len

end

function learn_omega_aon(e2n,weights,Z,kmax,d,n)
    """"
    Given a fixed clustering, this computes the maximum likelihood estimate
    for the hyperedge intensity parameters omega specifically for the
    all-or-nothing modularity objective.

	These esimates are obtained by directly deriving the AON modularity
	with respect to omega_{1,k} and omega_{0,k}, probabilities for turning size-k
	tuples into hypergraphs, if they are completely contained inside a cluster
	or not.

    Assumes kmin is 2.
    e2n is the edge-to-node list
    d is the degree vector
    n is number of nodes
    Z is the clustering
    """

    L = maximum(Z)
    m = length(e2n)
    ClusVol = zeros(L)			# volume of clusters
    for i = 1:n
        ClusVol[Z[i]] += d[i]
    end

	# top row is weight of edges; bottom row is weight of cut edges
    EdgesAndCuts = zeros(2,kmax)

	# initialize to small positive number as a heuristic, to avoid taking log of 0
    EdgesAndCuts[1,:] = 0.01*ones(1,kmax)

    # Compute the cut penalty for each hyperedge size
    for j = 1:m
        edge = e2n[j]
		we = weights[j]
        k = length(edge)
        EdgesAndCuts[1,k] += we     # just counting number of hyperedges
        iscut = notsame(edge,Z,Z[edge[1]])
        if iscut
            EdgesAndCuts[2,k] += we
        end
    end

    # top row is in probability, bottom is out probability
    omega = zeros(2,kmax)
    volV = sum(ClusVol)

    for k = 2:kmax
        # get the volume piece
        volsums = 0
        for l = 1:L
            volsums += ClusVol[l]^k
        end
        volVk = volV^k

        omega[1,k] = (EdgesAndCuts[1,k] - EdgesAndCuts[2,k])/volsums
        omega[2,k] = EdgesAndCuts[2,k]/(volVk - volsums)
    end

    β = zeros(kmax)
    γ = zeros(kmax)
    for k = 2:kmax
        β[k] = log(omega[1,k])-log(omega[2,k])
		# if β[k] == Inf
		# 	# In cases where β[k] = Inf,
		# 	# we replace with a high but finite value instead
		# 	β[k] = 1e3
		# end
        γ[k] = omega[1,k]-omega[2,k]
    end
    return β, γ,omega
end

function learn_omega_aon(H::hypergraph,Z::Vector{Int64})
	"""
	Same function as above, but takes input type Hypergraph.
	"""
	d = H.D
	n = length(d)
	kmax = maximum(keys(H.E))
	e2n,weights = hyperedge_formatting(H)
	return learn_omega_aon(e2n,weights,Z,kmax,d,n)
end

function Kaminski_default(H::hypergraph)
	kmax = maximum(keys(H.E))
	EdgeList,weights = hyperedge_formatting(H)
	m = length(EdgeList)
	Mvec = zeros(kmax)
	for i = 1:m
		edge = EdgeList[i]
		we = weights[i]
		k = length(edge)
		Mvec[k] += we
	end
	volH = sum(H.D)
	beta = ones(kmax)
	beta[1] = 0
	gamma = zeros(kmax)
	for k = 2:kmax
		gamma[k] = Mvec[k]/(volH)^k
	end
	return beta, gamma
end

function omega_to_betagamma(omega)
	"""
	Go from AON intensity parameters to parameters omega and gamma.
	Page 7 in text.
	"""
	kmax = size(omega,2)
	β = zeros(kmax)
	γ = zeros(kmax)
	for k = 1:kmax
		β[k] = log(omega[1,k]/omega[2,k])
		γ[k] = (omega[1,k]-omega[2,k])/β[k]
	end
	return β, γ
end

function betagamma_to_omega(β, γ)
	"""
	Go from parameters omega and gamma to AON intensity parameters.
	Page 7 in text.
	"""
	kmax = length(β)
	omega = zeros(2,kmax)
	for k = 1:kmax
		omega[2,k] = β[k]*γ[k]/(exp(β[k])-1)
		omega[1,k,] = exp(β[k])*omega[2,k]
	end

	return omega
end

function allsame(v::Vector{Float64})
	"""
	Checks whether all entries of a vector are the same.
	"""
    v1 = v[1]
    for j = 2:length(v)
        if v[j] != v1
            return false
        end
    end
    return true
end

function allsame(v::Vector{Int64})
	"""
	Checks whether all entries of a vector are the same.
	"""
    v1 = v[1]
    for j = 2:length(v)
        if v[j] != v1
            return false
        end
    end
    return true
end

function modularity_aon(H::hypergraph,Z::Vector{Int64},omega::Array{Float64,2};likelihood=false)
	"""
	Computes the all or nothing modularity function.
		Equation (4), function Q in manuscript.
	"""
	β, γ = omega_to_betagamma(omega)
	n = length(Z)
	kmax = maximum(keys(H.E))
	kmin = maximum(keys(H.E))
	if kmin == 1
		error("Function requires hyperedges to have at least size 2.")
	end

	# m_k = Mvec[k] = weighted number of hyperedges of size k
	mvec = zeros(kmax)
	for k = kmin:kmax
        El = H.E[k]
        for edge in keys(El)
            mvec[k] = El[edge]
        end
    end

	L = maximum(Z)
    ClusVol = zeros(L)			# volume of clusters
    for i = 1:n
        ClusVol[Z[i]] += H.D[i]
    end

	obj = 0
	for k = kmin:kmax
		# Compute cut penalty of the objective
		El = H.E[k]
		for edge in keys(El)
			if ~allsame(Z[edge])
				obj -= El[edge]*β[k]
			end
		end

		# Volume penalty of the objective
		for l = 1:L
			obj -= β[k]*γ[k]*ClusVol[l]^k
		end
	end

	# Additional term that is constant with respect to clustervec
	for k = kmin:kmax
		obj += β[k]*mvec[k]
		El = H.E[k]
		for edge in keys(El)
			obj += El[edge]*log(omega[2,k])
			p = partitionize(Z[edge])
			bR = Combinatorics.multinomial(p...)
			obj -= prod(H.D[edge])*omega[2,k]*bR
		end
	end

	if likelihood == true
		# extra part for log-likelihood computation
		loglikhood = obj
		for k = kmin:kmax
			El = H.E[k]
			for edge in keys(El)
				aR = El[edge]
				loglikhood += aR*log(prod(H.D[edge]))
				p = partitionize(Z[edge])
				bR = Combinatorics.multinomial(p...)
				loglikhood -=  aR*log(bR) - log(factorial(aR))
			end
		end
		return obj, loglikhood
	else
		return obj
	end
end

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg


text:
immutable RidgeReg <: RegERM
    X::Matrix               # n x m matrix of n m-dimensional training examples
    y::Vector               # 1 x n vector with training classes
    n::Int                  # number of training examples
    m::Int                  # number of features
    kernel::Symbol          # kernel function
    regression_type::Symbol # ordinal, binomial, multinomial
    params::Hyperparameters # hyperparameters (e.g. λ)
end

function RidgeReg(X::Matrix, y::Vector; kernel::Symbol=:linear, λ::Float64=0.1)
    check_arguments(X, y, :ordinal)
    RidgeReg(X, y, size(X)..., kernel, :ordinal, RegularizationParameters(λ))
end

methodname(::RidgeReg) = "Linear Regression"
loss(::RidgeReg) = SquaredLoss()
regularizer(RidgeReg::RidgeReg, w::Vector) = L2reg(w, RidgeReg.params.λ)

# closed-form solution
function optimize(RidgeReg::RidgeReg; optimizer::Symbol=:closed_form)
    check_hyperparameters(RidgeReg.params)

    if optimizer == :closed_form
        y, X = RidgeReg.y, RidgeReg.X
        model = Model(X, y, RidgeReg.regression_type, RidgeReg.kernel)
        model.theta = (X'*X + eye(RidgeReg.m)/RidgeReg.params.λ)\X'*y
        model
    else
        invoke(optimize, (RegERM, Symbol), RidgeReg, optimizer)
    end
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_Hit" begin
     plane = Plane()
     ray1 = Ray(Point(0.0, 0.0, 2), -VEC_Z)

     intersection1 = ray_intersection(plane, ray1)
     @test typeof(intersection1) == HitRecord
     @test HitRecord(
          Point(0.0, 0.0, 0.0),
          Normal(0.0, 0.0, 1.0),
          Vec2d(0.0, 0.0),
          2.0,
          ray1
          ) ≈ intersection1

     ray2 = Ray( Point(3.0, 2.0, -2), VEC_Z )
     intersection2 = ray_intersection(plane, ray2)
     @test typeof(intersection2) == HitRecord
     @test HitRecord(
          Point(3.0, 2.0, 0.0),
          Normal(0.0, 0.0, -1.0),
          Vec2d(0., 0.),
          2.0,
          ray2
          ) ≈ intersection2

     P = Point(0., 5., 2.)
     Q = Point(6., 5., 0)
     ray3 = Ray( P, Q-P )
     intersection3 = ray_intersection(plane, ray3)
     @test typeof(intersection3) == HitRecord
     @test HitRecord(
          Q,
          Normal(0.0, 0.0, 1.0),
          Vec2d(0., 0.),
          1.0,
          ray3
          ) ≈ intersection3
end

@testset "test_noHit" begin
     plane = Plane()

     ray1 = Ray(Point(0, 0, 1), VEC_X)
     intersection1 = ray_intersection(plane, ray1)
     @test isnothing(intersection1)

     ray2 = Ray(Point(3, 2, -1), VEC_X+VEC_Y)
     intersection2 = ray_intersection(plane, ray2)
     @test isnothing(intersection2)

     ray3 = Ray(Point(0, 0, -1), -VEC_Z-VEC_Y)
     intersection3 = ray_intersection(plane, ray3)
     @test isnothing(intersection3)
end

@testset "test_Transformation" begin
     plane = Plane(translation(Vec(0.0, 0.0, 3.0)))

     ray1 = Ray(Point(0, 0, 5), -VEC_Z)
     intersection1 = ray_intersection(plane, ray1)
     @test typeof(intersection1) == HitRecord
     @test HitRecord(
          Point(0.0, 0.0, 3.0),
          Normal(0.0, 0.0, 1.0),
          Vec2d(0.0, 0.0),
          2.0,
          ray1
          ) ≈ intersection1

     ray2 = Ray(Point(3.0, 2.0, 0.), VEC_Z)
     intersection2 = ray_intersection(plane, ray2)
     @test typeof(intersection2) == HitRecord
     @test HitRecord(
          Point(3.0, 2.0, 3.0),
          Normal(0.0, 0.0, -1.0),
          Vec2d(0.0, 0.0),
          3.0,
          ray2
          ) ≈ intersection2

     plane3 = Plane(translation(Vec(0., 1., 0.))*rotation_x(π/4))
     ray3 = Ray(Point(0.0, 0.0, √2), Vec(0., 1., 0.))
     intersection3 = ray_intersection(plane3, ray3)
     @test HitRecord(
          Point(0., 1+√2, √2),
          Normal(0.0, -1.0, 1.0),
          Vec2d(0.0, 0.0),
          1+√2,
          ray3
          ) ≈ intersection3

     # Check if the plane failed to move by trying to hit the untransformed shape
     @test isnothing( ray_intersection(plane, Ray( Point(0, 0, 2), -VEC_Z ) ) )
end

@testset "test_Normals" begin
     plane1 = Plane(rotation_y(-π/4))
     P = Point(0, 0, 1)
     Q = Point(1, 0, 0)
     ray1 = Ray(P, Q-P)
     intersection1 = ray_intersection(plane1, ray1)

     @test intersection1.normal ≈ Normal(-1.0, 0.0, 1.0)

     plane2 = Plane(rotation_y(π/2))
     ray2 = Ray(Point(-1.0, 0.0, 0.0), Vec(1., 0., 0.))
     intersection2 = ray_intersection(plane2, ray2)
     @test intersection2.normal ≈ Normal(-1.0, 0.0, 0.0)
end

@testset "test_Normal_direction" begin
     # Scaling a plane by -1 keeps the plane the same but reverses its
     # reference frame
     plane = Plane(scaling(Vec(-1.0, -1.0, -1.0)))

     ray = Ray(Point(0.0, 0.0, 2.0), -VEC_Z)
     intersection = ray_intersection(plane, ray)

     @test intersection.normal ≈ Normal(0.0, 0.0, 1.0)
end

@testset "test_UV_Coordinates" begin
     plane = Plane()

     ray1 = Ray(Point(0.0, 0.0, 1.0), -VEC_Z)
     @test ray_intersection(plane, ray1).surface_point ≈ Vec2d(0.0, 0.0)

     ray2 = Ray(Point(0.5, 0.25, 1.0), -VEC_Z)
     @test ray_intersection(plane, ray2).surface_point ≈ Vec2d(0.5, 0.25)

     ray3 = Ray(Point(-2.25, 1.6, -1.0), VEC_Z)
     @test ray_intersection(plane, ray3).surface_point ≈ Vec2d(0.75, 0.6)
end


text:
import CommutativeRings: pseudo_ideal

@testset "basics" begin
    R = ZZ{Int}
    RX = R[:x]
    RYZ = R[:y,:z]
    @test iszero(Ideal(0))
    @test isone(Ideal(-1))
    @test Ideal(15, 21) == Ideal(3)
    @test Ideal([3, 5]) == R
    @test R == Ideal(-1)
    @test Ideal(3, ZZ(15)) == Ideal(3)

    x, = generators(RX)
    y, z = generators(RYZ)
    @test Ideal(x, x^2-x) == Ideal(x)
    @test Ideal(y + z, y^2 - z^2) == Ideal(y + z)

    @test pseudo_ideal(R, 15) == ZZ(15)
    @test pseudo_ideal(R, ZZ(Int8(127))) == 127
    @test pseudo_ideal(R, [15, 21]) == Ideal(3)
    i = Ideal(ZZ(3))
    @test pseudo_ideal(R, i) === i
    j = Ideal(ZZ(Int16(3)))
    @test i == j
end

@testset "operations" begin
    P = QQ{Int}[:x,:y,:z]
    x, y, z = generators(P)
    @test Ideal(x, y) isa Ideal

    A = [x*z^2 - 3y*z; x^2-2y; x*y - 5z]
    B = [x*y]
    ida = Ideal(A)
    idb = Ideal(B)
    idint = intersect(ida, idb)
    idsum = ida + idb

    @test ida isa Ideal{P}
    @test Ideal(ida.base) == ida
    @test zero(Ideal{P}) == zero(ida)
    @test zero(Ideal{P}) == Ideal(P(0))
    @test one(Ideal{P}) == Ideal(P(1))
    @test isone(Ideal(QQ(2)))
    @test issubset(idint, ida)
    @test issubset(idint, idb)
    @test issubset(ida, idsum)
    @test issubset(idb, idsum)
    @test A[1] * y + A[2] in ida
    @test 2x*y*z^2 - 9x*y^2 in idint
    @test B[1] * x + A[1] * y in idsum

    i0 = zero(ida)
    i1 = one(ida)
    @test i0 ⊆ ida ⊆ i1
    @test i0 + ida == ida
    @test i1 + ida == i1
    @test i1 ∩ ida == ida
    @test i0 ∩ ida == i0
    @test zero(ida) != idint^2 ⊆ ida * idb ⊆ idint
    @test one(ida) != idsum ⊆ one(ida)

    @test isone(idb^0)
    @test idb^1 === idb
    @test idb^2 == idb * idb
    @test idb^3 == idb * idb * idb
    @test idb^4 == idb * idb * idb * idb

    @test idb + x == x + idb
    @test idb * x == x * idb
    @test rem((x + y)^2, idb) == x^2 + y^2
end


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

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

# ╔═╡ 8d4298e0-6c44-11eb-3406-31a39340d1aa
begin
	import Pkg; Pkg.add("Plotly")
	Pkg.add("PyPlot")
	Pkg.add("PlutoUI")
	Pkg.add("SpecialFunctions")
	Pkg.add("StatsBase")
	Pkg.add("Plots")
	Pkg.add("LinearAlgebra")
	Pkg.add("SparseArrays")
	using PlutoUI
	using Plots
	using LinearAlgebra
	using SparseArrays
	using SpecialFunctions
	using StatsBase
end

# ╔═╡ 9637fe40-6c44-11eb-1df3-eb0b0d36db34
md"""

# **The golden Nautilus wheel** - _Distributions and the law of large numbers_

**Captain Bayes and Captain Venn** are investigating the **golden Nautilus wheel** and are wondering how one can make predictions about the minimum free days ratio.

A **rough rule of thumb** says that for most distributions the standard deviation can be used as a measure of how much probability mass is **spread around the mean value**, namely that within the following intervals one can find:

$P(X ∈ [μ - σ, μ + σ])  ≈ 67\%,$
$P(X ∈ [μ - 2σ, μ + 2σ]) ≈ 95\%,$
$P(X ∈ [μ - 3σ, μ + 3σ]) ≈ 99\%.$

In the following diagram you see the **probability distribution of wages** for the the four tasks of the golden Nautilus wheel, You can **change the wages** Captain Venn pays for each task

"""

# ╔═╡ 9de6ada0-6c56-11eb-1055-1b29b78e6112
md"""
**Use the slider to change the number of measurements sample size $n$ Captain Venn performs and watch the histogram**
"""

# ╔═╡ 8b8cd540-6d86-11eb-339e-7b8a2df774f9
md"Choose the **number of measurements / the sample size $n$** 👉 $(@bind n_sample_size Slider(100:100:10000, default=1000, show_value = true)) "

# ╔═╡ 84bfd150-6d8a-11eb-011e-052c39c8db9c
x_rand_var = [1, 3, 4, 6]

# ╔═╡ b7a7b220-6c56-11eb-1e6c-6f7595f02439
md"""
In order to see the variation of the sample we **repeat the measurements with fixed sample size** ($n$) $N$ times. This way we can **sample the distribution of the sample mean value** (which is a random variable in contrast to the real sample mean). Use the slider to change the number of repeated measurements $N$ of samples of size $n$**

You can change between the **distributions of wages and free days**.
"""

# ╔═╡ 835da8c0-6c57-11eb-0a71-e776862dbfce
md"""
Example of some other distributions: **Poisson distriubtion**
"""

# ╔═╡ 1b51ce3e-6c58-11eb-2539-a34e0576c612
@bind lambda Slider(0.1:0.1:19, default=1, show_value = true)

# ╔═╡ 0e35a290-6c58-11eb-2004-0b4b90f9754e
k =0:0.1:20

# ╔═╡ 14661b90-6c58-11eb-28bd-4399c5896666
poisson = lambda.^k./gamma.(k.+1).*exp(-lambda)

# ╔═╡ b016de70-6c59-11eb-2ae2-b9ee6e5f2abd
begin
	plot(
    k, poisson,
    line = (0.0 , 2.0 , :bar),
    normalize = false,
	bar_width = 0.13,
    marker = (6, 0.5, :none),
    markerstrokewidth = 5.,
    color = [:steelblue],
    fill = 0.9,
    orientation = :v,
    title = "The Poisson distribution with lambda: "*string(lambda),
	ylabel = "probability mass function",
	xlabel = "Points",
	label = :none,
	ylim =[0,1],)

	plot!(
    k, poisson,
    line = ( 1, 0., :path),
    normalize = false,
    marker = (3, 0.5, :o),
    markerstrokewidth = 3.,
    color = [:steelblue],
    fill = 0.,
    orientation = :v,
    title = "The Poisson distribution with lambda: "*string(lambda),
	ylabel = "probability mass function",
	xlabel = "k - number of events",
	label = :none,
	ylim =[0,1],
)

end

# ╔═╡ 06a6c020-6c5a-11eb-3516-e5d4f4740be9
sum(poisson)

# ╔═╡ 3d66c0a0-6c5b-11eb-0fd0-0bd81c9672c9
begin
	@userplot CirclePlot
@recipe function f(cp::CirclePlot)
    x, y, i = cp.args
    n = length(x)
    inds = circshift(1:n, 1 - i)
    linewidth --> range(0, 10, length = n)
    seriesalpha --> range(0, 1, length = n)
    aspect_ratio --> 1
    label --> false
    x[inds], y[inds]
end

n = 150
t = range(0, 2π, length = n)
x = sin.(t)
y = cos.(t)

anim = @animate for i ∈ 1:n
    circleplot(x, y, i)
end
gif(anim, "anim_fps15.gif", fps = 15)
end

# ╔═╡ 48f7c480-6d84-11eb-1f24-e16d7f649c9c
phi = (1+sqrt(5))/2

# ╔═╡ 534e1970-6d84-11eb-3f16-bb6aaa46db6f
prob = phi.^(0:3)./sum(phi.^(0:3))

# ╔═╡ 1df828e0-6d86-11eb-0d3c-2398e6700066
sample = 5 .- sum(rand(1,n_sample_size) .<= cumsum(prob), dims=1)

# ╔═╡ cfdaa470-6d86-11eb-03a6-ef281daeebb2
sample_count = [count(==(i), sample) for i = 1:4]

# ╔═╡ 5e0bfb10-6d8a-11eb-06ef-9335eb90064c
begin
	plot(x_rand_var, sample_count,
    line = (0.0 , 2.0 , :bar),
    normalize = false,
	bar_width = 0.04,
    marker = (6, 0.5, :none),
    markerstrokewidth = 5.,
    color = [:steelblue],
    fill = 0.9,
    orientation = :v,
    title = "The Poisson distribution with lambda: ",
	ylabel = "probability mass function",
	xlabel = "Points",
	label = :none,
	ylim =[0,n_sample_size/2])

	plot!(
    x_rand_var, sample_count,
    line = ( 1, 0., :path),
    normalize = false,
    marker = (5, 1.0, :o),
    markerstrokewidth = 3.,
    color = [:steelblue],
    fill = 0.,
    orientation = :v,
    title = "The Poisson distribution with lambda: "*string(lambda),
	ylabel = "probability mass function",
	xlabel = "k - number of events",
	label = :none,
	ylim =[0,n_sample_size/2])
end

# ╔═╡ Cell order:
# ╠═8d4298e0-6c44-11eb-3406-31a39340d1aa
# ╟─9637fe40-6c44-11eb-1df3-eb0b0d36db34
# ╟─9de6ada0-6c56-11eb-1055-1b29b78e6112
# ╟─8b8cd540-6d86-11eb-339e-7b8a2df774f9
# ╠═84bfd150-6d8a-11eb-011e-052c39c8db9c
# ╟─5e0bfb10-6d8a-11eb-06ef-9335eb90064c
# ╟─b7a7b220-6c56-11eb-1e6c-6f7595f02439
# ╟─835da8c0-6c57-11eb-0a71-e776862dbfce
# ╟─1b51ce3e-6c58-11eb-2539-a34e0576c612
# ╟─b016de70-6c59-11eb-2ae2-b9ee6e5f2abd
# ╠═0e35a290-6c58-11eb-2004-0b4b90f9754e
# ╠═14661b90-6c58-11eb-28bd-4399c5896666
# ╠═06a6c020-6c5a-11eb-3516-e5d4f4740be9
# ╠═3d66c0a0-6c5b-11eb-0fd0-0bd81c9672c9
# ╠═48f7c480-6d84-11eb-1f24-e16d7f649c9c
# ╠═534e1970-6d84-11eb-3f16-bb6aaa46db6f
# ╠═1df828e0-6d86-11eb-0d3c-2398e6700066
# ╠═cfdaa470-6d86-11eb-03a6-ef281daeebb2


text:
"""
$(TYPEDEF)

Base for all types which explicitly represent a surreal number.
"""
abstract type Surreal <: Number end

Surreal(x::Surreal) = x
Surreal(x) = SurrealForm(x)

"""
	left(x)::SurrealSet

Get the left subset of the canonical surreal number form of `x`
"""
function left end

"""
	right(x)::SurrealSet

Get the right subset of the canonical surreal number form of `x`.
"""
function right end

"""
	$(FUNCTIONNAME)(x)

Equivalent to `(left(x), right(x))`.
"""
leftright(x) = (left(x), right(x))

"""
	issurreal(x::Number)::Bool

Whether the number `x` is a surreal number.
"""
issurreal(::Surreal) = true
issurreal(::Real) = true
issurreal(x::Complex) = isreal(x)

"""
	isordinal(x::Number)::Bool

Whether the number `x` is an ordinal number.
"""
isordinal(x::Real) = isinteger(x) && x > 0

"""
	generation(x::Number)

Get the generation of the surreal number `x`.
"""
function generation end


text:
@testset "Verify evaluation of Physicist hermite functions" begin

    # @testset "Verify evaluation of Physicist hermite functions" begin
    x = randn(100)
    m = 5
    dV = vander(PhyHermite(m; scaled = false), 0, x)

    for i=0:5
        @test norm(dV[:,i+1] - FamilyPhyPolyHermite[i+1].(x) .* exp.(-x.^2/2))<1e-8
    end

    dV = vander(PhyHermite(m; scaled = true), 0, x)

    for i=0:5
        @test norm(dV[:,i+1] - FamilyScaledPhyPolyHermite[i+1].(x).*exp.(-x.^2/2))<1e-8
    end
end

@testset "Verify 1st derivative of Physicist hermite functions" begin
    # Unscaled
    x = randn(100)
    m = 5
    dV = vander(PhyHermite(m; scaled = false), 1, x)

    for i=0:5
        F = PhyHermite(i; scaled = false)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(F, xi), x) )<1e-8
        @test norm(dV[:,i+1] - (FamilyPhyPolyHermite[i+1].(x) .* (-x) .* exp.(-x.^2/2) +
            derivative(FamilyPhyPolyHermite[i+1],1).(x) .* exp.(-x.^2/2)))<1e-8
    end

    # Scaled
    x = randn(100)
    m = 5
    dV = vander(PhyHermite(m; scaled = true), 1, x)

    for i=0:5
        F = PhyHermite(i; scaled = true)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(F, xi), x) )<1e-8

        @test norm(dV[:,i+1] - 1/Cphy(i)*(FamilyPhyPolyHermite[i+1].(x) .* (-x) .* exp.(-x.^2/2) +
            derivative(FamilyPhyPolyHermite[i+1],1).(x) .* exp.(-x.^2/2)))<1e-8
    end

end

@testset "Verify second derivative of Physicist hermite functions" begin
    # Unscaled
    x = randn(100)
    m = 5
    dV = vander(PhyHermite(m; scaled = false), 2, x)

    for i=0:5
        Hn = FamilyPhyPolyHermite[i+1].(x)
        Hnp = derivative(FamilyPhyPolyHermite[i+1],1).(x)
        Hnpp = derivative(FamilyPhyPolyHermite[i+1],2).(x)
        E = exp.(-x.^2/2)
        F = PhyHermite(i; scaled = false)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(y->ForwardDiff.derivative(z->F(z), y), xi), x))<1e-8

        @test norm(dV[:,i+1] - (Hnpp .*E -2 .*x .* Hnp.*E + Hn .* (x.^2 .- 1.0) .* E))<1e-8
    end

    # Scaled

    x = randn(100)
    m = 5
    dV = vander(PhyHermite(m; scaled = true), 2, x)

    for i=0:5
        Hn = FamilyPhyPolyHermite[i+1].(x)
        Hnp = derivative(FamilyPhyPolyHermite[i+1],1).(x)
        Hnpp = derivative(FamilyPhyPolyHermite[i+1],2).(x)
        E = exp.(-x.^2/2)

        F = PhyHermite(i; scaled = true)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(y->ForwardDiff.derivative(z->F(z), y), xi), x))<1e-8

        @test norm(dV[:,i+1] - (1/Cphy(i))*(Hnpp .*E -2 .*x .* Hnp.*E + Hn .* (x.^2 .- 1.0) .* E))<1e-8
    end
end

@testset "Verify third derivative of Physicist hermite functions" begin
    # Unscaled
    x = randn(100)
    m = 5
    dV = vander(PhyHermite(m; scaled = false), 3, x)

    for i=0:5
        F = PhyHermite(i; scaled = false)
        @test norm(dV[:,i+1] - map(xi->Diff(F, 3)(xi),x))<1e-8
    end

    # Scaled

    x = randn(100)
    m = 5
    dV = vander(PhyHermite(m; scaled = true), 3, x)

    for i=0:5
        F = PhyHermite(i; scaled = true)
        @test norm(dV[:,i+1] - map(xi->Diff(F, 3)(xi),x))<1e-8
    end
end

# @testset "Verify integration of Physicist hermite functions" begin
#     Nx = 10
#     x = randn(10)
#
#     # Unscaled
#     dV = vander(PhyHermite(8; scaled = false), -1, x)
#
#     for i=0:8
#         for j=1:Nx
#             @test abs(dV[j,i+1] - quadgk(y-> PhyHermite(i)(y), 0.0, x[j], rtol=1e-12)[1])<1e-8
#         end
#     end
#
#     # Scaled
#     dV = vander(PhyHermite(8; scaled = true), -1, x)
#
#     for i=0:8
#         for j=1:Nx
#             @test abs(dV[j,i+1] - quadgk(y-> PhyHermite(i; scaled=true)(y), 0.0, x[j], rtol=1e-12)[1])<1e-8
#         end
#     end
#
# end


text:
## Exercise 4-12
## Read about spirals at https://en.wikipedia.org/wiki/Spiral; then write a program that draws an Archimedan spiral as in Archimedan spiral.
## https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#fig04-4
println("Ans : ")

println("End.")


text:
function AbstractPlotting.plot!(plot::Plot(AbstractVector{<: Complex}))
    plot[:axis, :labels] = ("Re(x)", "Im(x)")
    lines!(plot, lift(im-> Point2f0.(real.(im), imag.(im)), x[1]))
end

"""
    showlibrary(lib::Symbol)::Scene

Shows all colour gradients in the given library.
Returns a Scene with these colour gradients arranged
as horizontal colourbars.
"""
function showlibrary(lib::Symbol)::Scene

    cgrads = sort(PlotUtils.cgradients(lib))

    PlotUtils.clibrary(lib)

    showgradients(cgrads)

end

"""
    showgradients(
        cgrads::AbstractVector{Symbol};
        h = 0.0, offset = 0.2, textsize = 0.7,
        resolution = (800, length(cgrads) * 84)
    )::Scene

Plots the given colour gradients arranged as horizontal colourbars.
If you change the offsets or the font size, you may need to change the resolution.
"""
function showgradients(
        cgrads::AbstractVector{Symbol};
        h = 0.0,
        offset = 0.4,
        textsize = 0.7,
        resolution = (800, length(cgrads) * 84),
        monospace = true
    )::Scene

    scene = Scene(resolution = resolution)

    map(collect(cgrads)) do cmap

         c = to_colormap(cmap)

         cbar = image!(
             scene,
             range(0, stop = 10, length = length(c)),
             range(0, stop = 1, length = length(c)),
             reshape(c, (length(c),1)),
             show_axis = false
         )[end]

         cmapstr = monospace ? UnicodeFun.to_latex("\\mono{$cmap}") : string(cmap, ":")

         text!(
             scene,
             cmapstr,
             position = Point2f0(-0.1, 0.5 + h),
             align = (:right, :center),
             show_axis = false,
             textsize = textsize
         )

         translate!(cbar, 0, h, 0)

         h -= (1 + offset)

    end

    scene

end


text:
function importcscmat(
    m::Integer,
    n::Integer,
    colptr,
    rowindices,
    values::Vector{T};
    jumbled::Bool = false,
    desc::Descriptor = Descriptors.NULL
) where {T}
A = Ref{libgb.GrB_Matrix}()
m = libgb.GrB_Index(m)
n = libgb.GrB_Index(n)
colsize = libgb.GrB_Index(sizeof(colptr))
rowsize = libgb.GrB_Index(sizeof(rowindices))
valsize = libgb.GrB_Index(sizeof(values))

col = Ptr{libgb.GrB_Index}(Libc.malloc(colsize))
unsafe_copyto!(col, Ptr{UInt64}(pointer(colptr .- 1)), length(colptr))
row = Ptr{libgb.GrB_Index}(Libc.malloc(rowsize))
unsafe_copyto!(row, Ptr{UInt64}(pointer(rowindices .- 1)), length(rowindices))
val = Ptr{T}(Libc.malloc(valsize))
unsafe_copyto!(val, pointer(values), length(values))
libgb.GxB_Matrix_import_CSC(
    A,
    toGrB_Type(T),
    m,
    n,
    Ref{Ptr{libgb.GrB_Index}}(col),
    Ref{Ptr{libgb.GrB_Index}}(row),
    Ref{Ptr{Cvoid}}(val),
    colsize,
    rowsize,
    valsize,
    false,
    jumbled,
    desc
)
return GBMatrix{T}(A[])
end

function toGB(S::SparseMatrixCSC)
    return importcscmat(S.m, S.n, S.colptr, S.rowval, S.nzval)
end

function importcscvec(
    n::Integer, vi, vx::Vector{T};
    jumbled::Bool = false, desc::Descriptor = Descriptors.NULL
) where {T}
    v = Ref{libgb.GrB_Vector}
    n = libgb.GrB_Index(n)
    vi_size = libgb.GrB_Index(sizeof(vi))
    vx_size = libgb.GrB_Index(sizeof(vx))

    indices = Ptr{libgb.GrB_Index}(Libc.malloc(vi_size))
    unsafe_copyto!(indices, Ptr{UInt64}(pointer(vi .- 1)), length(vi))

    values = Ptr{libgb.GrB_Index}(Libc.malloc(vx_size))
    unsafe_copyto!(values, pointer(vx), length(vx))
    libgb.GxB_Vector_import_CSC(
        v,
        toGrB_Type(T),
        n,
        Ref{Ptr{libgb.GrB_Index}}(indices),
        Ref{Ptr{Cvoid}}(values),
        vi_size,
        vx_size,
        false,
        length(vx),
        false,
        desc
    )
end

function toGB(v::SparseVector)
    return importcscvec(v.n, v.nzind, v.nzval)
end

function importdensematrix(
    m::Integer, n::Integer, A::Matrix{T};
    desc::Descriptor = Descriptors.NULL
) where {T}
    C = Ref{libgb.GrB_Matrix}()
    m = libgb.GrB_Index(m)
    n = libgb.GrB_Index(n)
    Asize = libgb.GrB_Index(sizeof(A))

    Ax = Ptr{T}(Libc.malloc(Asize))
    unsafe_copyto!(Ax, pointer(A), length(A))
    libgb.GxB_Matrix_import_FullC(
        C,
        toGrB_Type(T),
        m,
        n,
        Ref{Ptr{Cvoid}}(Ax),
        Asize,
        false,
        desc
    )
    return GBMatrix{T}(C[])
end

function toGB(M::Matrix)
    return importdensematrix(size(M, 1), size(M, 2), M)
end

function importdensevec(
    n::Integer, v::Vector{T};
    desc::Descriptor = Descriptors.NULL
) where {T}
    w = Ref{libgb.GrB_Vector}()
    n = libgb.GrB_Index(n)
    vsize = libgb.GrB_Index(sizeof(v))

    vx = Ptr{T}(Libc.malloc(vsize))
    unsafe_copyto!(vx, pointer(v), length(v))
    libgb.GxB_Vector_import_Full(
        w,
        toGrB_Type(T),
        n,
        Ref{Ptr{Cvoid}}(vx),
        vsize,
        false,
        desc
    )
    return GBVector{T}(w[])
end

function toGB(v::Vector)
    return importdensevec(size(v)..., v)
end


text:
using NeuroIO
using Base.Test

include("NCS.jl")
include("NEV.jl")
include("PLX.jl")


text:
using Revise

using ANOVADDPTest
using ANOVADDPTest: dp_mass, n_clusters,  cluster_sizes
using ANOVADDPTest: active_clusters, passive_clusters, cluster_capacity, update!
using Distributions
using SpecialFunctions
using StatsBase: mode
using Statistics: mean, var
using Random
using Test
using CSV
using Tables

X = CSV.File("Xmat.csv") |> Tables.matrix;
y = CSV.File("yvec.csv") |> Tables.matrix;
y = y[:];
ch, data0 = anova_bnp_normal(y, X);
simple_effect_probabilities(ch, data0)


text:
###
# Hcat
####
hcat(A::AbstractQuasiArray...) = apply(hcat, A...)

ApplyStyle(::typeof(hcat), ::Type{<:AbstractQuasiArray}...) = LazyQuasiArrayApplyStyle()

axes(f::ApplyQuasiMatrix{<:Any,typeof(hcat)}) = (axes(f.args[1],1), Base.OneTo(sum(size.(f.args,2))))

function getindex(f::ApplyQuasiMatrix{T,typeof(hcat)}, k::Number, j::Number) where T
    ξ = j
    for A in f.args
        n = size(A,2)
        ξ ≤ n && return T(A[k,ξ])::T
        ξ -= n
    end
    throw(BoundsError(f, (k,j)))
end


text:
# Clang#master + julia 1.6
using Clang.Generators

# other options, see Clang.JLLEnvs.JLL_ENV_TRIPLES
args = get_default_args("x86_64-linux-gnu")

# ctx = create_context("affinity.c", args)
ctx = create_context("affinity.c", args)

build!(ctx)


text:
module DTW
include("step_pattern.jl")
using Distances: PreMetric

"""
Given a Rule-pattern, computes its cost
"""
@inline function step_pattern_cost(step_pattern::StepPattern, matrix::Matrix{T}, i::Int64,
    j::Int64)::Tuple{T, Int64} where {T<:Number}

    step_costs = Vector{T}(undef, length(step_pattern.rules) + 1)
    for (k, rule) in enumerate(step_pattern.rules)
        rule_cost::T = 0.0
        stop = false
        for comp::Component in rule.components
            if i - comp.row > 0 && j - comp.col > 0
                val = matrix[i - comp.row, j - comp.col]
                if isinf(val)
                    # a value outside the window
                    step_costs[k] = Inf
                    stop = true
                    break
                end
                rule_cost += comp.weight * val
            else
                # this rule would connect an element outside the matrix
                step_costs[k] = Inf
                stop = true
                break
            end
        end
        if !stop
            step_costs[k] = rule.weight * rule_cost + rule.bias
        end
    end
    # if no rule is allowed, use Inf
    step_costs[end] = Inf
    return findmin(step_costs)
end

"""
DTW functions working with Vectors of Vectors, that is: each point can
have a variable numer of features.

Returns a matrix containing at each entry the dtw cost to reach that entry.
The total dtw_cost can be found in [end, end].
"""
function dtw_matrix(a::T, b::T, dist_fn::PreMetric; window_fn::Function=(i, j) -> true,
                    step_pattern::StepPattern=symmetric2) where {T}

    matrix = fill(Inf, length(a), length(b))
    matrix[1, 1] = dist_fn(a[1], b[1])
    for j in 1:size(matrix, 2), i in 1:size(matrix, 1)
        # only if in window and not (1, 1)
        if window_fn(i, j) && i + j > 2
            matrix[i, j] = dist_fn(a[i], b[j])  # Calculate distance of current signal indices
            matrix[i, j], _ = step_pattern_cost(step_pattern, matrix, i, j)
        end
    end
    return matrix
end

function dtw_matrix(dist_matrix::Matrix{Float64}; window_fn::Function=(i, j) -> true,
                    step_pattern::StepPattern=symmetric2)

    matrix = fill(Inf::Float64, size(dist_matrix)...)
    matrix[1, 1] = dist_matrix[1, 1]
    for j in 1:size(matrix, 2), i in 1:size(matrix, 1)
        # only if in window and not (1, 1)
        if window_fn(i, j) && i + j > 2
            matrix[i, j] = dist_matrix[i, j]
            matrix[i, j], _ = step_pattern_cost(step_pattern, matrix, i, j)
        end
    end

    return matrix
end

"""
Modifies dist_matrix in place to save ram
"""
function dtw_matrix!(dist_matrix::Matrix{Float64}; window_fn::Function=(i, j) -> true,
                     step_pattern::StepPattern=symmetric2)

    for j in 1:size(dist_matrix, 2), i in 1:size(dist_matrix, 1)
        # only if in window and not (1, 1)
        if window_fn(i, j) && i + j > 2
            dist_matrix[i, j], _ = step_pattern_cost(step_pattern, dist_matrix, i, j)
        end
    end

    return dist_matrix
end

"""
DTW functions working with Vectors of Vectors, that is: each point can
have a variable numer of features.

Returns a matrix containing at each entry the rule used for the best path.
"""
function dtw_matrix_directions(a::T, b::T, dist_fn::PreMetric;
                               window_fn::Function=(i, j) -> true,
                               step_pattern::StepPattern=symmetric2)::Union{Matrix{Int64},
                                                                             Nothing} where {T}

    directions = zeros(Int64, length(a), length(b))
    matrix = fill(Inf, length(a), length(b))
    matrix[1, 1] = dist_fn(a[1], b[1])
    for j in 1:size(matrix, 2), i in 1:size(matrix, 1)
        # only if in window and not (1, 1)
        if window_fn(i, j) && i + j > 2
            matrix[i, j] = dist_fn(a[i], b[j])  # Calculate distance of current signal indices
            matrix[i, j], directions[i, j] = step_pattern_cost(step_pattern, matrix, i, j)
        end
    end
    if isinf(matrix[end, end])
        # no valid path
        return nothing
    end

    return directions
end

function dtw_matrix_directions(dist_matrix::Matrix{Float64};
                               window_fn::Function=(i, j) -> true,
                               step_pattern::StepPattern=symmetric2)::Union{Matrix{Int64},
                                                                             Nothing}

    directions = zeros(Int64, size(dist_matrix)...)
    matrix = fill(Inf, size(dist_matrix)...)
    matrix[1, 1] = dist_matrix[1, 1]
    for j in 1:size(matrix, 2), i in 1:size(matrix, 1)
        # only if in window and not (1, 1)
        if window_fn(i, j) && i + j > 2
            matrix[i, j] = dist_matrix[i, j]
            matrix[i, j], directions[i, j] = step_pattern_cost(step_pattern, matrix, i, j)
        end
    end

    if isinf(matrix[end, end])
        # no valid path
        return nothing
    end

    return directions
end

"""
Modifies dist_matrix in place to save ram
"""
function dtw_matrix_directions!(dist_matrix::Matrix{Float64};
                                window_fn::Function=(i, j) -> true,
                                step_pattern::StepPattern=symmetric2)::Union{Matrix{Int64},
                                                                              Nothing}

    directions = zeros(Int64, size(dist_matrix)...)
    for j in 1:size(dist_matrix, 2), i in 1:size(dist_matrix, 1)
        # only if in window and not (1, 1)
        if window_fn(i, j) && i + j > 2
            dist_matrix[i, j], directions[i, j] = step_pattern_cost(step_pattern,
                                                                    dist_matrix, i, j)
        end
    end
    if isinf(dist_matrix[end, end])
        # no valid path
        return nothing
    end

    return directions
end

"""
Performs the trackback algorithm and returns the best path found starting
from the first index of the original input sequences (i.e. the first row
and column are discarded).
"""
@views function trackback(matrix::Matrix{Int64},
                          step_pattern::StepPattern=symmetric2)::Matrix{Int64}
    pos = collect(size(matrix))
    at_start_index = false
    shortest_path::Vector{Vector{Int64}} = []
    push!(shortest_path, copy(pos))  # We start at the end and work back
    L = length(step_pattern.rules)

    while !at_start_index

        # with direction matrix
        dir = matrix[pos...]

        best_rule = step_pattern.rules[dir]

        # computing the movements according to this rule
        for comp in best_rule.components[(end - 1):-1:1]
            end_comp = pos - [comp.row, comp.col]
            pushfirst!(shortest_path, end_comp)
        end

        # update pos
        pos -= [best_rule.components[1].row, best_rule.components[1].col]

        # Stop once you're at the beginning
        if pos == [1, 1]
            at_start_index = true
        end
    end
    return vcat(transpose(shortest_path)...)
end

"""
Perform DTW and path trackback; returns `Nothing` if no path is found.
"""
function dtw_path(a::T, b::T, dist_fn::PreMetric; window_fn::Function=(i, j) -> true,
                  step_pattern::StepPattern=symmetric2)::Union{Matrix{Int64},
                                                                Nothing} where {T}
    matrix = dtw_matrix_directions(a, b, dist_fn; window_fn, step_pattern)
    if matrix === nothing
        return nothing
    end
    return trackback(matrix, step_pattern)
end
function dtw_path(dist_matrix::Matrix{Float64}; window_fn::Function=(i, j) -> true,
                  step_pattern::StepPattern=symmetric2)::Union{Matrix{Int64},
                                                                Nothing}
    matrix = dtw_matrix_directions(dist_matrix; window_fn, step_pattern)
    if matrix === nothing
        return nothing
    end
    return trackback(matrix, step_pattern)
end
function dtw_path!(dist_matrix::Matrix{Float64}; window_fn::Function=(i, j) -> true,
                   step_pattern::StepPattern=symmetric2)::Union{Matrix{Int64},
                                                                 Nothing}
    matrix = dtw_matrix_directions!(dist_matrix; window_fn, step_pattern)
    if matrix === nothing
        return nothing
    end
    return trackback(matrix, step_pattern)
end
end # module


text:
using Interpolations
"""
    getTVDSS(md::Float64, mds::Array{Float64}, tvds::Array{Float64})
Interpolates an array of MD, TVDSS and gives the tvdss at a desired depth"""
function getTVDSS(md::Float64, mds::Array{Float64},tvds::Array{Float64} )
    ipt = LinearInterpolation(mds, tvds)
    ko = ipt(md)
    return ko
end

"""
    getTemperature(md::Float64, mds::Array{Float64}, temps::Array{Float64})
Interpolates an array of MD, Temperature and gives the temperature at a desired depth
"""
function getTemperature(md::Float64, mds::Array{Float64}, temps::Array{Float64})
    ipt = LinearInterpolation(mds, temps)
    ko = ipt(md)
    return ko
end

max_stars_repo_name: leytzher/Shot

max_issues_repo_name: leytzher/Shot

max_forks_repo_name: leytzher/Shot


text:
print("plot: execute functions ...\n")

J = 8
n = 2^J
x = testfunction(n, "Bumps")
y = dwt(x, wavelet(WT.cdf97, WT.Lifting))
d,l = wplotdots(y, 0.1, n)
A = wplotim(y)

x = randn(64,64)
L = 2
xts = wplotim(x, L, wavelet(WT.db3, WT.Filter))


text:
struct StepCG end

struct CyclicIterator
    start::Int
    n::Int
    rev::Bool
end
CyclicIterator(start, n; rev=false) = CyclicIterator(start, n, rev)
function Base.iterate(iter::CyclicIterator, i = iter.start)
    if iter.rev
        i >= 1 && return i, i - 1
        i + iter.n > iter.start && return i + iter.n, i - 1
    else
        i <= iter.n && return i, i + 1
        i - iter.n < iter.start && return i - iter.n, i + 1
    end
    return nothing
end

@params struct ApproxInverseHessian{T} <: AbstractMatrix{T}
    S::AbstractMatrix{T}
    Y::AbstractMatrix{T}
    p::AbstractVector{T}
    prev_x::AbstractVector{T}
    prev_g::AbstractVector{T}
    oldest::Base.RefValue{Int}
    activen::Base.RefValue{Int}
end
function ApproxInverseHessian(nvars::Int, n::Int)
    S = zeros(nvars, n)
    Y = zeros(nvars, n)
    p = zeros(n)
    prev_x = zeros(nvars)
    prev_g = zeros(nvars)
    oldest = 1
    return ApproxInverseHessian(S, Y, p, prev_x, prev_g, Ref(oldest), Ref(0))
end
Base.size(invH::ApproxInverseHessian) = (size(invH.S, 1), size(invH.S, 1))
function update!(invH::ApproxInverseHessian, g, x)
    @unpack S, Y, p, prev_x, prev_g, oldest, activen = invH
    n = size(S, 2)
    ind = oldest[]
    new_s = x .- prev_x
    new_y = g .- prev_g
    new_p = 1 / dot(new_s, new_y)
    if isfinite(new_p)
        S[:,ind] .= new_s
        Y[:,ind] .= new_y
        p[ind] = new_p
        oldest[] = oldest[] == n ? 1 : oldest[] + 1
        activen[] = min(activen[] + 1, n)
    else
        activen[] = 0
        oldest[] = 1
    end
    prev_x .= x
    prev_g .= g
    return invH
end
function reset!(invH::ApproxInverseHessian)
    invH.activen[] = 0
    invH.oldest[] = 1
    return invH
end
function LinearAlgebra.mul!(c::AbstractVector{T}, invH::ApproxInverseHessian{T}, b::AbstractVector{T}) where {T}
    @unpack S, Y, p, prev_x, prev_g, activen = invH
    if debugging[]
        println("    activen = $(activen[])")
    end
    n = min(activen[], size(S, 2))
    c .= b
    n == 0 && return c
    oldest = invH.oldest[]
    newest = oldest == 1 ? n : oldest - 1
    for i in CyclicIterator(newest, n, rev=true)
        @views c .-= dot(S[:,i], c) * p[i] .* Y[:,i]
    end
    for i in CyclicIterator(oldest, n, rev=false)
        @views c .-= dot(Y[:,i], c) * p[i] .* S[:,i]
    end
    @views c .+= dot(S[:,newest], b) * p[newest] .* S[:,newest]
    return c
end

struct StepLBFGS{TH <: Union{Nothing, ApproxInverseHessian}}
    n::Int
    invH::TH
end
StepLBFGS(n::Int) = StepLBFGS(n, nothing)

function apply!(opt::StepLBFGS, x, Δ)
    mask = (x .> 0)

    pgrad = zeros(length(Δ))
    pgrad[mask] = Δ[mask]

    line = zeros(length(Δ))
    line[mask] = (opt.invH * pgrad)[mask]
    line[.!mask] = Δ[.!mask]
    update!(opt.invH, Δ, x)

    Δ .= line
    return Δ
  end


text:
# # Getting started

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started.ipynb)
# **Adapted from**: SOSTOOLS' SOSDEMO1 (See Section 4.1 of [SOSTOOLS User's Manual](http://sysos.eng.ox.ac.uk/sostools/sostools.pdf)) and Example 2.4 of [PJ08]
#
# P. Parrilo and A. Jadbabaie
# *Approximation of the joint spectral radius using sum of squares*.
# Linear Algebra and its Applications, Elsevier (2008), 428, 2385-2402

using Test #src
using DynamicPolynomials
@polyvar x y
p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4

# We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.
# We use `SOSModel` instead of `Model` to be able to use the `>=` syntax for Sum-of-Squares constraints.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
con_ref = @constraint(model, p >= 0)
optimize!(model)
@test primal_status(model) == MOI.FEASIBLE_POINT #src
primal_status(model)

# We see above that the solver found a feasible solution.
# We now inspect this solution:

q = gram_matrix(con_ref)
Q = getmat(q) #src
@test isapprox(Q[1, 1], 2, rtol=1e-5) #src
@test isapprox(Q[1, 2], 1, rtol=1e-5) #src
@test isapprox(Q[3, 3], 5, rtol=1e-5) #src
@test abs(Q[2, 3]) < 1e-5 #src
@test isapprox(Q[2, 2] + 2Q[1, 3], -1, rtol=1e-5) #src

# We can get the SOS decomposition from the gram matrix as follows:

sosdec = SOSDecomposition(q)
@test isapprox(sosdec, sos_decomposition(con_ref)) #src
@test isapprox(sum(sosdec.ps.^2), p; rtol=1e-4, ztol=1e-6) #src

# We now seek for the SOS decomposition of the following polynomial:

p = 4*x^4*y^6 + x^2 - x*y^2 + y^2

# We build the same model as previously with this new polynomial.
# Here we can use `Model` instead of `SOSModel` as we explicitly constrain
# `p` to belong to the SOS cone with `p in SOSCone()`.

model = Model(solver)
con_ref = @constraint(model, p in SOSCone())
optimize!(model)
@test primal_status(model) == MOI.FEASIBLE_POINT #src
primal_status(model)

# We can query the SOS decomposition direction from the constraint reference
# as follows:

sos_decomposition(con_ref)
# p should be GramMatrix([1, 0, -1/2, 0, -1, 1, 0, -2/3, 0, 4/3, 0, 0, 2, 0, 4], [y, x, x*y, x*y^2, x^2*y^3]) #src

max_stars_repo_path: examples/Getting started.jl

max_issues_repo_path: examples/Getting started.jl

max_forks_repo_path: examples/Getting started.jl


text:
## This script generates a seris of directires in which there is a relax.in
## Author: Yuewen FANG; Contact: fyuewen@gmail.com; Date: 2020 July 2nd
####original cell parameter#####
a=[   2.512269491   1.450459467   4.397688755
  -2.512269491   1.450459467   4.397688755
   0.000000000  -2.900918933   4.397688755]

for i in 93:107
    ratio = i*0.01
    round_ratio = round(ratio, digits=2) # keep two digits
    rm(string(round_ratio), recursive=true)
    mkdir(string(round_ratio))
    filename = string(string(round_ratio),"/relax.in")
    println(filename)
    fout = open(filename, "w")
    ###the part before the CELLPARAMETER###
    f = open("relax.in")
    lines = readlines(f)
    for i in lines
        if startswith(i, "   2.512269")
            break
        else
                write(fout, i)
                write(fout, '\n')
        end
    end

    new_cell_parameter = cbrt(ratio)*a
    for row in eachrow(new_cell_parameter)
        for arr_element in row
            write(fout, string(arr_element))
            write(fout, ' ')
#             print(arr_element)
#             print(' ')
        end
        write(fout, '\n')
    end

    ###the part after ATOMIC_POSITIONS###
    f = open("relax.in")
    lines = readlines(f)
    line_n = 0
    for i in lines
#        global line_n
        line_n = line_n + 1
        if startswith(i, "ATOMIC_POSITIONS")
            break
        else
#                 write(fout, i)
#                 write(fout, '\n')
        end
    end
#    print(line_n)

    line_nn = 0
    for i in lines
#        local line_nn
        line_nn = line_nn+1
        if line_nn >= line_n
            write(fout, i)
            write(fout, '\n')
        end
    end

    close(fout)
end

max_stars_repo_name: yw-fang/epwtools

max_issues_repo_name: yw-fang/epwtools

max_forks_repo_name: yw-fang/epwtools