text:
using StaticArrays
using OnlineStats
using DelimitedFiles
using CPIMC
using CPIMC.PlaneWaves
using CPIMC.Estimators
using CPIMC.UniformElectronGas
using CPIMC.DefaultUpdates

function main()
    # MC options
    NMC = 10^5
    cyc = 50
    NEquil = 2*10^4
    # system parameters
    θ = 0.125
    rs = 2.0

    # use 7 particles
    S = sphere_with_same_spin(PlaneWave((0,0,0)),dk=1)
    N = length(S)
    ξ = fractional_spin_polarization(S)
    c = Configuration(S)
    d = dimension(c.occupations)

    println("parameters:")
    println("θ: ", θ)
    println("rs: ", rs)
    println("N: ", N)
    println("ξ: ", ξ)

    e = CEnsemble(λ(N, rs, d), β(θ, N, ξ, d), N)
    equilibrate_diagonal!(UEG(), e, c)

    updates = [move_particle!, add_type_B!, remove_type_B!, add_type_C!, remove_type_C!, add_type_D!, remove_type_D!, add_type_E!, remove_type_E!, shuffle_indices!]

    measurements = Dict(# TODO: type-specification in the construction of the statistic objects (use @code_warntype)
      :sign => (Variance(), signum)
    , :Ekin => (Variance(), Ekin)
    , :W_diag => (Variance(), W_diag)
    , :K => (Variance(), K)
    , :K_fermion => (Variance(), K)
    , :occs => (Group([Variance(Float64) for i in 1:100]), occupations)
    )

    println("Start MC process ... ")
    rates = sweep!(UEG(), e, c, updates, measurements, NMC, cyc, NEquil)
    println(" finished.")

    println("parameters:")
    println("θ: ", θ)
    println("rs: ", rs)
    println("N: ", N)

    print_rates(rates)

    print_results(measurements, e)
end

main()
#Juno.@run(main())


text:
@testset "bboptimize() single-objective methods smoketest" begin
    rosenbrock2d(x) = abs2(1.0 - x[1]) + 100.0 * abs2(x[2] - x[1]^2)

    for m in keys(BlackBoxOptim.SingleObjectiveMethods)
        @testset "$(m)" begin
            ctrl = bbsetup(rosenbrock2d; Method = m,
                SearchRange = [(-5.0, 5.0), (-2.0, 2.0)], TraceMode = :silent)
            # run first iteration before the main run to exclude compilation from timing
            bboptimize(ctrl, MaxSteps = 1)
            res = bboptimize(ctrl, MaxTime = 0.3)
            @test length(best_candidate(res)) == 2
            f = best_fitness(res)
            @test typeof(f) == Float64
            @test f < 100.0 # this can't be very tight since we give very little time for optimization...
        end
    end
end

@testset "bboptimize() multi-objective methods smoketest" begin
    schaffer1(x) = (sum(abs2, x), sum(xx -> abs2(xx - 2.0), x))

    for m in keys(BlackBoxOptim.MultiObjectiveMethods)
        @testset "$(m)" begin
            ctrl = bbsetup(schaffer1; Method = m,
                SearchRange = [(-10.0, 10.0), (-10.0, 10.0)], TraceMode = :silent,
                FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true), ϵ=0.01)
            # run first iteration before the main run to exclude compilation from timing
            bboptimize(ctrl, MaxSteps = 1)
            res = bboptimize(ctrl, MaxTime = 3.0)
            @test length(best_candidate(res)) == 2
            f = best_fitness(res)
            @test typeof(f) == NTuple{2,Float64}
            @test isapprox(f[1], 2.0; atol=8E-2)
            @test isapprox(f[2], 2.0; atol=8E-2)
        end
    end
end

ext: jl


text:
module LinearRegression

using LinearAlgebra
import StatsAPI.coef
import StatsAPI.islinear

export linregress
export coef
export SolveQR, SolveCholesky

include("linreg.jl")

end


text:
function find_sources(path::String, sources=String[])
    if isdir(path)
        for entry in readdir(path)
            find_sources(joinpath(path, entry), sources)
        end
    elseif endswith(path, ".jl")
        push!(sources, path)
    end
    sources
end

function examples_testsuite()
@testset "examples" begin
    examples_dir = joinpath(@__DIR__, "..", "examples")
    examples = find_sources(examples_dir)
    filter!(file -> readline(file) != "# EXCLUDE FROM TESTING", examples)

    @testset "$(basename(example))" for example in examples
        code = """
        $(Base.load_path_setup_code())
        include($(repr(example)))
        """
        cmd = `$(Base.julia_cmd()) --startup-file=no -e $code`
        @debug "Testing $example" Text(code) cmd
        @test success(pipeline(cmd, stderr=stderr))
    end

end
end


text:
module MLToolkit

export MLToolkit

greet() = print("Welcome to Kai's machine learning toolkit!")

using ArgCheck
import Distributions, DistributionsAD, Reexport, Requires, Tracker, Parameters, LinearAlgebra, Random

### Constants

const Float = Float32
const Double = Float64

export Float, Double

### Polluting functions

Base.sum(x::AbstractArray, drop::Val{:drop}; dims=:) = dropdims(sum(x; dims=dims); dims=dims)
Base.getproperty(ts::AbstractArray{<:NamedTuple}, k::Symbol) = getproperty.(ts, k)
Base.getindex(ts::AbstractArray{<:NamedTuple}, k::Symbol) = getindex.(ts, k)

