text:
module SKLearners

export SKLearner,transform!,fit!
export skkrun

using TSML.TSMLTypes
import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload
using TSML.Utils

using PyCall

const ENS=pyimport("sklearn.ensemble")
const LM=pyimport("sklearn.linear_model")
const DA=pyimport("sklearn.discriminant_analysis")
const NN=pyimport("sklearn.neighbors")
const SVM=pyimport("sklearn.svm")
const TREE=pyimport("sklearn.tree")
const ANN=pyimport("sklearn.neural_network")
const GP=pyimport("sklearn.gaussian_process")
const KR=pyimport("sklearn.kernel_ridge")
const NB=pyimport("sklearn.naive_bayes")
const ISO=pyimport("sklearn.isotonic")

# Available scikit-learn learners.
learner_dict = Dict(
  "AdaBoostClassifier" => ENS.AdaBoostClassifier,
  "BaggingClassifier" => ENS.BaggingClassifier,
  "ExtraTreesClassifier" => ENS.ExtraTreesClassifier,
  "VotingClassifier" => ENS.VotingClassifier,
  "GradientBoostingClassifier" => ENS.GradientBoostingClassifier,
  "RandomForestClassifier" => ENS.RandomForestClassifier,
  "LDA" => DA.LinearDiscriminantAnalysis,
  "QDA" => DA.QuadraticDiscriminantAnalysis,
  "LogisticRegression" => LM.LogisticRegression,
  "PassiveAggressiveClassifier" => LM.PassiveAggressiveClassifier,
  "RidgeClassifier" => LM.RidgeClassifier,
  "RidgeClassifierCV" => LM.RidgeClassifierCV,
  "SGDClassifier" => LM.SGDClassifier,
  "KNeighborsClassifier" => NN.KNeighborsClassifier,
  "RadiusNeighborsClassifier" => NN.RadiusNeighborsClassifier,
  "NearestCentroid" => NN.NearestCentroid,
  "SVC" => SVM.SVC,
  "LinearSVC" => SVM.LinearSVC,
  "NuSVC" => SVM.NuSVC,
  "MLPClassifier" => ANN.MLPClassifier,
  "GaussianProcessClassifier" => GP.GaussianProcessClassifier,
  "DecisionTreeClassifier" => TREE.DecisionTreeClassifier,
  "GaussianNB" => NB.GaussianNB,
  "MultinomialNB" => NB.MultinomialNB,
  "ComplementNB" => NB.ComplementNB,
  "BernoulliNB" => NB.BernoulliNB,
  "SVR" => SVM.SVR,
  "Ridge" => LM.Ridge,
  "RidgeCV" => LM.RidgeCV,
  "Lasso" => LM.Lasso,
  "ElasticNet" => LM.ElasticNet,
  "Lars" => LM.Lars,
  "LassoLars" => LM.LassoLars,
  "OrthogonalMatchingPursuit" => LM.OrthogonalMatchingPursuit,
  "BayesianRidge" => LM.BayesianRidge,
  "ARDRegression" => LM.ARDRegression,
  "SGDRegressor" => LM.SGDRegressor,
  "PassiveAggressiveRegressor" => LM.PassiveAggressiveRegressor,
  "KernelRidge" => KR.KernelRidge,
  "KNeighborsRegressor" => NN.KNeighborsRegressor,
  "RadiusNeighborsRegressor" => NN.RadiusNeighborsRegressor,
  "GaussianProcessRegressor" => GP.GaussianProcessRegressor,
  "DecisionTreeRegressor" => TREE.DecisionTreeRegressor,
  "RandomForestRegressor" => ENS.RandomForestRegressor,
  "ExtraTreesRegressor" => ENS.ExtraTreesRegressor,
  "AdaBoostRegressor" => ENS.AdaBoostRegressor,
  "GradientBoostingRegressor" => ENS.GradientBoostingRegressor,
  "IsotonicRegression" => ISO.IsotonicRegression,
  "MLPRegressor" => ANN.MLPRegressor
)

mutable struct SKLearner <: TSLearner
    model
    args

    function SKLearner(args=Dict())
        default_args=Dict(
           :output => :class,
           :learner => "LinearSVC",
           :impl_args => Dict()
        )
        new(nothing,mergedict(default_args,args))
    end
end

function fit!(skl::SKLearner, x::T, y::Vector) where {T<:Union{Vector,Matrix}}
  impl_args = copy(skl.args[:impl_args])
  learner = skl.args[:learner]
  py_learner = learner_dict[learner]

  # Assign CombineML-specific defaults if required
  if learner == "RadiusNeighborsClassifier"
    if get(impl_args, :outlier_label, nothing) == nothing
      impl_options[:outlier_label] = labels[rand(1:size(labels, 1))]
    end
  end

  # Train
  skl.model = py_learner(;impl_args...)
  #skl.model[:fit](x, y)
  skl.model.fit(x, y)
end

function transform!(skl::SKLearner, x::T) where {T<:Union{Vector,Matrix}}
  #return collect(skl.model[:predict](x))
  return collect(skl.model.predict(x))
end

function skkrun()
    iris=getiris()
    instances=iris[:,1:4] |> Matrix
    labels=iris[:,5] |> Vector
    model1 = SKLearner(Dict(:learner=>"LinearSVC",:impl_args=>Dict(:max_iter=>5000)))
    model2 = SKLearner(Dict(:learner=>"QDA"))
    model3 = SKLearner(Dict(:learner=>"MLPClassifier"))
    model = SKLearner(Dict(:learner=>"BernoulliNB"))
    fit!(model,instances,labels)
    println(sum(transform!(model,instances).==labels)/length(labels)*100)

    x=iris[:,1:3] |> Matrix
    y=iris[:,4] |> Vector
    #regmodel = SKLearner(Dict(:learner => "SVR",:impl_args=>Dict(:gamma=>"scale")))
    #regmodel = SKLearner(Dict(:learner => "RidgeCV"))
    regmodel = SKLearner(Dict(:learner => "GradientBoostingRegressor"))
    #regmodel = SKLearner(Dict(:learner => "MLPRegressor"))
    fit!(regmodel,x,y)
    println(sum(transform!(regmodel,x).-y)/length(labels)*100)
    xx=iris[:,1] |> Vector
    regmodel = SKLearner(Dict(:learner => "IsotonicRegression"))
    fit!(regmodel,xx,y)
    println(sum(transform!(regmodel,xx).-y)/length(labels)*100)
