text:
using Mosek
export solve_mosek

printstream(msg::String) = print(msg)

function get_SDPtriplets(problem::SDP_Problem; debug = false)

  nzc = 0
  nza = 0
  for ((objctr, block, var1, var2), coeff) in problem.matrices

      if objctr in problem.obj_keys
        nzc += 1
      else
        nza += 1
      end
  end
  # println("nza : ", nza)
  # println("nzc : ", nzc)
  barai = Int32[0 for i in 1:nza]
  baraj = Int32[0 for i in 1:nza]
  barak = Int32[0 for i in 1:nza]
  baral = Int32[0 for i in 1:nza]
  baraijkl = Float64[0 for i in 1:nza]

  barcj = Int32[0 for i in 1:nzc]
  barck = Int32[0 for i in 1:nzc]
  barcl = Int32[0 for i in 1:nzc]
  barcjkl = Float64[0 for i in 1:nzc]

  nzc=0
  nza=0
  for ((objctr, block, var1, var2), coeff) in problem.matrices
    sdp_block = problem.name_to_sdpblock[block]
    lower = min(sdp_block.var_to_id[var1], sdp_block.var_to_id[var2])
    upper = max(sdp_block.var_to_id[var1], sdp_block.var_to_id[var2])
    if objctr in problem.obj_keys
      nzc+=1
      barcj[nzc] = sdp_block.id
      barck[nzc] = upper
      barcl[nzc] = lower
      barcjkl[nzc] = coeff
    else
      nza+=1
      barai[nza] = problem.name_to_ctr[objctr][1]
      baraj[nza] = sdp_block.id
      barak[nza] = upper
      baral[nza] = lower
      baraijkl[nza] = coeff
    end
  end

  if debug
    println("*********************************************************************************")
    println("Debug -> Reading SDP_Problem")
    @printf("%5s  %5s  %5s  %s\n", "barcj", "barck", "barcl", "barcjkl")
    for i=1:length(barcj)
        @printf("%5i  %5i  %5i  %f\n", barcj[i], barck[i], barcl[i], barcjkl[i])
    end

    @printf("%5s  %5s  %5s  %5s  %s\n", "barai", "baraj", "barak", "baral", "baraijkl")
    for i=1:length(barai)
        @printf("%5i  %5i  %5i  %5i  %f\n", barai[i], baraj[i], barak[i], baral[i], baraijkl[i])
    end
  end
  return barcj, barck, barcl, barcjkl, barai, baraj, barak, baral, baraijkl
end

function get_linterms(problem; debug=debug)
  nza, nzc = 0, 0
  for (objctr, var) in keys(problem.linear)
    if objctr in problem.obj_keys
      nzc += 1
    else
      nza += 1
    end
  end

  ai = zeros(Int32, nza)
  aj = zeros(Int32, nza)
  aij = zeros(Float64, nza)
  cj = zeros(Int32, nzc)
  cjval = zeros(Float64, nzc)

  nza, nzc = 0, 0
  for ((objctr, var), coeff) in problem.linear
    if objctr in problem.obj_keys
      nzc += 1
      cj[nzc] = problem.scalvar_to_id[var]
      cjval[nzc] = coeff
    else
      nza += 1
      ai[nza] = problem.name_to_ctr[objctr][1]
      aj[nza] = problem.scalvar_to_id[var]
      aij[nza] = coeff
    end
  end

  if debug
    warn(LOGGER, "--- cj :")
    for i=1:length(cj)
      @printf("%i  %f\n", cj[i], cjval[i])
    end
    warn(LOGGER, "--- a :")
    for i=1:length(aj)
      @printf("%i  %i  %f\n", ai[i], aj[i], aij[i])
    end
  end
  return cj, cjval, ai, aj, aij
end

function get_ctrbounds(problem::SDP_Problem; debug = false)
  numcon=length(problem.name_to_ctr)
  MOSEK_KIND = Dict(["EQ"=>MSK_BK_FX, "GEQ"=>MSK_BK_LO, "LEQ"=>MSK_BK_UP, "RNG"=>MSK_BK_RA])
  bkc = Boundkey[ Mosek.Boundkey(1)  for kv in problem.name_to_ctr]
  # constraint
  buc = Float64[0 for i in 1:numcon]
  blc = Float64[0 for i in 1:numcon]
  for (ctrname, ctr) in problem.name_to_ctr
    # @printf("%10d%20s%10d%20s\n", ctr[1], ctrname, MOSEK_KIND[ctr[2]], ctr[2])
    id_ctr=ctr[1]
    lb = ctr[3]
    ub = ctr[4]
    cst = 0
    if haskey(problem.cst_ctr, ctrname)
      cst = problem.cst_ctr[ctrname]
    end
    bkc[id_ctr] = MOSEK_KIND[ctr[2]]
    if bkc[id_ctr] == MSK_BK_UP
      buc[id_ctr] = ub - cst
    elseif bkc[id_ctr] == MSK_BK_LO
      blc[id_ctr] = lb - cst
    elseif bkc[id_ctr] == MSK_BK_FX
      blc[id_ctr] = lb - cst
      buc[id_ctr] = ub - cst
    elseif bkc[id_ctr] == MSK_BK_RA
      blc[id_ctr] = lb - cst
      buc[id_ctr] = ub - cst
    else
      error(LOGGER, "get_ctrbounds() : Unknown constraint kind $(ctr[2]) $(bkc[id_ctr]) $(MSK_BK_FX[1])")
    end
  end
  if debug
    warn("get_ctrbounds(): done")
    @show numcon
    @show bkc
    @show blc
    @show buc
  end
  return numcon, bkc, blc, buc
end

function get_varbounds(problem::SDP_Problem)
  MSK_INFINITY = 1.0e30

  numvar = length(problem.scalvar_to_id)

  sub = [i for i in 1:numvar]
  bkx = [MSK_BK_FR for i in 1:numvar]
  blx = [-MSK_INFINITY for i in 1:numvar]
  bux = [MSK_INFINITY for i in 1:numvar]

  return sub, bkx, blx, bux
end