Base.getindex(arr::AbstractArray{<:Any,1}, ::Colon, i) = arr[i]
Base.getindex(arr::AbstractArray{<:Any,3}, ::Colon, i) = arr[:,:,i]
Base.getindex(arr::AbstractArray{<:Any,4}, ::Colon, i) = arr[:,:,:,i]

function Random.shuffle(nt::NamedTuple{K, V}) where {K, V<:Tuple{Vararg{<:AbstractArray}}}
    ns = [last(size(v)) for v in values(nt)]
    n = first(ns)
    @argcheck all(n .== ns)
    idcs = Random.shuffle(1:n)
    return nt[:,idcs]
end

function Base.getindex(nt::NamedTuple{K, V}, ::Colon, i) where {K, V<:Tuple{Vararg{<:AbstractArray}}}
    v = tuple([v[:, i] for v in values(nt)]...)
    return NamedTuple{K, typeof(v)}(v)
end

Base.:(+)(nt1::T, nt2::T) where {T<:NamedTuple} = T(tuple((values(nt1) .+ values(nt2))...))

function Base.union(nt1::NamedTuple{S1, T1}, nt2::NamedTuple{S2, T2}) where {S1, T1, S2, T2}
    return NamedTuple{(S1..., S2...), Tuple{T1.types..., T2.types...}}((nt1..., nt2...))
end

# TODO: clean-up `special`
include("special.jl")

# TODO: move `transformations` somewhere
include("transformations.jl")
export break_stick_ibp, break_logstick_ibp

include("PlotsX/PlotsX.jl")
Reexport.@reexport using .PlotsX
include("Datasets/Datasets.jl")
Reexport.@reexport using .Datasets
include("Scripting/Scripting.jl")
Reexport.@reexport using .Scripting

# TODO: merge `distributions` and `DistributionsX`
include("distributions/distributions.jl")
include("DistributionsX/DistributionsX.jl")
Reexport.@reexport using .DistributionsX
include("MonteCarlo/MonteCarlo.jl")
Reexport.@reexport using .MonteCarlo
# TODO: merge `neural.deprecated` (which contains codes for RAVE) and `Neural`; variational layers are also important
# include("neural.deprecated/neural.jl")

function __init__()
    Requires.@require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        include("gpu.jl")
        export seed!
    end

    Requires.@require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
        include("Neural/Neural.jl")
        Reexport.@reexport using .Neural
    end
end

end # module


text:
# parameters that govern counterfactual

ImpvolEquilibrium.remove_shock!(parameters, :global_sectoral_shock_njs)
ImpvolEquilibrium.remove_shock!(parameters, :idiosyncratic_shock_njs)


text:
"""
	MultipleRepresentation(extractors::Tuple)

Extractor extracts item to a `ProductNode` where each item is different extractor and item is
extracted by all extractors in multirepresentation.

# Examples

## Example of both categorical and string representation

One of usecases is to use string representation for strings and categorical variable representation for most frequent values.
This allows model to more easily learn frequent or somehow else significant values, which creating meaningful representation
for previously unseen inputs.

```jldoctest
julia> e = MultipleRepresentation((ExtractString(false), ExtractCategorical(["tcp", "udp", "dhcp"], false)));

julia> s1 = e("tcp")
ProductNode with 1 obs
  ├── e1: ArrayNode(2053×1 NGramMatrix with Int64 elements) with 1 obs
  └── e2: ArrayNode(4×1 OneHotArray with Bool elements) with 1 obs

julia> s1[:e1]
2053×1 Mill.ArrayNode{Mill.NGramMatrix{String, Vector{String}, Int64}, Nothing}:
 "tcp"

julia> s1[:e2]
4×1 Mill.ArrayNode{Flux.OneHotArray{Int64, 4, 1, 2, Vector{Int64}}, Nothing}:
 0
 1
 0
 0

julia> s2 = e("http")
ProductNode with 1 obs
  ├── e1: ArrayNode(2053×1 NGramMatrix with Int64 elements) with 1 obs
  └── e2: ArrayNode(4×1 OneHotArray with Bool elements) with 1 obs

julia> s2[:e1]
2053×1 Mill.ArrayNode{Mill.NGramMatrix{String, Vector{String}, Int64}, Nothing}:
 "http"

julia> s2[:e2]
4×1 Mill.ArrayNode{Flux.OneHotArray{Int64, 4, 1, 2, Vector{Int64}}, Nothing}:
 0
 0
 0
 1
```

## Example of irregular schema representation

The other usecase is to handle irregular schema, where extractor returns `missing` representation if it's unable to
extract it properly. Of course there do not have to be only leaf value extractors, some extractors may be ExtractDict,
while other are extracting leaves etc.

```jldoctest
julia> e = MultipleRepresentation((ExtractString(), ExtractScalar(Float32, 2, 3)));

julia> s1 = e(5)
ProductNode with 1 obs
  ├── e1: ArrayNode(2053×1 NGramMatrix with Union{Missing, Int64} elements) with 1 obs
  └── e2: ArrayNode(1×1 Array with Union{Missing, Float32} elements) with 1 obs

julia> s1[:e1]
2053×1 Mill.ArrayNode{Mill.NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}, Nothing}:
 missing

julia> s1[:e2]
1×1 Mill.ArrayNode{Matrix{Union{Missing, Float32}}, Nothing}:
 9.0f0

julia> s2 = e("hi")
ProductNode with 1 obs
  ├── e1: ArrayNode(2053×1 NGramMatrix with Union{Missing, Int64} elements) with 1 obs
  └── e2: ArrayNode(1×1 Array with Union{Missing, Float32} elements) with 1 obs

julia> s2[:e1]
2053×1 Mill.ArrayNode{Mill.NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}, Nothing}:
 "hi"

julia> s2[:e2]
1×1 Mill.ArrayNode{Matrix{Union{Missing, Float32}}, Nothing}:
 missing

```
"""
struct MultipleRepresentation{E<:Union{NamedTuple, Tuple}}
	extractors::E