end

end


text:
function nlfeast!(T, X::AbstractMatrix{ComplexF64}, nodes::Integer, iter::Integer;
    c=complex(0.0,0.0), r=1.0, debug=false, ϵ=10e-12, store=true, spurious=1e-5)

    N, m₀ = size(X)
    Λ = zeros(ComplexF64, m₀)
    θ = LinRange(π/nodes, 2*π-π/nodes, nodes)
    Q₀, Q₁ = zeros(ComplexF64, N, m₀), zeros(ComplexF64, N, m₀)
    Tinv, R = similar(X, ComplexF64), similar(X, ComplexF64)
    resolvent = similar(Λ)
    A, B = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)
    res = Array{Float64}(undef, m₀)
	if store
	    M = Array{ComplexF64}(undef, (nodes, N,N))
	    facts = Array{LU}(undef, nodes)
	end

    qt, rt = qr!(X)
    X .= Matrix(qt)

    for i=1:nodes
        z = (r*exp(θ[i]*im)+c)

		if store
	        M[i,:,:] .= T(z)
	        facts[i] = lu!(M[i,:,:])
	        ldiv!(Tinv, facts[i], X)
	        Tinv .*= (r*exp(θ[i]*im)/nodes)
		else
        	Tinv .= (T(z)\X) .* (r*exp(θ[i]*im)/nodes)
		end

        Q₀ .+= Tinv
        Q₁ .+= Tinv .* z

		if debug print(".") end
    end
	if debug println() end

    beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)

    update_R!(X, R, Λ, T)
    res .= residuals(R, Λ, T)

    if debug
        iter_debug_print(0, Λ, res, c, r, spurious)
    end

    for nit=1:iter

        Q₀ .= 0
        Q₁ .= 0

        for i=1:nodes
            z = (r*exp(θ[i]*im)+c)
            resolvent .= (1 ./(z .- Λ)) .* (r*exp(θ[i]*im)/nodes)

			if store
	            ldiv!(Tinv, facts[i], R)
	            Tinv .= X - Tinv
			else
				Tinv .= (X - T(z)\R)
			end

            rmul!(Tinv,  Diagonal(resolvent))
            Q₀ .+= Tinv
            Q₁ .+= Tinv .* z
			if debug print(".") end
        end
		if debug println() end

		beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)

        update_R!(X, R, Λ, T)
        res .= residuals(R, Λ, T)

        if debug
            iter_debug_print(nit, Λ, res, c, r, spurious)
        end

		res_inside = res[in_contour.(Λ, c, r)]
        if size(res_inside,1) > 0 && maximum(res_inside) < ϵ
            break
        end
		if nit>1 && sum(res_inside .< spurious) > 0 && maximum(res_inside[res_inside .< spurious]) < ϵ
			break
		end
    end

    normalize!(X)
    Λ, X, res
end

function nlfeast_it!(T, X::AbstractMatrix{ComplexF64}, nodes::Integer, iter::Integer;
    c=complex(0.0,0.0), r=1.0, debug=false, ϵ=0.05)

    N, m₀ = size(X)
    Λ = zeros(ComplexF64, m₀)
    θ = LinRange(π/nodes, 2*π-π/nodes, nodes)
    Q₀, Q₁ = zeros(ComplexF64, N, m₀), zeros(ComplexF64, N, m₀)
    R = similar(X, ComplexF64)
    Temp = similar(X, ComplexF64)
    resolvent = similar(Λ)
    Tinv = rand(ComplexF64, nodes, N, m₀)
    A, B = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)
    res = Array{Float64}(undef, size(R, 2))

    for i=1:nodes
        z = (r*exp(θ[i]*im)+c)
        # Temp .= (T(z)\X) .* (r*exp(θ[i]*im)/nodes)
        for j=1:m₀
            # Tinv[i,:,j], log = gmres(T(z), X[:,j]; tol=10e-8, log=true)
            Tinv[i,:,j] .= bicgstabl(T(z), X[:,j], 2; tol=1e-3)
        end
        Temp .= Tinv[i,:,:] .* (r*exp(θ[i]*im)/nodes)
        Q₀ .+= Temp
        Q₁ .+= Temp .* z
		if debug print(".") end
    end
	if debug println() end

    beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
    # beyn_qr_step!(Q₀, Q₁, X, Λ)

    update_R!(X, R, Λ, T)

    if (iter == 0)
        res .= residuals(R, Λ, T)
    end
    if debug
        iter_debug_print(0, Λ, res, c, r)
    end

    for nit=1:iter

        Q₀ .= 0
        Q₁ .= 0

        for i=1:nodes
            z = (r*exp(θ[i]*im)+c)
            resolvent .= (1 ./(z .- Λ))
            # Tinv[i,:,:] .= T(z)\R
            for j=1:m₀
                # Tinv[i,:,j], log = gmres!(Tinv[i,:,j], T(z), R[:,j]; log=true, maxiter=1000, tol=1e-3)
                # Tinv[i,:,j], log = gmres(T(z), R[:,j]; log=true, maxiter=1000, tol=1e-3)
                Tinv[i,:,j], log = bicgstabl!(Tinv[i,:,j], T(z), R[ :,j], 1; log=true, tol=1e-8)
                # Tinv[i,:,j], log = bicgstabl(T(z), R[:,j], 2; log=true, max_mv_products=1000, tol=1e-3)
                # if !log.isconverged println(log) end
                # print(log.iters)
                # print(":\t")
                # println(maximum(log[:resnorm]))
            end
            Temp .= X - Tinv[i,:,:]
            rmul!(Temp, (r*exp(θ[i]*im)/nodes) .* Diagonal(resolvent))
            Q₀ .+= Temp
            Q₁ .+= Temp .* z
			if debug print(".") end
        end
		if debug println() end

        beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
        # beyn_qr_step!(Q₀, Q₁, X, Λ)

        update_R!(X, R, Λ, T)
        res .= residuals(R, Λ, T)

        if debug
            iter_debug_print(nit, Λ, res, c, r)
        end

        if maximum(res[in_contour.(Λ, c, r)]) < ϵ
            break
        end
    end

    normalize!(X)
    Λ, X, residuals(R, Λ, T)