function solve_mosek(problem::SDP_Problem, primal::SortedDict{Tuple{String,String,String}, Float64},
                                           dual::SortedDict{Tuple{String, String, String}, Float64};
                                           debug = false,
                                           logname = "",
                                           printlog = true,
                                           msk_maxtime = -1,            #Default -1 means no time limit
                                           sol_info = OrderedDict())
  empty!(primal)
  empty!(dual)
  primobj = NaN
  dualobj = NaN

  # Get nb of variables, dimensions, nb of ctrs, bounds.
  nbarvar = length(problem.id_to_sdpblock)
  barvardim = [ length(problem.id_to_sdpblock[block].var_to_id) for block in 1:nbarvar ]

  numvar = length(problem.scalvar_to_id)

  numcon, bkc, blc, buc = get_ctrbounds(problem)
  sub, bkx, blx, bux = get_varbounds(problem)

  # Get objective and ctr matrix contribution
  barcj, barck, barcl, barcjkl, barai, baraj, barak, baral, baraijkl = get_SDPtriplets(problem, debug=debug)

  # Get obj and ctr linear contribution
  cj, cjval, ai, aj, aij = get_linterms(problem, debug=debug)

  # Create a task object and attach log stream printer
  maketask() do task
      logname != "" && linkfiletostream(task, MSK_STREAM_LOG, logname, 0)

      printlog && putstreamfunc(task,MSK_STREAM_LOG,printstream)

      # Append SDP matrix variables and scalar variables.
      # The variables will initially be fixed at zero.
      appendbarvars(task,barvardim)
      appendvars(task, numvar)

      putvarboundlist(task, sub, bkx, blx, bux)

      # Append 'numcon' empty constraints.
      # The constraints will initially have no bounds.
      appendcons(task,numcon)
      # Set the bounds on constraints.
      putconboundslice(task,1,numcon+1, bkc,blc,buc)

      # Minimize
      putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)

      # Set constraints SDP vars coeffs
      putbarablocktriplet(task, length(barai), barai, baraj, barak, baral, baraijkl)

      # Set contraints Sym vars coeffs as linear part
      putaijlist(task, ai, aj, aij)

      # Objective matrices and constant
      putbarcblocktriplet(task, length(barcj), barcj, barck, barcl, barcjkl)
      putclist(task, cj, cjval)

      obj_cst = 0
      for objmoment in problem.obj_keys
        if haskey(problem.cst_ctr, objmoment)
          obj_cst += problem.cst_ctr[objmoment]
        end
      end
      obj_cst != 0 && putcfix(task, obj_cst)

      # putintparam(task, MSK_IPAR_INTPNT_SCALING, MSK_SCALING_NONE)
      # putintparam(task, MSK_IPAR_INTPNT_SCALING, MSK_SCALING_AGGRESSIVE)

      # putdouparam(task, MSK_DPAR_INTPNT_CO_TOL_DFEAS, 1e-12)
      # putdouparam(task, MSK_DPAR_INTPNT_CO_TOL_INFEAS, 1e-12)
      # putdouparam(task, MSK_DPAR_INTPNT_CO_TOL_MU_RED, 1e-12)
      # putdouparam(task, MSK_DPAR_INTPNT_CO_TOL_PFEAS, 1e-12)
      # putdouparam(task, MSK_DPAR_INTPNT_CO_TOL_REL_GAP, 1e-1)

      # Set Mosek maxtime, default is infinity
      putdouparam(task, MSK_DPAR_OPTIMIZER_MAX_TIME, msk_maxtime)

      if debug
        dump_mosek_model(task)
      end

      optimize(task)
      solutionsummary(task,MSK_STREAM_MSG)

      # Get status information about the solution
      prosta = getprosta(task, MSK_SOL_ITR)
      solsta = getsolsta(task, MSK_SOL_ITR)

      save_SDPsolinfo!(task, sol_info)

      if solsta == MSK_SOL_STA_OPTIMAL || solsta == MSK_SOL_STA_NEAR_OPTIMAL
        # Get objective
        primobj = getprimalobj(task, MSK_SOL_ITR)
        dualobj = getdualobj(task, MSK_SOL_ITR)

        # Get primal and dual solutions
        get_primalsol!(task, problem, primal)
        get_dualsol!(task, problem, dual)

      elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
          println("Dual infeasibility.\n")
      elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
          println("Primal infeasibility.\n")
      elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
          println("Near dual infeasibility.\n")
      elseif  solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
          println("Near primal infeasibility.\n")
      elseif  solsta == MSK_SOL_STA_UNKNOWN
          println("Unknown solution status.\n")
      else
          println("Other solution status.\n")
      end
  end

  return primobj, dualobj
end

function save_SDPsolinfo!(task::Mosek.Task, sol_info::OrderedDict)
  prosta = getprosta(task, MSK_SOL_ITR)
  solsta = getsolsta(task, MSK_SOL_ITR)

  sol_info[:slv_prosta] = get_prosta_map()[prosta]
  sol_info[:slv_solsta] = get_solsta_map()[solsta]

  if solsta == MSK_SOL_STA_OPTIMAL || solsta == MSK_SOL_STA_NEAR_OPTIMAL
    sol_info[:slv_primobj] = getprimalobj(task, MSK_SOL_ITR)
    sol_info[:slv_dualobj] = getdualobj(task, MSK_SOL_ITR)
    sol_info[:slv_primfeas] = getdouinf(task, MSK_DINF_INTPNT_PRIMAL_FEAS)
    sol_info[:slv_dualfeas] = getdouinf(task, MSK_DINF_INTPNT_DUAL_FEAS)
    sol_info[:slv_solvetime] = getdouinf(task, MSK_DINF_OPTIMIZER_TIME)
  else
    sol_info[:slv_primobj] = -1.0
    sol_info[:slv_dualobj] = -1.0
    sol_info[:slv_primfeas] = -1.0
    sol_info[:slv_dualfeas] = -1.0
    sol_info[:slv_solvetime] = -1.0
  end
end

function get_primalsol!(task::Mosek.Task, problem::SDP_Problem,
                                          primal::SortedDict{Tuple{String,String,String}, Float64})
  for (id, block) in problem.id_to_sdpblock
    barx = getbarxj(task, MSK_SOL_ITR, id)
    all_variables = ["" for kv in block.var_to_id]
    for (var, varid) in block.var_to_id
      all_variables[varid] = var
    end
    n = 0
    for j in 1:length(all_variables)
      for i in j:length(all_variables)
        n+=1
        # @printf("%15s%15s%15s%25.10f\n", id_block[2].name, all_variables[i], all_variables[j], barx[n])
        primal[block.name, all_variables[i], all_variables[j]] = barx[n]
      end
    end
  end
end

function get_dualsol!(task::Mosek.Task, problem::SDP_Problem,
                                        dual::SortedDict{Tuple{String, String, String}, Float64})
  for (blockid, block) in problem.id_to_sdpblock
    bars = getbarsj(task, MSK_SOL_ITR, blockid)
    id_to_var = OrderedDict([id=>var for (var,id) in block.var_to_id])
    it = 0
    for j=1:length(id_to_var), i=j:length(id_to_var)
      it += 1
      # println("($(block.name), $(id_to_var[i]), $(id_to_var[j])) = $(bars[it])")
      dual[(block.name, id_to_var[i], id_to_var[j])] = bars[it]
    end
  end
end