end

MultipleRepresentation(vs::Vector) = MultipleRepresentation((;[Symbol("e$(i)") => v for (i,v) in enumerate(vs)]...))
MultipleRepresentation(vs::Tuple) = MultipleRepresentation((;[Symbol("e$(i)") => v for (i,v) in enumerate(vs)]...))
(m::MultipleRepresentation)(x; store_input=false) = ProductNode(map(e -> e(x; store_input), m.extractors))
Base.keys(e::MultipleRepresentation) = keys(e.extractors)

Base.hash(e::MultipleRepresentation, h::UInt) = hash(e.extractors, h)
Base.:(==)(e1::MultipleRepresentation, e2::MultipleRepresentation) = e1.extractors == e2.extractors
Base.getindex(e::MultipleRepresentation, i::Int) = e.extractors[i]


text:
@testset "1.2.2 Quartic" begin
    include("1.2.2.2 (d x)^m (a+b x^2+c x^4)^p.jl")
    include("1.2.2.3 (d+e x^2)^m (a+b x^2+c x^4)^p.jl")
    include("1.2.2.4 (f x)^m (d+e x^2)^q (a+b x^2+c x^4)^p.jl")
    include("1.2.2.5 P(x) (a+b x^2+c x^4)^p.jl")
    include("1.2.2.6 P(x) (d x)^m (a+b x^2+c x^4)^p.jl")
    include("1.2.2.7 P(x) (d+e x^2)^q (a+b x^2+c x^4)^p.jl")
    include("1.2.2.8 P(x) (d+e x)^q (a+b x^2+c x^4)^p.jl")
end

max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/.jl


text:
function calculaAnguloRotacion(x,y)

x=x[:,1];
y=y[:,1];

numLados=size(x,1);
largoLados=zeros(numLados,1);
theta=zeros(numLados,1);
cr=zeros(numLados,2);
for i=1:numLados
    if i<numLados
        x1=x[i];
        y1=y[i];
        x2=x[i+1];
        y2=y[i+1];
    else
        x1=x[end];
        y1=y[end];
        x2=x[1];
        y2=y[1];
    end
    d=sqrt((x2-x1)^2+(y2-y1)^2);
    largoLados[i,1]=d;
    theta1=atan((y1-y2)/(x1-x2));
    theta2=abs(theta1)-pi/2;
    if abs(theta1)<=abs(theta2)
        theta[i]=theta1;
    else
        theta[i]=theta2;
    end
    cr[i,:]=[x1 y1];
end
largoLados=largoLados[:,1];
isort = sortperm(largoLados, rev=true);
theta=theta[isort[1:numLados]];
cr=cr[isort[1:numLados],:];

theta, cr

end


text:
"""
svd_entropy(A::Matrix{Bool})
    A: Adjacency matrix of a unipartite directed simple network
Returns the SVD-entropy of the adjacency matrix of a deterministic network using a Singular Value Decomposition
"""
function svd_entropy(A::Matrix{Bool})
    # compute the SVD-entropy
    F = svd(A)
    Λ = F.S[1:rank(A)]
    λ = Λ ./ sum(Λ)
    return -sum(λ .* log.(λ)) * 1 / log(length(λ))
end

"""
network_maxent(N::UnipartiteNetwork, nsteps::Int64, T::Float64)
    N: Unipartite directed simple network
    nsteps: Number of steps
    T: Temperature of the simulating annealing algorithm
Returns the adjacency matrix of maximum SVD-entropy constrained by the joint degree sequence of N using a simulating annealing algorithm
"""
function network_maxent(N::UnipartiteNetwork, nsteps::Int64, T::Float64)
    A = convert(Matrix, N.edges) # convert UnipartiteNetwork to dense matrix format
    rmg = matrixrandomizer(A) # matrix generator object

    # initial random matrix that has the same joint degree sequence as N
    A0 = convert(Matrix{Bool}, rand(rmg))

    # simulating annealing algorithm
    for i in 1:nsteps
        # propose a new constrained random matrix and compute the difference in SVD-entropy
        A1 = convert(Matrix{Bool}, rand(rmg))
        delta_entropy = svd_entropy(A1) - svd_entropy(A0)
        # accept if the difference is positive or with a probability p if it's negative
        if delta_entropy > 0
            A0 = A1
        else
            p = exp(delta_entropy/T)
            P = rand(Uniform(0, 1))
            if P < p
                A0 = A1
            end
        end

        # update the temperature
        T <- T * exp(-100/nsim)
    end
    max_entropy = svd_entropy(A0)
    Amax = (A = A0, entropy = max_entropy)
    return Amax
end


text:
"""
    `NonAnalyticTest`
An internal model used primarly for testing. Any model passed to it will be interpreted
as not having an analytic Fourier transform.
"""
struct NonAnalyticTest{M} <: AbstractModel
    model::M
end