end

function nlfeast_moments!(T, X::AbstractMatrix{ComplexF64}, nodes::Integer, iter::Integer;
    c=complex(0.0,0.0), r=1.0, debug=false, ϵ=10e-12, moments=2, store=true, spurious=1e-5)

    N, m₀ = size(X)
    Λ = zeros(ComplexF64, moments*m₀)
    θ = LinRange(π/nodes, 2*π-π/nodes, nodes)
    Q = zeros(ComplexF64, 2*moments, N, m₀)
    Y = zeros(ComplexF64, N, m₀*moments)
    Q₀, Q₁ = zeros(ComplexF64, moments*N, moments*m₀), zeros(ComplexF64, moments*N, moments*m₀)
    Tinv, Temp, R = similar(X, ComplexF64), similar(X, ComplexF64), similar(Y, ComplexF64)
    resolvent = zeros(ComplexF64, m₀)
    A, B = zeros(ComplexF64, moments*m₀, moments*m₀), zeros(ComplexF64, moments*m₀, moments*m₀)
    res = Array{Float64}(undef, size(R, 2))
    Y = zeros(ComplexF64, N, m₀*moments)
	if store
	    M = Array{ComplexF64}(undef, (nodes, N,N))
	    facts = Array{LU}(undef, nodes)
	end

    for i=1:nodes
        z = (r*exp(θ[i]*im)+c)

		if store
			# M[i,:,:] .= ComplexF32.(T(z))
			M[i,:,:] .= T(z)
	        facts[i] = lu!(M[i,:,:])
	        Tinv .= X
	        ldiv!(facts[i], Tinv)
	        Temp .= (r*exp(θ[i]*im)/nodes) .* Tinv
		else
			Temp .= (T(z)\X) .* (r*exp(θ[i]*im)/nodes)
			# Tinv .= (ComplexF32.(T(z))\ComplexF32.(X)) .* (r*exp(θ[i]*im)/nodes)
		end

        Q[1,:,:] .+= Temp
        for j=2:2*moments
            Q[j,:,:] .+= Temp .* z^(j-1)
        end
    end
    for i=1:moments, j=1:moments
        Q₀[(i-1)*N+1:i*N, (j-1)*m₀+1:j*m₀] .= Q[i+j-1,:,:]
        Q₁[(i-1)*N+1:i*N, (j-1)*m₀+1:j*m₀] .= Q[i+j,:,:]
    end

    S = svd!(Q₀)
    mul!(A, S.U', Q₁)
    mul!(B, A, S.V)
    mul!(A, B, Diagonal(1 ./ S.S))
    F = eigen!(A)
    mul!(Y, S.U[1:N,:], F.vectors)
    Λ .= F.values

    # beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
    # beyn_qr_step!(Q₀, Q₁, X, Λ)

	update_R_moments!(Y, R, Λ, res, T, c, r)
	X .= Y[:,1:m₀]
    # update_R!(Y, R, Λ, T)
    # res .= residuals(R, Λ, T)
    # p = sortperm(res)
	# p = sortperm(abs.(Λ .- c))
    # res .= res[p]
    # X .= Y[:, p][:,1:m₀]
    # Λ .= Λ[p]
    # R .= R[:, p]

    if debug
        iter_debug_print(0, Λ[1:m₀], res[1:m₀], c, r, spurious)
		# iter_debug_print(0, Λ, res, c, r, spurious)
    end

    for nit=1:iter

        Q .= 0

        for i=1:nodes
            z = (r*exp(θ[i]*im)+c)
            resolvent .= (1 ./(z .- Λ[1:m₀])) .* (r*exp(θ[i]*im)/nodes)

			if store
	            Tinv .= R[:, 1:m₀]
	            ldiv!(facts[i], Tinv)
	            Temp .= X - Tinv
			else
				Temp .= (X - T(z)\R[:, 1:m₀])
				# Tinv .= (X - (ComplexF32.(T(z))\ComplexF32.(R[:, 1:m₀])))
			end
            rmul!(Temp, Diagonal(resolvent))

            Q[1,:,:] .+= Temp
            for j=2:2*moments
                Q[j,:,:] .+= Temp .* z^(j-1)
            end
        end

        for i=1:moments, j=1:moments
            Q₀[(i-1)*N+1:i*N, (j-1)*m₀+1:j*m₀] .= Q[i+j-1,:,:]
            Q₁[(i-1)*N+1:i*N, (j-1)*m₀+1:j*m₀] .= Q[i+j,:,:]
        end

        S = svd!(Q₀)
        mul!(A, S.U', Q₁)
        mul!(B, A, S.V)
        mul!(A, B, Diagonal(1 ./ S.S))
        F = eigen!(A)
        mul!(Y, S.U[1:N,:], F.vectors)
        Λ .= F.values

        # beyn_svd_step!(Q₀, Q₁, A, B, X, Λ)
        # beyn_rr_step!(Q₀, Q₁, A, B, X, Λ)
        # beyn_qr_step!(Q₀, Q₁, X, Λ)

        update_R_moments!(Y, R, Λ, res, T, c, r)
		X .= Y[:,1:m₀]
        # res .= residuals(R, Λ, T)
        # # p = sortperm(res)
        # # res .= res[p]
        # # X .= Y[:, p][:,1:m₀]
        # # Λ .= Λ[p]
        # # R .= R[:, p]

        if debug
            iter_debug_print(nit, Λ[1:m₀], res[1:m₀], c, r, spurious)
			# iter_debug_print(nit, Λ, res, c, r, spurious)
        end
		res_inside = res[1:m₀][in_contour.(Λ[1:m₀], c, r)]
        if size(res_inside, 1) > 0 && maximum(res_inside) < ϵ
            break
        end
		if nit>1 && sum(res_inside .< spurious) > 0 && maximum(res_inside[res_inside .< spurious]) < ϵ
			break
		end
    end

	p = sortperm(res)
	res .= res[p]
	X .= Y[:, p][:,1:m₀]
	Λ .= Λ[p]
	R .= R[:, p]

    normalize!(X)
    # Λ[1:m₀], X, res[1:m₀]
	Λ, Y, res
end

ext: jl


text:
module Callbacks

using CUDA
using Dates
using KernelAbstractions
using LinearAlgebra
using MPI
using Printf
using Statistics

using CLIMAParameters
using CLIMAParameters.Planet: day

using ..ClimateMachine
using ..BalanceLaws
using ..Courant
using ..Checkpoint
using ..DGMethods
using ..BalanceLaws: vars_state, Prognostic, Auxiliary
using ..Diagnostics
using ..GenericCallbacks
using ..Mesh.Grids: EveryDirection, HorizontalDirection, VerticalDirection
using ..MPIStateArrays
using ..ODESolvers
using ..TicToc
using ..VariableTemplates
using ..VTK

_sync_device(::Type{CuArray}) = synchronize()
_sync_device(::Type{Array}) = nothing

"""
    SummaryLogCallback([stimeend])

Log a summary of the current run.

The optional `stimeend` argument is used to print the total simulation time.
"""
mutable struct SummaryLogCallback
    stimeend::Union{Nothing, Real}
    wtimestart::DateTime # wall time at start of simulation
    function SummaryLogCallback(stimeend = nothing)
        new(stimeend, now())
    end
end

function GenericCallbacks.init!(cb::SummaryLogCallback, solver, Q, param, t)
    cb.wtimestart = now()
    return nothing
end
function GenericCallbacks.call!(cb::SummaryLogCallback, solver, Q, param, t)
    wallclock_time_ms = Dates.now() - cb.wtimestart
    wallclock = Dates.format(
        convert(Dates.DateTime, wallclock_time_ms),
        Dates.dateformat"HH:MM:SS",
    )
    normQ = norm(Q)
    if cb.stimeend === nothing
        simtime = @sprintf "%8.2f" t
    else
        simtime = @sprintf "%8.2f / %8.2f" t cb.stimeend
    end
    wallclock_s = wallclock_time_ms.value / 1000
    efficiency = @sprintf "%8.4f" t / wallclock_s

    estimated_wallclock_end = cb.stimeend * wallclock_s / t

    estimated_wallclock_end_ms =
        Dates.Millisecond(ceil(Int, estimated_wallclock_end * 1000))

    estimated_remaining = Dates.format(
        convert(Dates.DateTime, estimated_wallclock_end_ms),
        Dates.dateformat"HH:MM:SS",
    )
    @info @sprintf(
        """
        Update
            simtime = %s
            wallclock = %s
            efficiency (simtime / wallclock) = %s
            wallclock end (estimated) = %s
            norm(Q) = %.16e""",
        simtime,
        wallclock,
        efficiency,
        estimated_remaining,
        normQ,
    )
    if !isfinite(normQ)
        show_not_finite_fields(Q)
        error("norm(Q) is not finite")
    end
    return nothing
end
function GenericCallbacks.fini!(cb::SummaryLogCallback, solver, Q, param, t)
    return nothing
end

"""
    show_updates(show_updates_opt, solver_config, user_info_callback)

Return a callback function that shows simulation updates at the specified
interval and also invokes the user-specified info callback.
"""
function show_updates(show_updates_opt, solver_config, user_info_callback)
    timeend = solver_config.timeend

    cb_constr = CB_constructor(show_updates_opt, solver_config)
    if cb_constr !== nothing
        return cb_constr((
            SummaryLogCallback(solver_config.timeend),
            GenericCallbacks.AtInit(user_info_callback),
        ))
    else
        return nothing
    end
end

"""
    diagnostics(diagnostics_opt, solver_config, dgn_starttime, diagnostics_config)

Return callback functions that shows simulation updates at the specified
interval and also invokes the user-specified info callback.
"""
function diagnostics(
    diagnostics_opt,
    solver_config,
    dgn_starttime,
    diagnostics_config,
)
    if diagnostics_opt != "never" && diagnostics_config !== nothing
        # set up a callback for each diagnostics group
        dgncbs = ()
        for dgngrp in diagnostics_config.groups
            cb_constr =
                CB_constructor(diagnostics_opt, solver_config, dgngrp.interval)
            cb_constr === nothing && continue
            dgncbs = (dgncbs..., cb_constr(dgngrp))
        end
        return dgncbs
    else
        return nothing
    end
end

"""
    vtk(vtk_opt, solver_config)

Return a callback that saves state and auxiliary variables to a VTK
file.
"""
function vtk(vtk_opt, solver_config, output_dir, number_sample_points)
    cb_constr = CB_constructor(vtk_opt, solver_config)
    cb_constr === nothing && return nothing

    vtknum = Ref(1)

    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    bl = dg.balance_law
    Q = solver_config.Q
    FT = eltype(Q)

    cb_vtk = GenericCallbacks.AtInitAndFini() do
        # TODO: make an object
        vprefix = @sprintf(
            "%s_mpirank%04d_num%04d",
            solver_config.name,
            MPI.Comm_rank(mpicomm),
            vtknum[],
        )
        outprefix = joinpath(output_dir, vprefix)

        statenames = flattenednames(vars_state(bl, Prognostic(), FT))
        auxnames = flattenednames(vars_state(bl, Auxiliary(), FT))

        writevtk(
            outprefix,
            Q,
            dg,
            statenames,
            dg.state_auxiliary,
            auxnames;
            number_sample_points = number_sample_points,
        )

        # Generate the pvtu file for these vtk files
        if MPI.Comm_rank(mpicomm) == 0
            # name of the pvtu file
            pprefix = @sprintf("%s_num%04d", solver_config.name, vtknum[])
            pvtuprefix = joinpath(output_dir, pprefix)

            # name of each of the ranks vtk files
            prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
                @sprintf(
                    "%s_mpirank%04d_num%04d",
                    solver_config.name,
                    i - 1,
                    vtknum[],
                )
            end
            writepvtu(
                pvtuprefix,
                prefixes,
                (statenames..., auxnames...),
                eltype(Q),
            )
        end

        vtknum[] += 1
        nothing
    end
    return cb_constr(cb_vtk)
end

"""
    monitor_timestep_duration(mtd_opt, array_type, comm)

Returns a callback function that displays wall-clock time per time-step
statistics across MPI ranks in the communicator `comm`.  The times are
averaged over the `mtd_opt` time steps requested for output.  The
`array_type` is used to synchronize the compute device.
"""
function monitor_timestep_duration(mtd_opt, array_type, comm)
    if mtd_opt == "never"
        return nothing
    elseif !endswith(mtd_opt, "steps")
        @warn @sprintf(
            """
            monitor-timestep-duration must be in 'steps'; %s unrecognized; disabling""",
            mtd_opt,
        )
        return nothing
    end
    steps = parse(Int, mtd_opt[1:(end - 5)])

    _sync_device(array_type)
    before = time_ns()

    cb_mtd = GenericCallbacks.EveryXSimulationSteps(steps) do
        _sync_device(array_type)
        after = time_ns()

        time_per_timesteps = after - before

        times = MPI.Gather(time_per_timesteps, 0, comm)
        if MPI.Comm_rank(comm) == 0
            ns_per_s = 1e9
            times = times ./ ns_per_s ./ steps

            @info @sprintf(
                """
                Wall-clock time per time-step (statistics across MPI ranks)
                   maximum (s) = %25.16e
                   minimum (s) = %25.16e
                   median  (s) = %25.16e
                   std     (s) = %25.16e""",
                maximum(times),
                minimum(times),
                median(times),
                std(times),
            )
        end

        _sync_device(array_type)
        before = time_ns()

        nothing
    end

    return cb_mtd
end

"""
    monitor_courant_numbers(mcn_opt, solver_config)

Return a callback function that displays Courant numbers for the simulation
at `mcn_opt` intervals.
"""
function monitor_courant_numbers(mcn_opt, solver_config)
    cb_constr = CB_constructor(mcn_opt, solver_config)
    cb_constr === nothing && return nothing

    cb_cfl = cb_constr() do
        Δt = solver_config.dt
        c_v = courant(
            nondiffusive_courant,
            solver_config,
            direction = VerticalDirection(),
        )
        c_h = courant(
            nondiffusive_courant,
            solver_config,
            direction = HorizontalDirection(),
        )
        ca_v = courant(
            advective_courant,
            solver_config,
            direction = VerticalDirection(),
        )
        ca_h = courant(
            advective_courant,
            solver_config,
            direction = HorizontalDirection(),
        )
        cd_v = courant(
            diffusive_courant,
            solver_config,
            direction = VerticalDirection(),
        )
        cd_h = courant(
            diffusive_courant,
            solver_config,
            direction = HorizontalDirection(),
        )
        simtime = ODESolvers.gettime(solver_config.solver)
        @info @sprintf(
            """
            Courant numbers at simtime: %8.2f, Δt = %8.2f s
                Acoustic (vertical) Courant number    = %.2g
                Acoustic (horizontal) Courant number  = %.2g
                Advection (vertical) Courant number   = %.2g
                Advection (horizontal) Courant number = %.2g
                Diffusion (vertical) Courant number   = %.2g
                Diffusion (horizontal) Courant number = %.2g""",
            simtime,
            Δt,
            c_v,
            c_h,
            ca_v,
            ca_h,
            cd_v,
            cd_h,
        )
        nothing
    end
    return cb_cfl
end

function adapt_timestep(adp_opt, solver_config)
    cb_constr = CB_constructor(adp_opt, solver_config)
    cb_constr === nothing && return nothing
    cb_adp = cb_constr() do
        dt = solver_config.solver.dt
        dg = solver_config.dg
        bl = dg.balance_law
        Q = solver_config.Q
        t0 = solver_config.t0
        ode_solver_type = solver_config.ode_solver_type
        dtmodel = getdtmodel(ode_solver_type, bl)
        ndt = ClimateMachine.DGMethods.calculate_dt(
            dg,
            dtmodel,
            Q,
            solver_config.CFL,
            t0,
            solver_config.diffdir,
        )
        @info @sprintf("""Updating time step: %8.16f => %8.16f""", dt, ndt)
        updatedt!(solver_config.solver, ndt)
        nothing
    end
    return cb_adp
end

"""
    checkpoint(
        checkpoint_opt,
        checkpoint_keep_one,
        solver_config,
        checkpoint_dir,
    )

Return a callback function that runs at `checkpoint_opt` intervals
and stores a simulation checkpoint into `checkpoint_dir` identified
by a running number.
"""
function checkpoint(
    checkpoint_opt,
    checkpoint_keep_one,
    solver_config,
    checkpoint_dir,
)
    cb_constr = CB_constructor(checkpoint_opt, solver_config)
    cb_constr === nothing && return nothing

    cpnum = Ref(1)
    cb_checkpoint = cb_constr() do
        write_checkpoint(
            solver_config,
            checkpoint_dir,
            solver_config.name,
            solver_config.mpicomm,
            cpnum[],
        )
        if checkpoint_keep_one
            rm_checkpoint(
                checkpoint_dir,
                solver_config.name,
                solver_config.mpicomm,
                cpnum[] - 1,
            )
        end
        cpnum[] += 1
        nothing
    end
    return cb_checkpoint
end

"""
    ConsCallback(dg, mass, energy, show_cons)

Check mass and energy conservation against specified tolerances.
"""
mutable struct ConsCallback{FT}
    bl::BalanceLaw
    varname::String
    error_threshold::FT
    show::Bool
    Σvar₀::FT
end

function GenericCallbacks.init!(cb::ConsCallback, solver, Q, param, t)
    FT = eltype(Q)
    idx = varsindices(vars_state(cb.bl, Prognostic(), FT), cb.varname)
    cb.Σvar₀ = weightedsum(Q, idx)
    return nothing
end
function GenericCallbacks.call!(cb::ConsCallback, solver, Q, param, t)
    FT = eltype(Q)
    idx = varsindices(vars_state(cb.bl, Prognostic(), FT), cb.varname)
    Σvar = weightedsum(Q, idx)
    δvar = (Σvar - cb.Σvar₀) / cb.Σvar₀

    if abs(δvar) > cb.error_threshold
        error("abs(δ$(cb.varname)) > $(cb.error_threshold)")
    end

    if cb.show
        simtime = @sprintf "%8.2f" t
        @info @sprintf(
            """
            Conservation
                simtime = %s
                abs(δ%s) = %.5e""",
            simtime,
            cb.varname,
            abs(δvar),
        )
    end
    return nothing
end
function GenericCallbacks.fini!(cb::ConsCallback, solver, Q, param, t)
    return nothing
end

"""
    check_cons(
        check_cons,
        solver_config,
    )

Return a callback function for each element of `check_cons`, which must be
a tuple of [`ConservationCheck`s](@ref ClimateMachine.ConservationCheck).
"""
function check_cons(check_cons, solver_config)
    cbs = ()
    for cc in check_cons
        cb_constr = CB_constructor(cc.interval, solver_config)
        cb_constr === nothing && continue

        FT = eltype(solver_config.Q)
        cb = cb_constr((ConsCallback(
            solver_config.dg.balance_law,
            cc.varname,
            FT(cc.error_threshold),
            cc.show,
            FT(0),
        ),))

        cbs = (cbs..., cb)
    end

    return cbs
end

"""
    CB_constructor(interval, solver_config, default)

Parse the specified `interval` and return the appropriate `GenericCallbacks`
constructor. If `interval` is "default", then use `default` as the interval.
"""
function CB_constructor(interval::String, solver_config, default = nothing)
    mpicomm = solver_config.mpicomm
    solver = solver_config.solver
    dg = solver_config.dg
    bl = dg.balance_law
    param_set = parameter_set(bl)
    secs_per_day = day(param_set)

    if (
        m = match(r"^([0-9\.]+)(smonths|sdays|shours|smins|ssecs)$", interval)
    ) !== nothing
        n = parse(Float64, m[1])
        if m[2] == "smonths"
            ssecs = 30 * secs_per_day * n
        elseif m[2] == "sdays"
            ssecs = secs_per_day * n
        elseif m[2] == "shours"
            ssecs = 60 * 60 * n
        elseif m[2] == "smins"
            ssecs = 60 * n
        elseif m[2] == "ssecs"
            ssecs = n
        end
        return cb -> GenericCallbacks.EveryXSimulationTime(cb, ssecs)
    elseif (m = match(r"^([0-9]+)(steps)$", interval)) !== nothing
        steps = parse(Int, m[1])
        return cb -> GenericCallbacks.EveryXSimulationSteps(cb, steps)
    elseif (m = match(r"^([0-9\.]+)(hours|mins|secs)$", interval)) !== nothing
        n = parse(Float64, m[1])
        if m[2] == "hours"
            secs = 60 * 60 * n
        elseif m[2] == "mins"
            secs = 60 * n
        elseif m[2] == "secs"
            secs = n
        end
        return cb -> GenericCallbacks.EveryXWallTimeSeconds(cb, secs, mpicomm)
    elseif interval == "default"
        if default === nothing
            @warn "no default available; ignoring"
            return nothing
        elseif default === ""
            return nothing
        else
            return CB_constructor(default, solver_config, "")
        end
    elseif interval == "never"
        return nothing
    else
        @warn @sprintf(
            """
%s: unrecognized interval; ignoring""",
            interval,
        )
        return nothing
    end
end

function __init__()
    tictoc()
end

end # module


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 1.0.5
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# ### MRI reconstruction introductory notebook for ISBI 2019 tutorial
#
# This jupyter notebook is in a directory with other notebooks
# that demonstrate MRI reconstruction using the Julia language
# and the Michigan Image Reconstruction Toolbox (MIRT).
#
# 2019-03-14 Jeff Fessler, University of Michigan

# ### Initial steps
#
# See getting started instructions here:
# https://github.com/JeffFessler/MIRT.jl/blob/master/doc/start.md

# ### Test your installation
#
# The next cell loads the packages needed for this initial test.
# It might take a little bit of time to run the first time you try it
# because Julia is essentially a compiled language under the hood,
# even though it feels like an interactive language to the user,
# so it will be compiling things behind the scenes.

# Load all packages needed for this demo
using MIRT: ellipse_im, jim # https://github.com/JeffFessler/MIRT.jl
using FFTW: fft

# ### View the Shepp-Logan image and its spectrum
# The `ellipse_im` function in MIRT can generate phantom images consisting of ellipses such as the Shepp-Logan image.
#
# The `jim` function in MIRT is a jiffy image display routine.

image = ellipse_im(256, oversample=2)
jim(image, "Shepp-Logan Phantom")

# calculate k-space data using fft
kspace = fft(image) # this is a 2D FFT because "image" is a 2D array
jim(kspace, fft0=true) # show k-space with 0 at center, ala matlab's "fftshift"

# It is hard to see much on a linear scale, so now use a log scale
logger = (x ; min=-6) -> log10.(max.(abs.(x) / maximum(abs.(x)), (10.)^min))
jim(logger(kspace), "k-space data on log scale", fft0=true)

# ### Getting help
# In addition to the documentation linked above,
# to learn about any command in Julia just type `?` followed by the command name
# and Julia will return the documentation string
# that goes with that command.
# Here is an example.
# This example illustrates the very useful "multiple dispatch" feature of Julia;
# essentially there are multiple versions of this function
# that you can call in several different ways
# with different argument combinations,
# and with optimal keyword arguments.

#?jim # uncomment this line to see

# ### End of introduction!
#
# If you got this far you are ready for more interesting demos next.

max_stars_repo_name: cecroc/mirt-demo

max_issues_repo_name: cecroc/mirt-demo

max_forks_repo_name: cecroc/mirt-demo


text:
# @time using Dates
@time using CSV, XLSX, DataFrames
excuted_DIR = string(@__DIR__)
schedule = DataFrame(XLSX.readtable(excuted_DIR * "/schedule.xlsx", "schedule")...)

@time using Crayons
@time using Random, Distributions, Statistics
@time using Graphs, NearestNeighbors

@time include("src/lemma.jl")
@time include("src/main.jl")

# ------------------------------------------------------------------ directory

root = "D:/simulated/"
if !isdir(root)
    mkpath(root)
end
cd(root); pwd()

# ------------------------------------------------------------------ setting

export_type = :CSV # :both, :CSV, :XLSX
doing = parse(Int64, ARGS[1])

scenario = schedule[doing,:]
try
    mkpath(root * scenario.name)
    mkpath(root * scenario.name * "/00")
    mkpath(root * scenario.name * "/0V")
    mkpath(root * scenario.name * "/M0")
    mkpath(root * scenario.name * "/MV")
    CSV.write(root * scenario.name * "/cnfg.csv", DataFrame(scenario), bom = true)
    cd(root * scenario.name)
catch
    @warn "some error occur!"
end

# ------------------------------------------------------------------ parameters

global number_of_host = 1
global θ = 10

global network = scenario.network
global σ = scenario.σ
global e_σ = scenario.e_σ
global p_V = scenario.p_V
global e_V = scenario.e_V
global β = scenario.β
global n = scenario.n
global δ = scenario.δ
global end_time = scenario.end_time
global first_seed = scenario.first_seed
global last_seed = scenario.last_seed

global ID = 1:n
global brownian = MvNormal(2, 0.01) # moving process

if network == 0
    global N = scenario.N
    global m = 3 # number of network link
    backbone = barabasi_albert(N, m, seed = 0)
    global NODE = Dict(1:N .=> backbone.fadjlist)
    CSV.write(root * "ntwk.csv", DataFrame(
        deg = degree(backbone),
        betweenness = betweenness_centrality(backbone),
        closeness = closeness_centrality(backbone),
        stress = stress_centrality(backbone)
        ))
elseif network == -1
    using JLD
    realnetwork = load(excuted_DIR * "/world_wide_flight_network.jld")
    global NODE = Dict(realnetwork["nodes"] .=> realnetwork["adj"])
    global N = length(NODE)
end
global NODE_ID = keys(NODE)

global latent_period = Weibull(3, 7.17) # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7014672/#__sec2title
global recovery_period = Weibull(3, 7.17)

@time for v ∈ [true, false], m ∈ [true, false]
    for seed_number ∈ first_seed:last_seed
        simulation(
        seed_number,
        moving = m,
        vaccin = v)
    end
end

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
# ---
# title: 994. Rotting Oranges
# id: problem994
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/rotting-oranges/description/>
# hidden: true
# ---
#
# In a given grid, each cell can have one of three values:
#
#   * the value `0` representing an empty cell;
#   * the value `1` representing a fresh orange;
#   * the value `2` representing a rotten orange.
#
# Every minute, any fresh orange that is adjacent (4-directionally) to a rotten
# orange becomes rotten.
#
# Return the minimum number of minutes that must elapse until no cell has a
# fresh orange.  If this is impossible, return `-1` instead.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/02/16/oranges.png)**
#
#
#
#     Input: [[2,1,1],[1,1,0],[0,1,1]]
#     Output: 4
#
#
# **Example 2:**
#
#
#
#     Input: [[2,1,1],[0,1,1],[1,0,1]]
#     Output: -1
#     Explanation: The orange in the bottom left corner (row 2, column 0) is never rotten, because rotting only happens 4-directionally.
#
#
# **Example 3:**
#
#
#
#     Input: [[0,2]]
#     Output: 0
#     Explanation: Since there are already no fresh oranges at minute 0, the answer is just 0.
#
#
#
#
# **Note:**
#
#   1. `1 <= grid.length <= 10`
#   2. `1 <= grid[0].length <= 10`
#   3. `grid[i][j]` is only `0`, `1`, or `2`.
#
#
## @lc code=start
using LeetCode

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


text:
using PyPlot, BinningAnalysis

function make_plot(xs, label)
    figure()
    title(label)
    plot(eachindex(xs), xs)
    xlabel("Binnging level")
    ylabel(label)
end

# Create Binning Analysis
BA = LogBinner()

# Get some data
N_corr = 16
N_blocks = 131_072
for _ in 1:N_blocks
    x = rand()
    for __ in 1:N_corr
        push!(BA, x)
    end
end

means = all_means(BA)
vars = all_vars(BA)
varNs = all_varNs(BA)
errors = all_std_errors(BA)
taus = all_taus(BA)

# Means should be constant across all levels
make_plot(means, "mean")

# The variance decays with increasing lvl, because the sample size decreases.
make_plot(vars, "Var")

# varN is the Variance normalized to the sample size. It increases for correlated
# measurements and stays roughly constant for uncorrelated ones. It fluctuates
# more with decreasing sample size (i.e. increasing binning level)
make_plot(varNs, "Var/N")

# The standard error gives the error of the mean. It behaves like varN.
make_plot(errors, "Standard Error")

# The auto correlation time should be 16... ?
make_plot(taus, "Autocorrelation Time")

# Convergence:
# Green dots are accepted as converged, red dots are not. In general convergence
# should be accepted when the Var/N becomes roughly constant
make_plot(varNs, "Var/N")
for lvl in 1:20
    tag = if has_converged(BA, lvl)
        "go"
    else
        "ro"
    end
    plot([lvl], [varNs[lvl]], tag)
end


text:
radius2(x, y) = x^2 + y^2
function loop_radius2(n)
    s = 0
    for i = 1:n
        s += radius2(1, i)
    end
    s
end

@testset "Breakpoints" begin
    breakpoint(radius2)
    stack = JuliaStackFrame[]
    frame = JuliaInterpreter.enter_call(loop_radius2, 2)
    bp = JuliaInterpreter.finish_and_return!(stack, frame)
    @test isa(bp, Breakpoints.BreakpointRef)
    @test length(stack) == 2
    @test stack[end].code.scope == @which radius2(0, 0)
    bp = JuliaInterpreter.finish_stack!(stack)
    @test isa(bp, Breakpoints.BreakpointRef)
    @test length(stack) == 2
    @test JuliaInterpreter.finish_stack!(stack) == loop_radius2(2)

    # Conditional breakpoints
    function runsimple()
        stack = JuliaStackFrame[]
        frame = JuliaInterpreter.enter_call(loop_radius2, 2)
        bp = JuliaInterpreter.finish_and_return!(stack, frame)
        @test isa(bp, Breakpoints.BreakpointRef)
        @test length(stack) == 2
        @test stack[end].code.scope == @which radius2(0, 0)
        @test JuliaInterpreter.finish_stack!(stack) == loop_radius2(2)
    end
    remove()
    breakpoint(radius2, :(y > x))
    runsimple()
    remove()
    @breakpoint radius2(0,0) y>x
    runsimple()
    # Demonstrate the problem that we have with scope
    local_identity(x) = identity(x)
    remove()
    @breakpoint radius2(0,0) y>local_identity(x)
    @test_broken @interpret loop_radius2(2)

    # Conditional breakpoints on local variables
    remove()
    halfthresh = loop_radius2(5)
    @breakpoint loop_radius2(10) 5 s>$halfthresh
    stack, bp = @interpret loop_radius2(10)
    @test isa(bp, Breakpoints.BreakpointRef)
    frame = stack[end]
    s_extractor = eval(Breakpoints.prepare_slotfunction(frame.code, :s))
    @test s_extractor(frame) == loop_radius2(6)
    JuliaInterpreter.finish_stack!(stack)
    @test s_extractor(frame) == loop_radius2(7)
    disable(bp)
    @test JuliaInterpreter.finish_stack!(stack) == loop_radius2(10)

    # Next line with breakpoints
    function outer(x)
        inner(x)
    end
    function inner(x)
        return 2
    end
    breakpoint(inner)
    stack = JuliaStackFrame[]
    frame = JuliaInterpreter.enter_call(outer, 2)
    bp = JuliaInterpreter.next_line!(stack, frame)
    @test isa(bp, Breakpoints.BreakpointRef)
    @test JuliaInterpreter.finish_stack!(stack) == 2

    # Direct return
    @breakpoint gcd(1,1) a==5
    @test @interpret(gcd(10,20)) == 10
    # FIXME: even though they pass, these tests break Test!
    # stack, bp = @interpret gcd(5, 20)
    # @test length(stack) == 1 && isa(stack[1], JuliaStackFrame)
    # @test isa(bp, Breakpoints.BreakpointRef)
    remove()

    # break on error
    try
        JuliaInterpreter.break_on_error[] = true

        inner(x) = error("oops")
        outer() = inner(1)
        stack = JuliaStackFrame[]
        frame = JuliaInterpreter.enter_call(outer)
        bp = JuliaInterpreter.finish_and_return!(stack, frame)
        @test bp.err == ErrorException("oops")
        @test length(stack) >= 2
        @test stack[1].code.scope.name == :outer
        @test stack[2].code.scope.name == :inner

        f_catch() = try error(); catch; return 2; end
        @test @interpret f_catch() == 2

    finally
        JuliaInterpreter.break_on_error[] = false
    end

    # Breakpoint display
    io = IOBuffer()
    frame = JuliaInterpreter.enter_call(loop_radius2, 2)
    bp = JuliaInterpreter.BreakpointRef(frame.code, 1)
    show(io, bp)
    @test String(take!(io)) == "breakpoint(loop_radius2(n) in $(@__MODULE__) at $(@__FILE__):3, 3)"
    bp = JuliaInterpreter.BreakpointRef(frame.code, 0)  # fictive breakpoint
    show(io, bp)
    @test String(take!(io)) == "breakpoint(loop_radius2(n) in $(@__MODULE__) at $(@__FILE__):3, %0)"
    bp = JuliaInterpreter.BreakpointRef(frame.code, 1, ArgumentError("whoops"))
    show(io, bp)
    @test String(take!(io)) == "breakpoint(loop_radius2(n) in $(@__MODULE__) at $(@__FILE__):3, 3, ArgumentError(\"whoops\"))"
end