function dump_mosek_model(task)
  println("*********************************************************************************")
  println("Debug -> Reading Mosek params")
  println("MSK_DPAR_INTPNT_CO_TOL_DFEAS,  $(getdouparam(task, MSK_DPAR_INTPNT_CO_TOL_DFEAS))")
  println("MSK_DPAR_INTPNT_CO_TOL_INFEAS,   $(getdouparam(task, MSK_DPAR_INTPNT_CO_TOL_INFEAS))")
  println("MSK_DPAR_INTPNT_CO_TOL_MU_RED,   $(getdouparam(task, MSK_DPAR_INTPNT_CO_TOL_MU_RED))")
  println("MSK_DPAR_INTPNT_CO_TOL_PFEAS,  $(getdouparam(task, MSK_DPAR_INTPNT_CO_TOL_PFEAS))")
  println("MSK_DPAR_INTPNT_CO_TOL_REL_GAP,  $(getdouparam(task, MSK_DPAR_INTPNT_CO_TOL_REL_GAP))")

  println("\nDebug -> Reading Mosek problem")
  num, subcj, subck, subcl, valcjkl = getbarcblocktriplet(task)
  @printf("%5s  %5s  %5s  %s\n", "subcj", "subck", "subcl", "valcjkl")
  for i=1:num
      @printf("%5i  %5i  %5i  %f\n", subcj[i], subck[i], subcl[i], valcjkl[i])
  end

  num, subai, subaj, subak, subal, valajkl = getbarablocktriplet(task)
  @printf("\n%5s  %5s  %5s  %5s  %s\n", "subai", "subaj", "subak", "subal", "valajkl")
  for i=1:num
      @printf("%5i  %5i  %5i  %5i  %f\n", subai[i], subaj[i], subak[i], subal[i], valajkl[i])
  end

  boundkeys, lbs, ubs = getconboundslice(task, 1, numcon+1)
  @show boundkeys
  @show lbs
  @show ubs
  println("*********************************************************************************")
end

function get_solsta_map()
  return OrderedDict(MSK_SOL_STA_UNKNOWN => "Status of the solution is unknown.",
                     MSK_SOL_STA_OPTIMAL => "The solution is optimal.",
                     MSK_SOL_STA_PRIM_FEAS => "The solution is primal feasible.",
                     MSK_SOL_STA_DUAL_FEAS => "The solution is dual feasible.",
                     MSK_SOL_STA_PRIM_AND_DUAL_FEAS => "The solution is both primal and dual feasible.",
                     MSK_SOL_STA_NEAR_OPTIMAL => "The solution is nearly optimal.",
                     MSK_SOL_STA_NEAR_PRIM_FEAS => "The solution is nearly primal feasible.",
                     MSK_SOL_STA_NEAR_DUAL_FEAS => "The solution is nearly dual feasible.",
                     MSK_SOL_STA_NEAR_PRIM_AND_DUAL_FEAS => "The solution is nearly both primal and dual feasible.",
                     MSK_SOL_STA_PRIM_INFEAS_CER => "The solution is a certificate of primal infeasibility.",
                     MSK_SOL_STA_DUAL_INFEAS_CER => "The solution is a certificate of dual infeasibility.",
                     MSK_SOL_STA_NEAR_PRIM_INFEAS_CER => "The solution is almost a certificate of primal infeasibility.",
                     MSK_SOL_STA_NEAR_DUAL_INFEAS_CER => "The solution is almost a certificate of dual infeasibility.",
                     MSK_SOL_STA_PRIM_ILLPOSED_CER => "The solution is a certificate that the primal problem is illposed.",
                     MSK_SOL_STA_DUAL_ILLPOSED_CER => "The solution is a certificate that the dual problem is illposed.",
                     MSK_SOL_STA_INTEGER_OPTIMAL => "The primal solution is integer optimal.",
                     MSK_SOL_STA_NEAR_INTEGER_OPTIMAL => "The primal solution is near integer optimal.")
end

function get_prosta_map()
  return OrderedDict(MSK_PRO_STA_UNKNOWN => "Unknown problem status.",
                     MSK_PRO_STA_PRIM_AND_DUAL_FEAS => "The problem is primal and dual feasible.",
                     MSK_PRO_STA_PRIM_FEAS => "The problem is primal feasible.",
                     MSK_PRO_STA_DUAL_FEAS => "The problem is dual feasible.",
                     MSK_PRO_STA_NEAR_PRIM_AND_DUAL_FEAS => "The problem is at least nearly primal and dual feasible.",
                     MSK_PRO_STA_NEAR_PRIM_FEAS => "The problem is at least nearly primal feasible.",
                     MSK_PRO_STA_NEAR_DUAL_FEAS => "The problem is at least nearly dual feasible.",
                     MSK_PRO_STA_PRIM_INFEAS => "The problem is primal infeasible.",
                     MSK_PRO_STA_DUAL_INFEAS => "The problem is dual infeasible.",
                     MSK_PRO_STA_PRIM_AND_DUAL_INFEAS => "The problem is primal and dual infeasible.",
                     MSK_PRO_STA_ILL_POSED => "The problem is ill-posed. For example, it may be primal and dual feasible but have a positive duality gap.",
                     MSK_PRO_STA_PRIM_INFEAS_OR_UNBOUNDED => "The problem is either primal infeasible or unbounded. This may occur for mixed-integer problems.")
end


text:
using Pkg; Pkg.activate("~/.julia/dev/Tullio/benchmarks/02/")
pkg"add Tullio, LoopVectorization, Compat, PaddedMatrices, StructArrays, LinearAlgebra, BenchmarkTools, Libdl, MKL_jll, OpenBLAS_jll"
pkg"up"

Threads.nthreads() == 6

using Tullio, LoopVectorization, Compat

TVERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(pkgdir(Tullio), "Project.toml"))["version"])

tmul!(C,A,B) = @tullio C[i,j] := A[i,k] * B[k,j]

# Adapted from:
# https://github.com/chriselrod/PaddedMatrices.jl/blob/master/benchmark/blasbench.jl

using PaddedMatrices, StructArrays, LinearAlgebra, BenchmarkTools, Libdl

randa(::Type{T}, dim...) where {T} = rand(T, dim...)
randa(::Type{T}, dim...) where {T <: Signed} = rand(T(-100):T(200), dim...)

using MKL_jll, OpenBLAS_jll

const libMKL = Libdl.dlopen(MKL_jll.libmkl_rt)
const DGEMM_MKL = Libdl.dlsym(libMKL, :dgemm)
const SGEMM_MKL = Libdl.dlsym(libMKL, :sgemm)
const DGEMV_MKL = Libdl.dlsym(libMKL, :dgemv)
const MKL_SET_NUM_THREADS = Libdl.dlsym(libMKL, :MKL_Set_Num_Threads)

const libOpenBLAS = Libdl.dlopen(OpenBLAS_jll.libopenblas)
const DGEMM_OpenBLAS = Libdl.dlsym(libOpenBLAS, :dgemm_64_)
const SGEMM_OpenBLAS = Libdl.dlsym(libOpenBLAS, :sgemm_64_)
const DGEMV_OpenBLAS = Libdl.dlsym(libOpenBLAS, :dgemv_64_)
const OPENBLAS_SET_NUM_THREADS = Libdl.dlsym(libOpenBLAS, :openblas_set_num_threads64_)