ComradeBase.visanalytic(::Type{<:NonAnalyticTest}) = NotAnalytic()
ComradeBase.isprimitive(::Type{<:NonAnalyticTest}) = IsPrimitive()

@inline radialextent(m::NonAnalyticTest) = radialextent(m.model)
@inline intensity_point(m::NonAnalyticTest, x, y) = intensity_point(m.model, x, y)


text:
module NeXLDatabase

using SQLite
using Dates
using DataFrames
using IntervalSets
using Reexport
using Requires

@reexport using NeXLSpectrum

export openNeXLDatabase

function openNeXLDatabase(filename::AbstractString)::SQLite.DB
    function stripcomment(ss)
        i=findfirst("--",ss)
        return i == nothing ? ss : ss[1:i.start-1]
    end
    function buildTable(tbl)
        path = dirname(pathof(@__MODULE__))
        cmds = reduce(*, map(stripcomment, readlines("$(path)/sql/$(tbl).sql")))
        for cmd in strip.(split(cmds,";"))
            if length(cmd)>0
                SQLite.execute(db, cmd)
            end
        end
    end
    db = SQLite.DB(filename)
    existing = SQLite.tables(db)
    tables = (
        "material", "massfraction", #
        "person", "laboratory", "labmember", #
        "sample", "project", #
        "instrument", "detector", #
        "artifact", "spectrum", "fitspectrum", #
        "referencespectrum", "campaign", #
        "kratio", "standardfor", "campaigncomment" #
    )
    for tbl in tables
        if (!haskey(existing, :name)) || (!(uppercase(tbl) in existing.name))
            # @info "Creating database table $(tbl)."
            buildTable(tbl)
        end
    end
    return db
end

include("material.jl")
include("person.jl")
include("laboratory.jl")
include("instrument.jl")
include("sample.jl")
include("artifact.jl")
include("project.jl")
include("spectrum.jl")
include("referencespectrum.jl")
include("fitspectrum.jl")
include("campaign.jl")
include("campaigncomment.jl")
include("kratio.jl")
include("standardfor.jl")

export DBPerson
export DBLaboratory
export DBInstrument
export DBDetector
export DBArtifact
export DBSpectrum
export DBMember
export DBSample
export DBProject
export DBProjectSpectrum
export find
export DBCampaign
export DBFitSpectrum
export DBReference
export DBKRatio
export DBStandardFor
export DBFSComment

export measured, references, dbreferences
export constructFitSpectra
export disposition, disposition!

function __init__()
    @require Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004" include("gadflysupport.jl")
end

end # module


text:
#
struct LinearSolveFunction{F} <: AbstractSolveFunction
    solve_func::F
end

function SciMLBase.solve(cache::LinearCache, alg::LinearSolveFunction,
                         args...; kwargs...)
    @unpack A,b,u,p,isfresh,Pl,Pr,cacheval = cache
    @unpack solve_func = alg

    u = solve_func(A,b,u,p,isfresh,Pl,Pr,cacheval;kwargs...)
    cache = set_u(cache, u)

    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
end


text:
module TestTextAnalysis
using SparseArrays
using Test
using Languages
using TextAnalysis
using Compat

# @testset "TextAnalysis" begin

println("Running tests:")
println(typeof(Compat.String))

include("tokenizer.jl")
include("ngramizer.jl")
include("document.jl")
include("metadata.jl")
include("corpus.jl")
include("preprocessing.jl")
include("dtm.jl")
include("stemmer.jl")
include("tf_idf.jl")
include("lda.jl")
include("summarizer.jl")
include("sentiment.jl")

# end
end


text:
using Printf

function grouppoints(Xcopy)

    X = deepcopy(Xcopy)

    K = length(X)

    #------------------------------------------------
    # find overall earliest and latest point in time
    #------------------------------------------------

    earliestK = 0
    earliestT = Inf

    latestK = 0
    latestT = -Inf

    for k in 1:K

        # make sure it is sorted
        @assert(issorted(vec(X[k][:,1])))

        if X[k][1,1] < earliestT
            earliestT = X[k][1,1]
            earliestK = k
        end

        if X[k][end,1] > latestT
            latestT = X[k][end,1]
            latestK = k
        end

    end

    @printf("Earliest time is %5.3f\n", earliestT)
    @printf("Latest   time is %5.3f\n", latestT)

    #------------------------------------------------
    # group points
    #------------------------------------------------

    gridpoints = ceil(Int, (latestT - earliestT) * 10)

    tol = 1.0

    Grouped = zeros(0, K)
    T       = zeros(0)

    for t in LinRange(earliestT, latestT, gridpoints)

        closestintime = [minimum(abs.(X[k][:,1] .- t)) for k in 1:K]

        if all(closestintime .< tol)

            idx = [argmin(abs.(X[k][:,1] .- t)) for k in 1:K]

            @show t idx

            Grouped = [Grouped; reshape([X[k][idx[k], 2] for k in 1:K], 1, K)]

            T = [T;t]

            # delete used elements from respective arrays

            for k in 1:K
                X[k] = X[k][setdiff(1:size(X[k],1), idx[k]), :]

                if isempty(X[k])
                    return T, Grouped
                end

            end

        end

    end

    return T, Grouped

end


text:
println("Master node pid: $(getpid())"); flush(stdout)

# Load libraries on master node.
# NOTE: In theory, this step is not needed. In practice, this step guarantees
# that libraries are compiled and cached before they are used by the workers.
# If this is not done, workers may try to compile the same library
# simultaneously and write the compiled cache to the same locatoin.
println("Compile libraries on master node ..."); flush(stdout)
include("imports.jl")
println("Finished compiling libraries on master node ..."); flush(stdout)

using Distributed

# Parse command line args
if length(ARGS) == 0
  println("Running in debug mode! This does not reproduce results in paper!")
  flush(stdout)

  num_proc = 4
  ps = [0.85, 0.99]
  SEEDS = [2, 4]
  NITERS = 2000
  BATCHSIZE = 200
  RESULTS_DIR = "results/test"
else
  num_proc = parse(Int, ARGS[1])  # 40
  ps = [0.85, 0.9, 0.95, 0.99]
  SEEDS = collect(1:10)
  NITERS = 20000
  # BATCHSIZE = 2000
  BATCHSIZE = parse(Int, ARGS[2])  # 2000, 500
  RESULTS_DIR = ARGS[3]
end

# Remove old workers
rmprocs(filter(w -> w > 1, workers()))

# Add processeors
println("Adding $(num_proc) workers ..."); flush(stdout)
addprocs(num_proc)

@everywhere RESULTS_DIR = $(RESULTS_DIR)
@everywhere NITERS = $(NITERS)
@everywhere BATCHSIZE = $(BATCHSIZE)
@everywhere function separatesamples(Y)
  sample_ids = sort(unique(Y.sample_id))
  return Matrix.(Y[Y.sample_id .== i, DF.Not(:sample_id)]
                 for i in sample_ids)
end

# Print worker process ids
worker_pids = pmap(id -> getpid(), workers())
println("Worker pids: $(join(worker_pids, " "))"); flush(stdout)

# Load the same libraries on worker nodes.
println("Loading libraries on worker nodes ..."); flush(stdout)
@everywhere include("imports.jl")
println("Finished loading libraries on workers node ..."); flush(stdout)

@everywhere function fit(settings)
  seed = settings[:seed]
  p = settings[:p]

  # Read data
  y = let
    path_to_data = "../data/cb_transformed_reduced_p$(p).csv"
    Y = coalesce.(CSV.read(path_to_data), NaN)
    separatesamples(Y)
  end

  # Results directory
  results_dir = "$(RESULTS_DIR)/p$(p)_seed$(seed)"
  mkpath(results_dir)

  # Redirect output
  CytofResearch.Model.redirect_stdout_to_file("$(results_dir)/log.txt") do
    # Set random seed for reproducibiliy
    Random.seed!(seed)
    println("seed: $(seed)")

    # Generate model constants
    c = CytofResearch.VB.Constants(y=y, K=30, L=Dict(false => 5, true => 3),
                                   yQuantiles=[.0, .25, .5], pBounds=[.05, .8, .05],
                                   use_stickbreak=false, tau=.001)
    c.priors.eps = Beta(1, 99)
    c.priors.sig2 = Gamma(.1, 1)

    # Fit model
    out = CytofResearch.VB.fit(y=y, niters=NITERS, batchsize=BATCHSIZE, c=c,
                               nsave=30, seed=seed)

    # Save results
    BSON.bson("$(results_dir)/output.bson", out)
  end  # of redirect

  return
end

# Run the thing in parallel
settings_vec = [(seed=seed, p=p) for seed in SEEDS for p in ps]
println("Number of settings: $(length(settings_vec))")
status = pmap(fit, settings_vec, on_error=identity)
status_sanitized = map(s -> s == nothing ? "success" : s, status)

# Printing success / failure status of runs.
for (settings, s) in zip(settings_vec, status_sanitized)
  println("$(settings) => status: $(s)")
end

println("ALL DONE!")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
"""
# Summary
    mutable struct ZNorm <: Normalizer
# Fields
    γ        :: VarOrNil                        # scaling params
    β        :: VarOrNil                        # shifting params
    μ        :: Union{AbstractArray,Nothing}    # running average
    σ        :: Union{AbstractArray,Nothing}    # running variance otherwise standard deviation
    views    :: NTuple                          # views to collect elements for mean and var
    training :: Bool                            # if traning or not
    epsilion :: AbstractFloat                   # prevent dividing by zero, 1e-10 for default
    momentum :: AbstractFloat                   # smoothing const, or called historical inertia coefficient

Applies mean and scaling normalization over a N-dimensional input, like BatchNorm LayerNorm and InstanceNorm.

"""
mutable struct ZNorm <: Normalizer
    γ::VarOrNil                        # scaling params
    β::VarOrNil                        # shifting params
    μ::Union{AbstractArray,Nothing}    # running average
    σ::Union{AbstractArray,Nothing}    # running variance otherwise standard deviation
    views::NTuple                      # views to collect elements for mean and var
    training::Bool                     # if traning or not
    epsilion::AbstractFloat            # prevent dividing by zero, 1e-10 for default
    momentum::AbstractFloat            # inertia coefficient
    function ZNorm(;ndims::Int,        # how many dimentions the input data has
                   keptdims::Union{Tuple,Int},     # must be unique and sorted and positive
                   keptsize::Union{Tuple,Int},     # must be positive
                   epsilion::AbstractFloat=1e-10,  # stability const
                   momentum::AbstractFloat=0.900,  # smoothing const or historical inertia
                   type::Type=Array{Float32})

        shape, views = ShapeAndViews(ndims, keptdims, keptsize);
        γ = Variable{type}( Ones(type, shape), true, true, true);
        β = Variable{type}(Zeros(type, shape), true, true, true);
        μ = Zeros(type, shape);
        σ =  Ones(type, shape);
        new(γ, β, μ, σ, views, true, epsilion, momentum)
    end
    function ZNorm(views, training, epsilion, momentum)
        new(nothing, nothing, nothing, nothing, views, training, epsilion, momentum)
    end
end

function clone(this::ZNorm; type::Type=Array{Float32})
    cloned = ZNorm(this.views, this.training, this.epsilion, this.momentum)
    cloned.γ = clone(this.γ, type=type)
    cloned.β = clone(this.β, type=type)
    cloned.μ = type(this.μ)
    cloned.σ = type(this.σ)
    return cloned
end

function Base.show(io::IO, m::ZNorm)
    SIZE = size(m.β.value)
    TYPE = typeof(m.β.value)
    print(io, "ZNorm(statistic vars' size=$SIZE; type=$TYPE)")
end

function paramsof(m::ZNorm)
    params = Vector{Variable}(undef,2)
    params[1] = m.γ
    params[2] = m.β
    return params
end

function xparamsof(m::ZNorm)
    xparams = Vector{XVariable}(undef,2)
    xparams[1] = ('w', m.γ)
    xparams[2] = ('b', m.β)
    return xparams
end

function nparamsof(model::ZNorm)
    return 4*length(model.β)
end

function bytesof(model::ZNorm, unit::String="MB")
    n = nparamsof(model) * sizeof(eltype(model.β))
    return blocksize(n, uppercase(unit))
end

function forward(b::ZNorm, x::Variable{T}) where T
    γ = b.γ
    β = b.β
    ϵ = b.epsilion
    ρ = b.momentum
    v = b.views
    μ = mean(ᵛ(x), dims=v)
    σ =  std(ᵛ(x), dims=v, mean=μ, corrected=false)
    𝐗 = (ᵛ(x) .- μ) ./ (σ .+ ϵ)
    y = Variable{T}(𝐗 .* ᵛ(γ) .+ ᵛ(β), x.backprop)

    if y.backprop
        Σ = σ .* σ
        @. b.μ = ρ * b.μ + (1 - ρ) * μ    # running mean
        @. b.σ = ρ * b.σ + (1 - ρ) * Σ    # running var
        function ZNormBackward()
            if need2computeδ!(x)
                n     = length(σ)/length(x)
                σ¯¹   = 1 ./ (σ .+ ϵ)
                σ¯³   = (σ¯¹).^3
                Δ     = ᵛ(x) .- μ
                ∂𝐋∂𝐗  = δ(y) .* ᵛ(γ)
                Δ∂𝐋∂𝐗 = Δ .* ∂𝐋∂𝐗
                SumΔ∂𝐋∂𝐗  = sum(Δ∂𝐋∂𝐗, dims=v)

                δ(x) .+= σ¯¹ .* ∂𝐋∂𝐗
                δ(x) .-= σ¯³ .* n   .* SumΔ∂𝐋∂𝐗 .* Δ
                δ(x) .+= σ¯³ .* n^2 .* SumΔ∂𝐋∂𝐗 .* sum(Δ, dims=v) .- σ¯¹ .* n .* sum(∂𝐋∂𝐗, dims=v)

                if need2computeδ!(γ) δ(γ) .+= sum(δ(y) .* 𝐗, dims=v) end
                if need2computeδ!(β) δ(β) .+= sum(δ(y),      dims=v) end
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, ZNormBackward)
    end
    return y
end

function predict(b::ZNorm, x::AbstractArray)
    ϵ = b.epsilion
    γ = ᵛ(b.γ)
    β = ᵛ(b.β)
    μ = b.μ
    σ = b.σ
    return @. (x - μ) / sqrt(σ + ϵ) * γ + β
end

function BatchNorm1d(nchannels::Int;
            epsilion::AbstractFloat=1e-10,
            momentum::AbstractFloat=0.900,
            type::Type=Array{Float32})
    return ZNorm(ndims=3,
                 keptdims=1,
                 keptsize=nchannels,
                 epsilion=epsilion,
                 momentum=momentum,
                 type=type)
end

function BatchNorm0d(nchannels::Int;
            epsilion::AbstractFloat=1e-10,
            momentum::AbstractFloat=0.900,
            type::Type=Array{Float32})
    return ZNorm(ndims=2,
                 keptdims=1,
                 keptsize=nchannels,
                 epsilion=epsilion,
                 momentum=momentum,
                 type=type)
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# Autogenerated wrapper script for YASM_jll for x86_64-w64-mingw32
export vsyasm, yasm, ytasm

using NASM_jll
JLLWrappers.@generate_wrapper_header("YASM")
JLLWrappers.@declare_executable_product(vsyasm)
JLLWrappers.@declare_executable_product(yasm)
JLLWrappers.@declare_executable_product(ytasm)
function __init__()
    JLLWrappers.@generate_init_header(NASM_jll)
    JLLWrappers.@init_executable_product(
        vsyasm,
        "bin\\vsyasm.exe",
    )

    JLLWrappers.@init_executable_product(
        yasm,
        "bin\\yasm.exe",
    )

    JLLWrappers.@init_executable_product(
        ytasm,
        "bin\\ytasm.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("../../common/numder/src/numder_module.jl")
# include("testfunctions.jl")

module opti
using ..numder
using LinearAlgebra

include("opti_1d.jl")
include("line_search.jl")
include("gradient_descent.jl")

mutable struct BFGS <: DescentMethod
    B
end
function BFGS(x::Vector{T}) where {T}
    BFGS( Matrix(one(T)*I, length(x), length(x)) )
end

function init!(m::BFGS, x)
    m.B = Matrix(one(eltype(x))*I, length(x), length(x))
end

function bfgs_matrix_update(B, g, g_prev, x, x_prev)
    dg = g - g_prev
    dx = x - x_prev

    dxdg = dot(dx, dg)

    if dxdg < 0
        println("No update happens, some accuracy errors")
        println(" dxdg = $(dxdg) ")
        return B
    end

    t1 = (dx ./ dxdg) * (dg' * B)
    t1 = t1 + t1'

    t2f = 1 + (dg' * (B * dg)) / dxdg
    t2 = (dx .* (t2f / dxdg)) * dx'

    # B = B - t1 + t2
    B += (t2 - t1)
    return B
end

function quasi_newton(
    fcn,
    x0;
    tol_options::ToleranceOptions       = ToleranceOptions(),
    search_options::LineSearchOptions   = LineSearchOptions(),
    gradient_fcn                        = numder.gradient_fw,
    line_search_fcn                     = line_search,
    log_path                            = false)

    m = BFGS(x0)

    x = x0
    f = fcn(x)
    g = gradient_fcn(fcn, x)

    step_direction  = -m.B * g
    step_size       = line_search_fcn(fcn, x, step_direction, f, search_options)
    dx              = step_size * step_direction
    g_new           = gradient_fcn(fcn, x + dx)
    dg              = g_new - g
    m.B             = m.B * dot(dg, dx) / dot(dg, dg)

    x_prev  = x
    x       = x_prev + dx
    f_prev  = f
    f       = fcn(x)
    g_prev  = g
    g       = g_new

    # out parameters
    x_sol = x
    y_sol = f
    stopping_crit = unknown::GradientDescentStoppingCrit

    iter = 0
    for outer iter = 1:tol_options.max_iter

        # update B
        # dg = g - g_prev
        # dx = x - x_prev

        # dxdg = dot(dx, dg)

        # t1 = (dx ./ dxdg) * (dg' * m.B)
        # t1 = t1 + t1'

        # t2f = 1 + (dg' * (m.B * dg)) / dxdg
        # t2 = (dx .* (t2f / dxdg)) * dx'

        # m.B = m.B - t1 + t2

        m.B = bfgs_matrix_update(m.B, g, g_prev, x, x_prev)

        # step dir
        step_direction = - m.B * g
        step_size      = line_search_fcn(fcn, x, step_direction, f, search_options)

        dx = step_direction * step_size

        # stopping criterion based on the step size
        abs_step_size = norm( dx )
        if abs_step_size <= tol_options.step_size_tol
            x_sol = x
            y_sol = f
            stopping_crit = small_step_size::GradientDescentStoppingCrit
            break
        end

        x_prev = x
        x      = x_prev + dx
        f_prev = f
        f      = fcn(x)
        g_prev = g

        # stopping criteria for function value
        if abs(f - f_prev) <= tol_options.f_abs_tol || abs(f - f_prev) / ( abs(f_prev) == 0 ? mach_eps : abs(f_prev)) < tol_options.f_rel_tol
            x_sol = x
            y_sol = f
            stopping_crit = small_function_value_delta::GradientDescentStoppingCrit
            break
        end

        # calculate the gradient at the new point
        g = gradient_fcn(fcn, x)

        # stopping criterion based on the gradient norm
        norm_g = norm(g)
        if norm_g <= tol_options.g_abs_tol
            x_sol = m.x
            y_sol = m.f
            stopping_crit = small_gradient::GradientDescentStoppingCrit
            break
        end

    end

    return (x_sol, y_sol, iter, stopping_crit)

end

function bfgs(
	fcn,
	x0;
	grad_fcn = numder.gradient_fw,
	options = ToleranceOptions() )

    n_dim = length(x0)
    H     = Matrix{Float64}(I, n_dim, n_dim)
    eye   = Matrix{Float64}(I, n_dim, n_dim)

    # first step is the steepest gradient descent
    g0        = grad_fcn(fcn, x0)
    step_size = line_search(fcn, x0, -g0)
    x1        = x0 - step_size * g0
    g1        = grad_fcn(fcn, x1)

    y            = g1 - g0
    s            = -step_size * g0
    scale_factor = dot(y, s) / dot(y, y)
    for ii = 1:n_dim
        H[ii, ii] *= scale_factor
    end

    # main loop
    iter = 0
    for outer iter = 1:options.max_iter

        # update of the matrix
        s       = x1 - x0
        y       = g1 - g0

		sdy     = dot(s, y)
		Hs 		= H * s
		sHs     = dot(s, Hs)
		theta   = 1.0
		if sdy < 0.2 * sHs
			theta = 0.8 * sHs / (sHs - sdy)
		end

		r = theta * y + (1 - theta) * Hs

        Hy      = H * y
		rho     = 1 / sdy
        rhos    = rho * s
        Hyst    = Hy * rhos'
        sytHyst = rhos * (y' * Hyst)
        # sst     = rhos * s'
        sst 	= r * r' / dot(s, r)

		Hnew    = H - (Hyst + Hyst') + sytHyst + sst
        # Hnew   = (eye - rho * s * y') * H * (eye - rho * y * s') + rho * s * s'

        # dH = norm(Htest - Hnew)
        # println(dH)

        x0        = x1
        p         = -(Hnew * g1)
        step_size = line_search(fcn, x0, p)
        dx        = p * step_size
        x1        = x0 + dx

        if norm(dx) <= options.step_size_tol
            break
        end

        # options.f_abs_tol
        # options.f_rel_tol

        g0        = g1
        g1        = grad_fcn(fcn, x1)

        if norm(g1) <= options.g_abs_tol
            break
        end

        H = Hnew
    end

    return (x1, iter)
end

end # opti

#=

function nonlinear_conjugate_gradient_descent_template(beta_method, fcn, x0, options, wolfe_condition_options)
	x_act = x0
	y_act = fcn(x_act)
	x_next = x_act
	x_sol = x_act

	# first step is a gradient descent step
	g_act = NumDiff.gradient_fw(fcn, x_act)
	direction_act = -g_act
	alpha = strong_wolfe_line_search(fcn, x_act, direction_act,
	 								wolfe_condition_options.alpha_0,
									wolfe_condition_options.beta,
									wolfe_condition_options.sigma)
	step = alpha * direction_act
	x_act = x_act + step
	y_act = fcn(x_act)
	g_prev = g_act
	direction_prev = direction_act

	# main loop
	println("x: $x_act")
	iter = 1
	for outer iter = 1:options.max_iter

		g_act = NumDiff.gradient_fw(fcn, x_act)

		## specific method is applied
		beta = beta_method(g_act, g_prev, direction_prev)

		direction_act = -g_act + beta .* direction_prev

		alpha = strong_wolfe_line_search(fcn, x_act, direction_act)
		step = alpha * direction_act
		x_next = x_act + step
		y_next = fcn(x_next)

		# convergence checks
		if norm(g_act) <= options.gradient_tolerance
			x_sol = x_next
			break
		end
		if norm(step) <= options.step_tolerance
			x_sol = x_next
			break
		end
		if abs(y_next - y_act) <= options.function_value_change_tolerance
			x_sol = x_next
			break
		end

		# updates
		y_act = y_next
		x_act = x_next
		direction_prev = direction_act
		println("x: $x_act")
	end

	if iter == options.max_iter
		x_sol = x_next
	end
	return x_sol
end

function beta_polak_ribiere(g_act, g_prev, direction_prev)
	y_hat = g_act - g_prev
	beta = dot(g_act, y_hat) / dot( g_prev, g_prev )
	beta = max(beta, 0)
	return beta
end

function beta_dai_yuan(g_act, g_prev, direction_prev)
	y_hat = g_act - g_prev
	beta = dot(g_act, g_act) / dot(y_hat, direction_prev)
	beta = max(beta, 0)
	return beta
end

function beta_hestenes_stiefel(g_act, g_prev, direction_prev)
	y_hat = g_act - g_prev
	beta = dot(y_hat, g_act) / dot(y_hat, direction_prev)
	beta = max(beta, 0)
	return beta
end

function beta_hager_zhang(g_act, g_prev, direction_prev)
	# probably implementation error
	y_hat = g_act - g_prev
	numer = dot(y_hat, direction_prev)
	t1 = y_hat - 2 .* direction_prev .* dot(y_hat, y_hat) ./ numer
	t2 = g_act ./ numer
	beta = dot(t1, t2)
	# beta = max(beta, 0)
	return beta
end

function nonlinear_conjugate_gradient_descent_polak_ribiere(
	fcn,
	x0,
	options::GradientDescentOptions = default_options(GradientDescentOptions),
	wolfe_condition_options = default_strong_wolfe_conditions_options()
	)
	x_sol = nonlinear_conjugate_gradient_descent_template(beta_polak_ribiere, fcn, x0, options, wolfe_condition_options)
end

function nonlinear_conjugate_gradient_descent_dai_yuan(
	fcn,
	x0,
	options::GradientDescentOptions = default_options(GradientDescentOptions),
	wolfe_condition_options = default_strong_wolfe_conditions_options()
	)
	x_sol = nonlinear_conjugate_gradient_descent_template(beta_dai_yuan, fcn, x0, options, wolfe_condition_options)
end

function nonlinear_conjugate_gradient_descent_hestenes_stiefel(
	fcn,
	x0,
	options::GradientDescentOptions = default_options(GradientDescentOptions),
	wolfe_condition_options = default_strong_wolfe_conditions_options()
	)
	x_sol = nonlinear_conjugate_gradient_descent_template(beta_hestenes_stiefel, fcn, x0, options, wolfe_condition_options)
end

function nonlinear_conjugate_gradient_descent_hager_zhang(
	fcn,
	x0,
	options::GradientDescentOptions = default_options(GradientDescentOptions),
	wolfe_condition_options = default_strong_wolfe_conditions_options()
	)
	x_sol = nonlinear_conjugate_gradient_descent_template(beta_hager_zhang, fcn, x0, options, wolfe_condition_options)
end

=#

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
function getindex(cur_matrix::FourBlockMatrix,i::Int,j::Int)
  n = cur_matrix.n
  m = cur_matrix.m

  @assert i <= n + m
  @assert j <= n + m

  if i <= n
    if j <= n
      cur_value = cur_matrix.A[i,j]
    else
      cur_value = cur_matrix.B[i,j-n]
    end
  else
    if j <= n
      cur_value = cur_matrix.C[i-n,j]
    else
      cur_value = cur_matrix.D[i-n,j-n]
    end
  end

  cur_value
end