istransposed(x) = 'N'
istransposed(x::Adjoint{<:Real}) = 'T'
istransposed(x::Adjoint) = 'C'
istransposed(x::Transpose) = 'T'

for (lib,f) ∈ [(:GEMM_MKL,:gemmmkl!), (:GEMM_OpenBLAS,:gemmopenblas!)]
    for (T,prefix) ∈ [(Float32,:S),(Float64,:D)]
        fm = Symbol(prefix, lib)
        @eval begin
            function $f(C::AbstractMatrix{$T}, A::AbstractMatrix{$T}, B::AbstractMatrix{$T})
                transA = istransposed(A)
                transB = istransposed(B)
                M, N = size(C); K = size(B, 1)
                pA = parent(A); pB = parent(B)
                ldA = stride(pA, 2)
                ldB = stride(pB, 2)
                ldC = stride(C, 2)
                α = one($T)
                β = zero($T)
                ccall(
                    $fm, Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{$T}, Ref{$T},
                     Ref{Int64}, Ref{$T}, Ref{Int64}, Ref{$T}, Ref{$T}, Ref{Int64}),
                    transA, transB, M, N, K, α, pA, ldA, pB, ldB, β, C, ldC
                )
            end
        end
    end
end
mkl_set_num_threads(N::Integer) = ccall(MKL_SET_NUM_THREADS, Cvoid, (Int32,), N % Int32)
# mkl_set_num_threads(1)
openblas_set_num_threads(N::Integer) = ccall(OPENBLAS_SET_NUM_THREADS, Cvoid, (Int64,), N)
# openblas_set_num_threads(1)

function benchmark_fun!(f!, C, A, B, force_belapsed = false, reference = nothing)
    tmin = @elapsed f!(C, A, B)
    if force_belapsed || 2tmin < BenchmarkTools.DEFAULT_PARAMETERS.seconds
        tmin = min(tmin, @belapsed $f!($C, $A, $B))
    elseif tmin < BenchmarkTools.DEFAULT_PARAMETERS.seconds
        tmin = min(tmin, @elapsed f!(C, A, B))
        if tmin < 2BenchmarkTools.DEFAULT_PARAMETERS.seconds
            tmin = min(tmin, @elapsed f!(C, A, B))
        end
    end
    isnothing(reference) || @assert C ≈ reference
    tmin
end

function runbench(::Type{T}, sizes = [2:255..., round.(Int, range(57.16281374121401, length=200) .^ 1.3705658916944428)...]) where {T}
    (StructVector ∘ map)(sizes) do sz
        n, k, m = sz, sz, sz
        C1 = Matrix{T}(undef, n, m)
        C2 = similar(C1);
        C3 = similar(C1);
        C4 = similar(C1);
        # C5 = similar(C1);
        # C6 = similar(C1);
        A  = randa(T, n, k)
        B  = randa(T, k, m)
        # BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.05

        tmlt = benchmark_fun!(tmul!, C1, A, B, sz == first(sizes))

        jmlt = benchmark_fun!(PaddedMatrices.jmul!, C1, A, B, sz == first(sizes))
        # res = if T <: Integer
        #     (matrix_size=sz, MaBLAS_24x9=ma24, MaBLAS_32x6=ma32, MaBLAS_40x5=ma40, PaddedMatrices=jmlt)
        # else
            opbt = benchmark_fun!(gemmopenblas!, C2, A, B, sz == first(sizes), C1)
            mklbt= benchmark_fun!(gemmmkl!, C3, A, B, sz == first(sizes), C1)

            res = (matrix_size=sz, OpenBLAS=opbt, MKL=mklbt, PaddedMatrices=jmlt, Tullio=tmlt)
        # end
        @show res
    end
end

b5 = runbench(Float32, [3,10,30,100,300])

using Plots

function makeplot(res, title="")
    plt = plot()
    for lab in propertynames(res)[2:end]
        times = getproperty(res, lab)
        flops = 2e-9 * res.matrix_size.^3 ./ times
        str = lab==:MKL ? "MKL $(Compat.get_num_threads())" :
            lab==:OpenBLAS ? "OpenBLAS $(Compat.get_num_threads())" :
            lab==:Tullio ? "Tullio $(Threads.nthreads())" :
            string(lab)
        lab==:Tullio ?
            plot!(res.matrix_size, flops, lab=str, m=:circle) :
            plot!(res.matrix_size, flops, lab=str)
    end
    plot!(yaxis=("gigaflops", ([12.5,25,50,100,200,400],["12.5",25,50,100,200,400]), :log10), xaxis=("size", :log10), legend=:bottomright)
    # plot!(1:0, 1:0, c=:white, lab="i7-8700 + $VERSION")
    plot!(1:0, 1:0, c=:white, lab="Intel " * split(Sys.cpu_info()[1].model, " ")[3], title=title * "Julia " * string(VERSION))
end

makeplot(b5, "warmup ")

# Threads.nthreads() <= 6 || error("expected to run with at most 6 threads")
mkl_set_num_threads(Threads.nthreads())
openblas_set_num_threads(Threads.nthreads())

for Ty in [Float64, Float32]
    global b36

    b36 = runbench(Ty, [10,11,12, 20,21, 30,31,32,33, 49,50,51, 63,64,65,66, 77,78, 100,101,102, 127,128,129, 200, 255,256,257, 300, 400, 500, 511,512,513, 600, 700, 800, 999,1000,1024,1025, 1600,1601, 1999,2000])

    p36 = makeplot(b36, "$Ty, Tullio $TVERSION, ")
    savefig(p36, joinpath("~/.julia/dev/Tullio", "benchmarks/02/matmul-$TVERSION-$Ty-$VERSION.png"))

end

# Summary:

# it's fine at small sizes, sometimes beating OpenBLAS, thanks entirely to @avx
# Threading helps at large sizes, but it still ends up about half the speed.
# Which is OK, the goal isn't to replace BLAS, it's to do other things!
# This is just a test to see how much is left on the table.

ext: jl


text:
# ---
# title: 546. Remove Boxes
# id: problem546
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Depth-first Search
# link: <https://leetcode.com/problems/remove-boxes/description/>
# hidden: true
# ---
#
# Given several boxes with different colors represented by different positive
# numbers.
# You may experience several rounds to remove boxes until there is no box left.
# Each time you can choose some continuous boxes with the same color (composed
# of k boxes, k >= 1), remove them and get `k*k` points.
# Find the maximum points you can get.
#
#
#
# **Example 1:**
#
#
#
#     Input: boxes = [1,3,2,2,2,3,4,3,1]
#     Output: 23
#     Explanation:
#     [1, 3, 2, 2, 2, 3, 4, 3, 1]
#     ----> [1, 3, 3, 4, 3, 1] (3*3=9 points)
#     ----> [1, 3, 3, 3, 1] (1*1=1 points)
#     ----> [1, 1] (3*3=9 points)
#     ----> [] (2*2=4 points)
#
#
#
#
# **Constraints:**
#
#   * `1 <= boxes.length <= 100`
#   * `1 <= boxes[i] <= 100`
#
#
## @lc code=start
using LeetCode

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


text:
abstract type AbstractParams end
abstract type AbstractModel end
abstract type AbstractSetup end

struct Setup <: AbstractSetup
    n_polymers::Int64
    n_monomers::Int64
    n_microbes::Int64
    n_enzymes::Int64
    n_minerals::Int64
    dim::Int64
end

function Setup(n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
    if n_minerals == 0
        dim = n_polymers + n_monomers + 3*n_microbes + n_enzymes
    else
        dim = n_polymers + 2*n_monomers + 3*n_microbes + 2*n_enzymes
    end
    return Setup(n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals, dim)
end

@selectable struct Params{SE,ME,AS,DEP,TU} <: AbstractParams
    # Field                     | Selectable Types
    setup_pars::SE              | AbstractSetup
    metabolism_pars::ME         | AbstractMetabolism
    assimilation_pars::AS       | AbstractAssimilation
    depolymerization_pars::DEP  | Union{nothing, AbstractDepolymerization}
    turnover_pars::TU           | AbstractTurnover
end

for fn in fieldnames(Params)
    @eval $fn(p::Params) = p.$fn
end

function split_state_batch(u::AbstractVector{<:Real}, p::AbstractParams)
    n_polymers = p.setup_pars.n_polymers
    n_monomers = p.setup_pars.n_monomers
    n_microbes = p.setup_pars.n_microbes
    n_enzymes  = p.setup_pars.n_enzymes

    D    = u[1+n_polymers:n_polymers+n_monomers]
    E    = u[1+n_polymers+n_monomers:n_polymers+n_monomers+n_microbes]
    V    = u[1+n_polymers+n_monomers+n_microbes:n_polymers+n_monomers+2*n_microbes]
    X    = u[1+n_polymers+n_monomers+2*n_microbes:n_polymers+n_monomers+2*n_microbes+n_enzymes]
    CO2  = u[1+n_polymers+n_monomers+2*n_microbes+n_enzymes:n_polymers+n_monomers+2*n_microbes+n_enzymes+n_microbes]
    return D, E, V, X, CO2
end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
"""
    FuzzyART.jl

Description:
    Includes all of the structures and logic for running a Gamma-Normalized Fuzzy ART module.

References:
[1] G. Carpenter, S. Grossberg, and D. Rosen, "Fuzzy ART: Fast stable learning and categorization of analog patterns by an adaptive resonance system," Neural Networks, vol. 4, no. 6, pp. 759-771, 1991.
"""

# --------------------------------------------------------------------------- #
# OPTIONS
# --------------------------------------------------------------------------- #

"""
    opts_FuzzyART(;kwargs)

Gamma-Normalized Fuzzy ART options struct.

# Keyword Arguments
- `rho::Float`: vigilance value, [0, 1], default 0.6.
- `alpha::Float`: choice parameter, alpha > 0, default 1e-3.
- `beta::Float`: learning parameter, (0, 1], default 1.0.
- `gamma::Float`: "pseudo" kernel width, gamma >= 1, default 3.0.
- `gamma_ref::Float`: "reference" kernel width, 0 <= gamma_ref < gamma, default 1.0.
- `display::Bool`: display flag, default true.
- `max_epoch::Int`: maximum number of epochs during training, default 1.
- `gamma_normalization::Bool`: normalize the threshold by the feature dimension, default false.
"""
@with_kw mutable struct opts_FuzzyART <: ARTOpts @deftype Float
    # Vigilance parameter: [0, 1]
    rho = 0.6; @assert rho >= 0.0 && rho <= 1.0
    # Choice parameter: alpha > 0
    alpha = 1e-3; @assert alpha > 0.0
    # Learning parameter: (0, 1]
    beta = 1.0; @assert beta > 0.0 && beta <= 1.0
    # "Pseudo" kernel width: gamma >= 1
    gamma = 3.0; @assert gamma >= 1.0
    # "Reference" gamma for normalization: 0 <= gamma_ref < gamma
    gamma_ref = 1.0; @assert 0.0 <= gamma_ref && gamma_ref <= gamma
    # Display flag
    display::Bool = true
    # Maximum number of epochs during training
    max_epochs::Int = 1
    # Normalize the threshold by the feature dimension
    gamma_normalization::Bool = false
end # opts_FuzzyART

# --------------------------------------------------------------------------- #
# STRUCTS
# --------------------------------------------------------------------------- #

"""
    FuzzyART <: ART

Gamma-Normalized Fuzzy ART learner struct

For module options, see [`AdaptiveResonance.opts_FuzzyART`](@ref).

# Option Parameters
- `opts::opts_FuzzyART`: FuzzyART options struct.
- `config::DataConfig`: data configuration struct.

# Working Parameters
- `threshold::Float`: operating module threshold value, a function of the vigilance parameter.
- `labels::IntegerVector`: incremental list of labels corresponding to each F2 node, self-prescribed or supervised.
- `T::RealVector`: activation values for every weight for a given sample.
- `M::RealVector`: match values for every weight for a given sample.
- `W::RealMatrix`: category weight matrix.
- `n_instance::IntegerVector`: number of weights associated with each category.
- `n_categories::Int`: number of category weights (F2 nodes).
- `epoch::Int`: current training epoch.

# References
1. G. Carpenter, S. Grossberg, and D. Rosen, "Fuzzy ART: Fast stable learning and categorization of analog patterns by an adaptive resonance system," Neural Networks, vol. 4, no. 6, pp. 759-771, 1991.
"""
mutable struct FuzzyART <: ART
    # Assign numerical parameters from options
    opts::opts_FuzzyART
    config::DataConfig

    # Working variables
    threshold::Float
    labels::IntegerVector
    T::RealVector
    M::RealVector

    # "Private" working variables
    W::RealMatrix
    n_instance::IntegerVector
    n_categories::Int
    epoch::Int
end # FuzzyART <: ART

# --------------------------------------------------------------------------- #
# CONSTRUCTORS
# --------------------------------------------------------------------------- #

"""
    FuzzyART()

Implements a Gamma-Normalized Fuzzy ART learner.

# Examples
```julia-repl
julia> FuzzyART()
FuzzyART
    opts: opts_FuzzyART
    ...
```
"""
function FuzzyART()
    opts = opts_FuzzyART()
    FuzzyART(opts)
end # FuzzyART()

"""
    FuzzyART(;kwargs...)

Implements a Gamma-Normalized Fuzzy ART learner with keyword arguments.

# Examples
```julia-repl
julia> FuzzyART(rho=0.7)
FuzzyART
    opts: opts_FuzzyART
    ...
```
"""
function FuzzyART(;kwargs...)
    opts = opts_FuzzyART(;kwargs...)
    FuzzyART(opts)
end # FuzzyART(;kwargs...)

"""
    FuzzyART(opts::opts_FuzzyART)

Implements a Gamma-Normalized Fuzzy ART learner with specified options.

# Examples
```julia-repl
julia> FuzzyART(opts)
FuzzyART
    opts: opts_FuzzyART
    ...
```
"""
function FuzzyART(opts::opts_FuzzyART)
    FuzzyART(opts,                          # opts
        DataConfig(),                   # config
        0.0,                            # threshold
        Array{Int}(undef,0),            # labels
        Array{Float}(undef, 0),         # T
        Array{Float}(undef, 0),         # M
        Array{Float}(undef, 0, 0),      # W
        Array{Int}(undef, 0),           # n_instance
        0,                              # n_categories
        0                               # epoch
    )
end # FuzzyART(opts::opts_FuzzyART)

"""
    FuzzyART(opts::opts_FuzzyART, sample::RealVector)

Create and initialize a FuzzyART with a single sample in one step.
"""
function FuzzyART(opts::opts_FuzzyART, sample::RealVector ; preprocessed::Bool=false)
    art = FuzzyART(opts)
    init_train!(sample, art, preprocessed)
    initialize!(art, sample)
    return art
end # FuzzyART(opts::opts_FuzzyART, sample::RealVector)

# --------------------------------------------------------------------------- #
# ALGORITHMIC METHODS
# --------------------------------------------------------------------------- #

function set_threshold!(art::FuzzyART)
    if art.opts.gamma_normalization
        art.threshold = art.opts.rho*(art.config.dim^art.opts.gamma_ref)
    else
        art.threshold = art.opts.rho
    end
end # set_threshold!(art::FuzzyART)

"""
    initialize!(art::FuzzyART, x::Vector{T} ; y::Integer=0) where {T<:RealFP}

Initializes a FuzzyART learner with an intial sample 'x'.

# Examples
```julia-repl
julia> my_FuzzyART = FuzzyART()
FuzzyART
    opts: opts_FuzzyART
    ...
julia> initialize!(my_FuzzyART, [1 2 3 4])
```
"""
function initialize!(art::FuzzyART, x::Vector{T} ; y::Integer=0) where {T<:RealFP}
    # Initialize the instance and categories counters
    art.n_instance = [1]
    art.n_categories = 1

    # Set the threshold
    set_threshold!(art)

    # Fast commit the weight
    art.W = Array{T}(undef, art.config.dim_comp, 1)

    # Assign the contents, valid this way for 1-D or 2-D arrays
    art.W[:, 1] = x

    # Set the label to either the supervised label or 1 if unsupervised
    label = !iszero(y) ? y : 1

    # Add the label to the label list
    push!(art.labels, label)
end # initialize!(art::FuzzyART, x::Vector{T} ; y::Integer=0) where {T<:RealFP}

# FuzzyART incremental training method
function train!(art::FuzzyART, x::RealVector ; y::Integer=0, preprocessed::Bool=false)
    # Flag for if training in supervised mode
    supervised = !iszero(y)

    # Run the sequential initialization procedure
    sample = init_train!(x, art, preprocessed)

    # Initialization if weights are empty; fast commit the first sample
    if isempty(art.W)
        y_hat = supervised ? y : 1
        initialize!(art, sample, y=y_hat)
        return y_hat
    end

    # If we have a new supervised category, create a new category
    if supervised && !(y in art.labels)
        create_category(art, sample, y)
        return y
    end

    # Compute activation/match functions
    activation_match!(art, sample)
    # Sort activation function values in descending order
    index = sortperm(art.T, rev=true)
    # Initialize mismatch as true
    mismatch_flag = true

    # Loop over all categories
    for j = 1:art.n_categories
        # Best matching unit
        bmu = index[j]
        # Vigilance check - pass
        if art.M[bmu] >= art.threshold
            # If supervised and the label differed, force mismatch
            if supervised && (art.labels[bmu] != y)
                break
            end
            # Learn the sample
            learn!(art, sample, bmu)
            # Save the output label for the sample
            y_hat = art.labels[bmu]
            # No mismatch
            mismatch_flag = false
            break
        end
    end

    # If there was no resonant category, make a new one
    if mismatch_flag
        # Get the correct label for the new category
        y_hat = supervised ? y : art.n_categories + 1
        # Create a new category
        create_category(art, sample, y_hat)
    end

    return y_hat
end # train!(art::FuzzyART, x::RealVector ; y::Integer=0, preprocessed::Bool=false)

"""
    create_category(art::FuzzyART, x::RealVector, y::Integer)
"""
function create_category(art::FuzzyART, x::RealVector, y::Integer)
    # Increment the number of categories
    art.n_categories += 1
    # Fast commit
    art.W = hcat(art.W, x)
    # Increment number of samples associated with new category
    push!(art.n_instance, 1)
    # Add the label for the ategory
    push!(art.labels, y)
end # create_category(art::FuzzyART, x::RealVector, y::Integer)

# FuzzyART incremental classification method
function classify(art::FuzzyART, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)
    # Preprocess the data
    x = init_classify!(x, art, preprocessed)

    # Compute activation and match functions
    activation_match!(art, x)
    # Sort activation function values in descending order
    index = sortperm(art.T, rev=true)
    # Default is mismatch
    mismatch_flag = true
    y_hat = -1
    for jx in 1:art.n_categories
        bmu = index[jx]
        # Vigilance check - pass
        if art.M[bmu] >= art.threshold
            # Current winner
            y_hat = art.labels[bmu]
            mismatch_flag = false
            break
        end
    end
    # If we did not find a match
    if mismatch_flag
        # Create new weight vector
        @debug "Mismatch"
        # Report either the best matching unit or the mismatch label -1
        y_hat = get_bmu ? art.labels[index[1]] : -1
    end
    return y_hat
end # classify(art::FuzzyART, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)

"""
    activation_match!(art::FuzzyART, x::RealVector)

Computes the activation and match functions of the art module against sample x.

# Examples
```julia-repl
julia> my_FuzzyART = FuzzyART()
FuzzyART
    opts: opts_FuzzyART
    ...
julia> x, y = load_data()
julia> train!(my_FuzzyART, x)
julia> x_sample = x[:, 1]
julia> activation_match!(my_FuzzyART, x_sample)
```
"""
function activation_match!(art::FuzzyART, x::RealVector)
    art.T = zeros(art.n_categories)
    art.M = zeros(art.n_categories)
    for i = 1:art.n_categories
        W_norm = norm(art.W[:, i], 1)
        numerator = norm(element_min(x, art.W[:, i]), 1)
        art.T[i] = (numerator/(art.opts.alpha + W_norm))^art.opts.gamma
        if art.opts.gamma_normalization
            art.M[i] = (W_norm^art.opts.gamma_ref)*art.T[i]
        else
            art.M[i] = numerator/norm(x, 1)
        end
    end
end # activation_match!(art::FuzzyART, x::RealVector)

"""
    learn(art::FuzzyART, x::RealVector, W::RealVector)

Return the modified weight of the art module conditioned by sample x.
"""
function learn(art::FuzzyART, x::RealVector, W::RealVector)
    # Update W
    return art.opts.beta .* element_min(x, W) .+ W .* (1 - art.opts.beta)
end # learn(art::FuzzyART, x::RealVector, W::RealVector)

"""
    learn!(art::FuzzyART, x::RealVector, index::Integer)

In place learning function with instance counting.
"""
function learn!(art::FuzzyART, x::RealVector, index::Integer)
    # Update W
    art.W[:, index] = learn(art, x, art.W[:, index])
    art.n_instance[index] += 1
end # learn!(art::FuzzyART, x::RealVector, index::Integer)

"""
    stopping_conditions(art::FuzzyART)

Stopping conditions for a FuzzyART module.
"""
function stopping_conditions(art::FuzzyART)
    return art.epoch >= art.opts.max_epochs
end # stopping_conditions(art::FuzzyART)


text:
using JSON
const include_re = r"^include\(\"(.+?)\"\)$"
const module_re = r"^module\s.+"

function strip_module_end(content)
    content = strip(content)
    if endswith(content, "end  # module")
        content = strip(content[1:end-length("end  # module")])
    end
    content
end

function process_cell(cell_body)
    if startswith(cell_body, "#=") && endswith(cell_body, "=#")
        return md_cell(strip(match(r"#=(.+)=#"s, cell_body)[1]))
    end
    if startswith(cell_body, "\"\"\"") && endswith(cell_body, "\"\"\"")
        return md_cell(strip(match(r"\"\"\"(.+)\"\"\""s, cell_body)[1]))
    end
    return code_cell(strip_module_end(cell_body))
end

myd(;kwargs...) = Dict(string(k) => v for (k, v) in kwargs)

function md_cell(content)
    source = [line * "\n" for line in split(content, "\n")]
    source[end] = rstrip(source[end])
    myd(
        cell_type="markdown",
        metadata=myd(outputExpanded=false),
        source=source
    )
end

function code_cell(content)
    source = [line * "\n" for line in split(content, "\n")]
    source[end] = rstrip(source[end])
    myd(
        outputs=Any[],
        cell_type="code",
        source=source,
        metadata=myd(outputExpanded=false),
        execution_count=nothing
    )
end

# from a string that contains txt, build a cell
function _build_cell(txt::String)
    out = Dict[]
    pos = 1
    while true
        try
            cell_body, pos = get_cell(txt, pos)
            cell = process_cell(cell_body)
            if isa(cell, Dict)
                push!(out, process_cell(cell_body))
            else
                break
            end catch e
            if isa(e, BoundsError)
                break
            else
                rethrow(e)
            end
        end
    end
    out
end

function build_cells(fn::AbstractString)
    @assert isfile(fn)
    out = Dict[]

    # read in the file as a vector of lines
    lines = readlines(fn)

    # track which lines are important (start a cell or include a file)
    important_lines = Tuple{Int,Symbol}[]
    for (i, l) in enumerate(lines)
        if occursin(module_re, l)
            push!(important_lines, (i, :skip))
        end
        if l == "#%% cell"
            push!(important_lines, (i, :cell))
        end
        if l == "#%% skip"
            push!(important_lines, (i, :skip))
        end
        if occursin(include_re, l)
            push!(important_lines, (i, :include))
        end
    end

    # add the end of the file
    push!(important_lines, (length(lines)+1, :end))

    for (i, (l, typ)) in enumerate(important_lines[1:end-1])
        lp, typ_p = important_lines[i+1]
        if typ == :cell
            txt = strip(join(lines[l+1:lp-1], '\n'))
            push!(out, process_cell(txt))
        elseif typ == :include
            append!(out, build_cells(match(include_re, lines[l])[1]))
        elseif typ == :skip
            continue
        else
            error("Handle this type of important line")
        end
    end

    return out

end

function create_notebook(fn::String)
    @assert isfile(fn)
    cells = build_cells(fn)
    if endswith(basename(fn), ".jl")
        return nb = myd(
            nbformat=4,
            nbformat_minor=1,
            metadata = myd(
                kernelspec = myd(
                    name = "julia-1.1",
                    display_name = "Julia 1.1",
                    language="julia"
                ),
                language_info = myd(
                    file_extension = ".jl",
                    mimetype = "application/julia",
                    name = "julia",
                    version = "1.1"
                )
            ),
            cells=build_cells(fn)
        )
    elseif endswith(basename(fn), ".py")
        return myd(
            nbformat=4,
            nbformat_minor=1,
            metadata = myd(
                kernelspec = myd(
                    name = "python3",
                    display_name = "Python 3",
                    language="python"
                ),
                language_info = myd(
                    file_extension = ".py",
                    mimetype = "text/x-python",
                    name = "python",
                    nbconvert_exporter="python",
                    pygments_lexer="ipython3",
                    codemirror_mode=myd(
                        name="ipython",
                        version=3
                    ),
                    version="3.6.0"
                )
            ),
            cells=build_cells(fn)
        )
    end
end

function save_nb(fn, dest)
    open(dest, "w") do f
        print(f, JSON.json(create_notebook(fn)))
    end
end

function main()
    # save julia file
    this_dir = dirname(@__FILE__)
    nb_dir = joinpath(dirname(this_dir), "notebooks")
    py_dir = joinpath(dirname(this_dir), "python")

    save_nb(joinpath(this_dir, "main.jl"), joinpath(nb_dir, "julia.ipynb"))
    save_nb(joinpath(py_dir, "main.py"), joinpath(nb_dir, "python.ipynb"))
end

max_stars_repo_name: sglyon/CLMMJuliaPythonMatlab

max_issues_repo_name: sglyon/CLMMJuliaPythonMatlab

max_forks_repo_name: sglyon/CLMMJuliaPythonMatlab


text:
module GeoJSONTables

import JSON3, Tables, GeoInterface

struct FeatureCollection{T} <: AbstractVector{eltype(T)}
    json::T
end

function read(source)
    fc = JSON3.read(source)
    features = get(fc, :features, nothing)
    if get(fc, :type, nothing) == "FeatureCollection" && features isa JSON3.Array
        FeatureCollection{typeof(features)}(features)
    else
        throw(ArgumentError("input source is not a GeoJSON FeatureCollection"))
    end
end

Tables.istable(::Type{<:FeatureCollection}) = true
Tables.rowaccess(::Type{<:FeatureCollection}) = true
Tables.rows(fc::FeatureCollection) = fc

Base.IteratorSize(::Type{<:FeatureCollection}) = Base.HasLength()
Base.length(fc::FeatureCollection) = length(json(fc))
Base.IteratorEltype(::Type{<:FeatureCollection}) = Base.HasEltype()

# read only AbstractVector
Base.size(fc::FeatureCollection) = size(json(fc))
Base.getindex(fc::FeatureCollection, i) = Feature(json(fc)[i])
Base.IndexStyle(::Type{<:FeatureCollection}) = IndexLinear()

miss(x) = ifelse(x === nothing, missing, x)

struct Feature{T}
    json::T
end

# these features always have type="Feature", so exclude that
# the keys in properties are added here for direct access
Base.propertynames(f::Feature) = keys(properties(f))

"Access the properties JSON3.Object of a Feature"
properties(f::Feature) = json(f).properties
"Access the JSON3.Object that represents the Feature"
json(f::Feature) = getfield(f, :json)
"Access the JSON3.Array that represents the FeatureCollection"
json(f::FeatureCollection) = getfield(f, :json)
"Access the JSON3.Object that represents the Feature's geometry"
geometry(f::Feature) = json(f).geometry

"""
Get a specific property of the Feature

Returns missing for null/nothing or not present, to work nicely with
properties that are not defined for every feature. If it is a table,
it should in some sense be defined.
"""
function Base.getproperty(f::Feature, nm::Symbol)
    props = properties(f)
    val = get(props, nm, missing)
    miss(val)
end

@inline function Base.iterate(fc::FeatureCollection)
    st = iterate(json(fc))
    st === nothing && return nothing
    val, state = st
    return Feature(val), state
end

@inline function Base.iterate(fc::FeatureCollection, st)
    st = iterate(json(fc), st)
    st === nothing && return nothing
    val, state = st
    return Feature(val), state
end

Base.show(io::IO, fc::FeatureCollection) = println(io, "FeatureCollection with $(length(fc)) Features")
function Base.show(io::IO, f::Feature)
    println(io, "Feature with geometry type $(geometry(f).type) and properties $(propertynames(f))")
end
Base.show(io::IO, ::MIME"text/plain", fc::FeatureCollection) = show(io, fc)
Base.show(io::IO, ::MIME"text/plain", f::Feature) = show(io, f)

include("geointerface.jl")

end # module


text:
using CALiPPSO
using Statistics, Random
using CALiPPSO.Distributions

Random.seed!(456) # Just for reproducibility

#= UNcomment the following lines if you want to use Gurobi Solver, for running the tests.
    Or adjust them to the solver of your choice
=#
# using Gurobi
# const solver = Gurobi
# const tol_overlap=1e-8 # tolerance for identifying an Overlap. Given that Gurobi's precision is 10^-9, a larger value is needed.
# const tol_optimality = 1e-9; # optimality tolerances. This is the most precise value allowed by Gurobi
# const solver_args = Gurobi.Env()
# const solver_attributes = Dict("OutputFlag" => 0, "FeasibilityTol" => tol_optimality, "OptimalityTol" => tol_optimality, "Method" => 3, "Threads" =>  CALiPPSO.max_threads)

#= Comment the following lines if you want to use CALiPPSO with a different solver than GLPK; e.g., Gurobi (see above)=#
const solver = CALiPPSO.default_solver
const tol_overlap = CALiPPSO.default_tol_overlap
const tol_optimality = CALiPPSO.default_tol_optimality
const solver_args = CALiPPSO.default_args
const solver_attributes = CALiPPSO.default_solver_attributes

precompile_main_function(solver, solver_attributes, solver_args)

const Nrand = 512 # size of configurations to test
const L = 1.0 # size of each side of the system's volume

##################################################
##################################################
# First a test on a 2d BIdisperse configuration
##################################################
##################################################
const d_bi = 2 # dimension for testing BIdisperse packings
const N1 = ceil(Int, Nrand/2)
const ϕ_bi = 0.4
r1, r2, Xs0_bi = generate_random_configuration(d_bi, N1, N1, ϕ_bi, 1.4, L); # this creates a bidisperse configuration  with radii ratio of r2/r1= 1.4
ℓ0_bi = min(L, 6*r2)

printstyled("\n\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:blue)
printstyled("Using CALiPPSO to jam a BIdisperse system of N = ", 2N1, " in d = ", d_bi, "\t Initial (φ, R1, R2) = ", [ϕ_bi, r1, r2], "\n", bold=:true, color=:blue)
printstyled("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:blue)

Rs_bi = [r1*ones(N1); r2*ones(N1)]

@time jammed_packing_bi, info_convergence_bi, Γs_vs_t_bi, smax_vs_t_bi, iso_vs_t_bi = produce_jammed_configuration!(Xs0_bi, Rs_bi; ℓ0=ℓ0_bi, sqrΓ0=1.5,  max_iters=1000, initial_monitor=20, solver=solver, solver_attributes=solver_attributes, solver_args=solver_args)
println("_______________________________________________________________________________________\n\n")
times = info_convergence_bi.times_LP_optim
println("Info about LP times")
println("(min, avg ± std, max) LP times:\t", minimum(times), ", \t", mean(times), " ± ", std(times), ", \t", maximum(times))
println("_______________________________________________________________________________________\n\n\n")

network_of_contacts(jammed_packing_bi)

##################################################
##################################################
# Now fully polydisperse packings with the radii drawn from a Log-normal distribution
##################################################
##################################################
const ds = [2, 3, 4] # dimensions to use with POLYdisperse packings
const σRs = sqrt(log(1+0.2^2)) # this parameter determines the width of the radii distribution
const distr_Rs = LogNormal(0, σRs) # radii will be drawn from a LogNormal distribution
const scl_facs = [0.017, 0.05, 0.09]; # scale factors to obtain a feasible and good initial value of φ

for (nd, d) in enumerate(ds)
    # Random.seed!(456) # Just for reproducibility

    Rs = sort!(scl_facs[nd].*rand(distr_Rs, Nrand), rev=true) # by sorting by the largest values, it's easier to obtain a higher density, initial configuration. However, this could introduce unwanted correlations. So use it with care
    ϕ, Xs0 =generate_random_configuration(d, Rs, L)

    μR = mean(Rs); Rmin = minimum(Rs); Rmax = maximum(Rs)
    ℓ0 = min(L, 6*Rmax)

    printstyled("\n\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:blue)
    printstyled("Using CALiPPSO to jam a POLYdisperse system of N = ", Nrand, " in d = ", d, "\t Initial (φ, Rₘᵢₙ, ⟨R⟩, Rₘₐₓ) = ", [ϕ, Rmin, μR, Rmax], "\n", bold=:true, color=:blue)
    printstyled("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:blue)

    @time jammed_packing, info_convergence, Γs_vs_t, smax_vs_t, iso_vs_t = produce_jammed_configuration!(Xs0, Rs; ℓ0=ℓ0, sqrΓ0=1.5,  max_iters=3000, initial_monitor=20, solver=solver, solver_attributes=solver_attributes, solver_args=solver_args)
    println("_______________________________________________________________________________________\n\n")
    local times = info_convergence.times_LP_optim
    println("Info about LP times")
    # println(times)
    println("(min, avg ± std, max) LP times:\t", minimum(times), ", \t", mean(times), " ± ", std(times), ", \t", maximum(times))
    println("_______________________________________________________________________________________\n\n\n")

end

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO