text:
using Base.Threads
using LinearAlgebra
using HDF5
using PrettyTables
using Printf

const do_continue_print = false
const print_eri = false

function rhf_energy(mol::Molecule, basis::Basis,
  scf_flags::Union{Dict{String,Any},Dict{Any,Any}}; output)

  debug::Bool = haskey(scf_flags, "debug") ? scf_flags["debug"] : false
  niter::Int = haskey(scf_flags, "niter") ? scf_flags["niter"] : 50
  guess::String = haskey(scf_flags, "guess") ? scf_flags["guess"] : "sad"

  ndiis::Int = haskey(scf_flags, "ndiis") ? scf_flags["ndiis"] : 10
  dele::Float64 = haskey(scf_flags, "dele") ? scf_flags["dele"] : 1E-6
  rmsd::Float64 = haskey(scf_flags, "rmsd") ? scf_flags["rmsd"] : 1E-6
  load::String = haskey(scf_flags, "load") ? scf_flags["load"] : "static"
  fdiff::Bool = haskey(scf_flags, "fdiff") ? scf_flags["fdiff"] : false

  return rhf_kernel(mol,basis; output=output, debug=debug,
    niter=niter, guess=guess, ndiis=ndiis, dele=dele, rmsd=rmsd, load=load,
    fdiff=fdiff)
end

"""
	 rhf_kernel(FLAGS::RHF_Flags, basis::Basis, read_in::Dict{String,Any},
       type::T)
Summary
======
Perform the core RHF SCF algorithm.

Arguments
======
FLAGS = Input flags

basis = Generated basis set

read_in = file required to read in from input file

type = Precision of variables in calculation
"""
function rhf_kernel(mol::Molecule,
  basis::Basis;
  output::Int64, debug::Bool, niter::Int, guess::String, ndiis::Int,
  dele::Float64, rmsd::Float64, load::String, fdiff::Bool)

  comm=MPI.COMM_WORLD
  calculation_status = Dict([])

  #== read in some variables from scf input ==#
  debug_output = debug ? h5open("debug.h5","w") : nothing

  #== compute nuclear repulsion energy ==#
  E_nuc = compute_enuc(mol)

  jeri_oei_engine = JERI.OEIEngine(mol.mol_cxx,
    basis.basis_cxx, 0)

  #== compute one-electron integrals and Hamiltonian ==#
  S = zeros(Float64, (basis.norb, basis.norb))
  compute_overlap(S, basis, jeri_oei_engine)

  #for i in 1:basis.norb, j in 1:i
  #  println("OVR($i,$j): ", S[i,j])
  #end

  T = zeros(Float64, (basis.norb, basis.norb))
  compute_ke(T, basis, jeri_oei_engine)

  V = zeros(Float64, (basis.norb, basis.norb))
  compute_nah(V, mol, basis, jeri_oei_engine)

  H = T .+ V

  #== compute initial guess ==#
  guess_matrix = guess == "sad" ? sad_guess(mol, basis) : deepcopy(H)

  if debug && MPI.Comm_rank(comm) == 0
    h5write("debug.h5","RHF/Iteration-None/E_nuc", E_nuc)
    h5write("debug.h5","RHF/Iteration-None/S", S)
    h5write("debug.h5","RHF/Iteration-None/T", T)
    h5write("debug.h5","RHF/Iteration-None/V", V)
    h5write("debug.h5","RHF/Iteration-None/H", H)
    h5write("debug.h5","RHF/Iteration-None/Guess", guess_matrix)
  end

  #== build the initial matrices ==#
  D = guess == "sad" ? guess_matrix : zeros(size(H))
  F = guess == "hcore" ? guess_matrix : zeros(size(H))

  F_eval = Vector{Float64}(undef,basis.norb)
  F_evec = similar(F)

  W = similar(F)

  C = similar(F)

  #== allocate workspace matrices ==#
  workspace_a = similar(F)
  workspace_b = similar(F)
  workspace_c = [ similar(F) ]

  #== check for linear dependencies ==#
  minimum_eval_threshold = 1.0E-6
  LinearAlgebra.BLAS.blascopy!(length(S), S, 1, workspace_b, 1)
  #workspace_b .= S
  S_eval_diag, workspace_a[:,:] = eigen!(LinearAlgebra.Hermitian(workspace_b))

  #fill!(workspace_b, 0.0)
  LinearAlgebra.BLAS.scal!(length(workspace_b), 0.0, workspace_b, 1)
  for i in 1:basis.norb
    workspace_b[i,i] = S_eval_diag[i]
  end

  S_good_evals = findall(eval -> eval >= minimum_eval_threshold, S_eval_diag)
  #ortho = similar(F)
  #LinearAlgebra.BLAS.gemm!('N', 'T', 1.0,
  #  workspace_b, workspace_a, 0.0, ortho)
  #LinearAlgebra.BLAS.gemm!('N', 'N', 1.0, workspace_a, ortho, 0.0, ortho)

  #== build the orthogonalization matrix ==#
  @views ortho = workspace_a[1:end,S_good_evals]*
    (LinearAlgebra.Diagonal(workspace_b[S_good_evals,S_good_evals])^-0.5)*
    transpose(workspace_a[1:end, S_good_evals])

  if debug && MPI.Comm_rank(comm) == 0
    h5write("debug.h5","RHF/Iteration-None/X", ortho)
  end

  if MPI.Comm_rank(comm) == 0 && output >= 2
    println("----------------------------------------          ")
    println("       Starting RHF iterations...                 ")
    println("----------------------------------------          ")
    println(" ")
    println("Iter        Energy                ΔE                Drms")
  end

  E_elec = 0.0
  F_eval = zeros(size(F)[1])
  if guess == "hcore"
    E_elec, F_eval[:] = iteration(F, D, C, H, F_eval, F_evec, workspace_a,
      workspace_b, ortho, basis, 0, debug)
  end

  F_old = deepcopy(F)

  E = 0.0
  E_old = E

  #if MPI.Comm_rank(comm) == 0 && output >= 2
  #  @printf("0     %.10f\n", E)
  #end

  #=============================#
  #== start scf cycles: #7-10 ==#
  #=============================#
  F, D, W, C, E, converged = scf_cycles(F, D, W, C, E, H, ortho, S,
    F_eval, F_evec, F_old, workspace_a, workspace_b, workspace_c,
    E_nuc, E_elec, E_old, basis;
    output=output, debug=debug, niter=niter, ndiis=ndiis, dele=dele,
    rmsd=rmsd, load=load, fdiff=fdiff)

  if !converged
    if MPI.Comm_rank(comm) == 0 && output >= 1
      println(" ")
      println("----------------------------------------")
      println(" The SCF calculation did not converge.  ")
      println("      Restart data is being output.     ")
      println("----------------------------------------")
      println(" ")
    end

    calculation_fail = Dict(
    "success" => false,
    "error" => Dict(
      "error_type" => "convergence_error",
      "error_message" => " SCF calculation did not converge within $niter
        iterations. "
      )
    )

    merge!(calculation_status, calculation_fail)

  else
    if MPI.Comm_rank(comm) == 0 && output >= 1
      println(" ")
      println("----------------------------------------")
      println("   The SCF calculation has converged!   ")
      println("----------------------------------------")
      #println("Total SCF Energy: ",E," h")
      @printf("Total SCF Energy: %.10f h\n",E)
      println(" ")

      calculation_success = Dict(
      "return_result" => E,
      "success" => true,
      "properties" => Dict(
        "return_energy" => E,
        "nuclear_repulsion_energy" => E_nuc,
        #"scf_iterations" => iter,
        "scf_total_energy" => E
        )
      )

      merge!(calculation_status, calculation_success)
    end
  end

  if debug close(debug_output) end

  scf = Dict("Fock" => F,
             "Density" => D,
             "Energy-Weighted Density" => W,
             "MO Coeff" => C,
             "Overlap" => S,
             "Energy" => E,
             "Converged?" => converged
            )

  return scf
end

function scf_cycles(F::Matrix{Float64}, D::Matrix{Float64},
  W::Matrix{Float64}, C::Matrix{Float64},
  E::Float64, H::Matrix{Float64}, ortho::Matrix{Float64},
  S::Matrix{Float64}, F_eval::Vector{Float64},
  F_evec::Matrix{Float64},  F_old::Matrix{Float64},
  workspace_a::Matrix{Float64}, workspace_b::Matrix{Float64},
  workspace_c::Vector{Matrix{Float64}}, E_nuc::Float64, E_elec::Float64,
  E_old::Float64, basis::Basis;
  output::Int64, debug::Bool, niter::Int, ndiis::Int,
  dele::Float64, rmsd::Float64, load::String, fdiff::Bool)

  #== read in some more variables from scf flags input ==#
  nsh = length(basis)
  nindices = (muladd(nsh,nsh,nsh)*(muladd(nsh,nsh,nsh) + 2)) >> 3

  #== build DIIS arrays ==#
  F_array = fill(similar(F), ndiis)

  e_array = fill(similar(F), ndiis)
  e_array_old = fill(similar(F), max(0,ndiis-1))

  F_array_old = fill(similar(F), max(0,ndiis-1))

  #FD = similar(F)
  FDS = similar(F)
  #SDF = similar(F)

  #== build arrays needed for post-fock build iteration calculations ==#
  #F_temp = similar(F)
  ΔF = similar(F)
  F_cumul = zeros(size(F))

  D_old = zeros(size(F))
  ΔD = deepcopy(D)
  D_input = similar(F)

  #== build matrix of Cauchy-Schwarz upper bounds ==#
  schwarz_bounds = zeros(Float64,(nsh,nsh))
  compute_schwarz_bounds(schwarz_bounds, basis, nsh)

  Dsh = similar(schwarz_bounds)

  #== build eri batch arrays ==#
  #eri_sizes::Vector{Int64} = load("tei_batch.jld",
  #  "Sizes/$quartet_batch_num_old")
  #length_eri_sizes::Int64 = length(eri_sizes)

  #@views eri_starts::Vector{Int64} = [1, [ sum(eri_sizes[1:i])+1 for i in 1:(length_eri_sizes-1)]... ]

  #eri_batch::Vector{Float64} = load("tei_batch.jld",
  #  "Integrals/$quartet_batch_num_old")

  #eri_sizes = []
  #eri_starts = []
  #eri_batch = []

  #== execute convergence procedure ==#
  scf_converged = true

  E = scf_cycles_kernel(F, D, W, C, E, H, ortho, S, E_nuc,
    E_elec, E_old, basis, F_array, e_array, e_array_old,
    F_array_old, F_eval, F_evec, F_old, workspace_a,
    workspace_b, workspace_c, ΔF, F_cumul,
    D_old, ΔD, D_input, scf_converged, FDS,
    schwarz_bounds, Dsh;
    output=output, debug=debug, niter=niter, ndiis=ndiis, dele=dele,
    rmsd=rmsd, load=load, fdiff=fdiff)

  #== we are done! ==#
  if debug
    h5write("debug.h5","RHF/Iteration-Final/F", F)
    h5write("debug.h5","RHF/Iteration-Final/D", D)
    h5write("debug.h5","RHF/Iteration-Final/C", C)
    h5write("debug.h5","RHF/Iteration-Final/E", E)
    h5write("debug.h5","RHF/Iteration-Final/converged", scf_converged)
  end

  return F, D, W, C, E, scf_converged
end

function scf_cycles_kernel(F::Matrix{Float64}, D::Matrix{Float64},
  W::Matrix{Float64}, C::Matrix{Float64}, E::Float64, H::Matrix{Float64},
  ortho::Matrix{Float64}, S::Matrix{Float64}, E_nuc::Float64,
  E_elec::Float64, E_old::Float64, basis::Basis,
  F_array::Vector{Matrix{Float64}},
  e_array::Vector{Matrix{Float64}}, e_array_old::Vector{Matrix{Float64}},
  F_array_old::Vector{Matrix{Float64}},
  F_eval::Vector{Float64}, F_evec::Matrix{Float64},
  F_old::Matrix{Float64}, workspace_a::Matrix{Float64},
  workspace_b::Matrix{Float64}, workspace_c::Vector{Matrix{Float64}},
  ΔF::Matrix{Float64},
  F_cumul::Matrix{Float64}, D_old::Matrix{Float64},
  ΔD::Matrix{Float64}, D_input::Matrix{Float64}, scf_converged::Bool,
  FDS::Matrix{Float64},
  schwarz_bounds::Matrix{Float64}, Dsh::Matrix{Float64};
  output, debug, niter, ndiis, dele, rmsd, load, fdiff)

  #== initialize a few more variables ==#
  comm=MPI.COMM_WORLD

  B_dim = 1
  D_rms = 1.0
  ΔE = 1.0
  cutoff = fdiff ? 5E-11 : 1E-10

  #length_eri_sizes = length(eri_sizes)

  #=================================#
  #== now we start scf iterations ==#
  #=================================#
  iter = 1
  iter_converged = false

  nthreads = Threads.nthreads()

  max_am = max_ang_mom(basis)
  eri_quartet_batch_thread = [ Vector{Float64}(undef,
    eri_quartet_batch_size(max_am))
    for thread in 1:nthreads ]

  F_thread = [ zeros(size(F)) for thread in 1:nthreads ]
  jeri_engine_thread = [ JERI.TEIEngine(basis.basis_cxx, basis.shpdata_cxx)
    for thread in 1:nthreads ]

  while !iter_converged
    #== reset eri arrays ==#
    #if quartet_batch_num_old != 1 && iter != 1
    #  resize!(eri_sizes,length_eri_sizes)
    #  resize!(eri_starts,length_eri_sizes)

    #  eri_sizes[:] = load("tei_batch.jld",
  #      "Sizes/$quartet_batch_num_old")

    #  @views eri_starts[:] = [1, [ sum(eri_sizes[1:i])+1 for i in 1:(length_eri_sizes-1)]... ]
      #eri_starts[:] = load("tei_batch.jld",
      #  "Starts/$quartet_batch_num_old")
      #@views eri_starts[:] = eri_starts[:] .- (eri_starts[1] - 1)

    #  resize!(eri_batch,sum(eri_sizes))
    #  eri_batch[:] = load("tei_batch.jld","Integrals/$quartet_batch_num_old")
    #end

    #== determine input D and F ==#
    if fdiff
      LinearAlgebra.BLAS.blascopy!(length(ΔD), ΔD, 1, D_input, 1)
      LinearAlgebra.BLAS.blascopy!(length(ΔF), ΔF, 1, workspace_b, 1)
      #D_input .= fdiff ? ΔD : D
      #workspace_b .= fdiff ? ΔF : F
    else
      LinearAlgebra.BLAS.blascopy!(length(D), D, 1, D_input, 1)
      LinearAlgebra.BLAS.blascopy!(length(F), F, 1, workspace_b, 1)
      #D_input .= fdiff ? ΔD : D
      #workspace_b .= fdiff ? ΔF : F
    end

    #== compress D into shells in Dsh ==#
    for ish in 1:length(basis), jsh in 1:ish
      ipos = basis[ish].pos
      ibas = basis[ish].nbas

      jpos = basis[jsh].pos
      jbas = basis[jsh].nbas

      max_value = 0.0
      for i in ipos:(ipos+ibas-1), j in jpos:(jpos+jbas-1)
        max_value = max(max_value, Base.abs_float(D_input[i,j]))
      end
      Dsh[ish, jsh] = max_value
      Dsh[jsh, ish] = Dsh[ish, jsh]
    end

    #== build new Fock matrix ==#
    workspace_a .= fock_build(workspace_b, F_thread, D_input, H, basis,
      schwarz_bounds, Dsh, eri_quartet_batch_thread, jeri_engine_thread,
      iter, cutoff, debug, load)

    workspace_b .= MPI.Allreduce(workspace_a,MPI.SUM,comm)
    MPI.Barrier(comm)

    if debug && MPI.Comm_rank(comm) == 0
      h5write("debug.h5","RHF/Iteration-$iter/F/Skeleton", workspace_b)
    end

    if fdiff
      ΔF .= workspace_b
      F_cumul .+= ΔF
      F .= F_cumul .+ H
    else
      LinearAlgebra.BLAS.axpy!(1.0, H, workspace_b)
      LinearAlgebra.BLAS.blascopy!(length(workspace_b), workspace_b, 1,
        F, 1)

      #F .= workspace_b .+ H
    end

    if debug && MPI.Comm_rank(comm) == 0
      h5write("debug.h5","RHF/Iteration-$iter/F/Total", F)
    end

    #== do DIIS ==#
    if ndiis > 0
      BLAS.symm!('L', 'U', 1.0, F, D, 0.0, workspace_a)
      BLAS.gemm!('N', 'N', 1.0, workspace_a, S, 0.0, FDS)

      transpose!(workspace_b, FDS)

      #== compute error vector ==#
      LinearAlgebra.BLAS.blascopy!(length(FDS), FDS, 1,
        workspace_a, 1)
      axpy!(-1.0, workspace_b, workspace_a)

      e_array_old = view(e_array,1:(ndiis-1))
      workspace_c[1] = deepcopy(workspace_a)
      e_array = vcat(workspace_c, e_array_old)
      F_array_old = view(F_array,1:(ndiis-1))
      workspace_c[1] = deepcopy(F)
      F_array = vcat(workspace_c, F_array_old)

      if iter > 1
        B_dim += 1
        B_dim = min(B_dim,ndiis)
        try
          DIIS(F, e_array, F_array, B_dim)
        catch
          println("Faulty DIIS!")
          B_dim = 2
          DIIS(F, e_array, F_array, B_dim)
        end
      end
    end

    #== dynamic damping of Fock matrix ==#
    x = ΔE >= 1.0 ? 1.0/log(50,50*ΔE) : 1.0
    LinearAlgebra.BLAS.axpby!(1.0-x, F_old, x, F)

    F_old .= F

    LinearAlgebra.BLAS.blascopy!(length(F), F, 1,
      F_old, 1)

    #== obtain new F,D,C matrices ==#
    LinearAlgebra.BLAS.blascopy!(length(D), D, 1,
      D_old, 1)

    E_elec, F_eval[:] = iteration(F, D, C, H, F_eval, F_evec, workspace_a,
      workspace_b, ortho, basis, iter, debug)

    #== check for convergence ==#
    #LinearAlgebra.BLAS.blascopy!(length(FDS), FDS, 1,
    #  workspace_a, 1)
    #axpy!(-1.0, workspace_b, workspace_a)

    ΔD .= D .- D_old
    D_rms = √(LinearAlgebra.BLAS.dot(length(ΔD), ΔD, 1, ΔD, 1))

    E = E_elec + E_nuc
    ΔE = E - E_old

    if MPI.Comm_rank(comm) == 0 && output >= 2
      #println(iter,"     ", E,"     ", ΔE,"     ", D_rms)
      @printf("%d      %.10f      %.10f      %.10f\n", iter, E, ΔE, D_rms)
    end

    iter_converged = Base.abs_float(ΔE) <= dele && D_rms <= rmsd
    iter += 1
    if iter > niter
      scf_converged = false
      break
    end

    #== if not converged, replace old D and E values for next iteration ==#
    E_old = E
  end

  #== build energy-weighted density matrix ==#
  #fill!(W, 0.0)
  LinearAlgebra.BLAS.scal!(length(W), 0.0, W, 1)

  nocc = basis.nels >> 1
  for i in 1:basis.norb, j in 1:basis.norb
    for iocc in 1:nocc
      W[i,j] += 2.0 * F_eval[iocc] * C[i, iocc] * C[j, iocc]
    end
  end

  return E
end

#=
"""
	 fock_build(F::Array{Float64}, D::Array{Float64}, tei::Array{Float64}, H::Array{Float64})
Summary
======
Perform Fock build step.

Arguments
======
F = Current iteration's Fock Matrix

D = Current iteration's Density Matrix

tei = Two-electron integral array

H = One-electron Hamiltonian Matrix
"""
=#

@inline function fock_build(F::Matrix{Float64},
  F_thread::Vector{Matrix{Float64}}, D::Matrix{Float64},
  H::Matrix{Float64}, basis::Basis,
  schwarz_bounds::Matrix{Float64}, Dsh::Matrix{Float64},
  eri_quartet_batch_thread::Vector{Vector{Float64}},
  jeri_engine_thread, iter::Int64,
  cutoff::Float64, debug::Bool, load::String)

  comm = MPI.COMM_WORLD

  #fill!(F,zero(Float64))
  LinearAlgebra.BLAS.scal!(length(F), 0.0, F, 1)
  #fill!.(F_thread,zero(Float64))
  for ithread_fock in F_thread
    LinearAlgebra.BLAS.scal!(length(ithread_fock), 0.0, ithread_fock, 1)
  end

  nsh = length(basis)
  nindices = (muladd(nsh,nsh,nsh)*(muladd(nsh,nsh,nsh) + 2)) >> 3

  batches_per_thread = nsh
  batch_size = ceil(Int,nindices/(MPI.Comm_size(comm)*
    Threads.nthreads()*batches_per_thread))

  #== use static task distribution for multirank runs if selected... ==#
  if MPI.Comm_size(comm) == 1  || load == "static"
    #== set up initial indices ==#
    stride = MPI.Comm_size(comm)*batch_size
    top_index = nindices - (MPI.Comm_rank(comm)*batch_size)

    #== execute kernel of calculation ==#
    @sync for ijkl_index in top_index:-stride:1
      Threads.@spawn begin
        thread = Threads.threadid()

        eri_quartet_batch_priv = eri_quartet_batch_thread[thread]
        jeri_tei_engine_priv = jeri_engine_thread[thread]
        F_priv = F_thread[thread]

        for ijkl in ijkl_index:-1:(max(1,ijkl_index-batch_size+1))
          fock_build_thread_kernel(F_priv, D,
            H, basis, eri_quartet_batch_priv,
            ijkl, jeri_tei_engine_priv,
            schwarz_bounds, Dsh,
            cutoff, debug)
        end
      end
    end

    #== reduce into Fock matrix ==#
    for (thread,ithread_fock) in enumerate(F_thread)
      if debug && MPI.Comm_rank(comm) == 0
        h5write("debug.h5","RHF/Iteration-$iter/F/Thread-$thread", ithread_fock)
      end

      axpy!(1.0, ithread_fock, F)
    end

  #== ..else use dynamic task distribution ==#
  elseif MPI.Comm_size(comm) > 1 && load == "dynamic"
    #== master rank ==#
    if MPI.Comm_rank(comm) == 0
      #== send out initial tasks to slaves ==#
      task = [ nindices ]
      initial_task = 1

      recv_mesg_master = [ 0 ]

      #println("Start sending out initial tasks")
      while initial_task < MPI.Comm_size(comm)
        for thread in 1:Threads.nthreads()
          #println("Sending task $task to rank $initial_task")
          sreq = MPI.Isend(task, initial_task, thread, comm)
          #println("Task $task sent to rank $initial_task")

          task[1] -= batch_size
        end
        initial_task += 1
      end
      #println("Done sending out intiial tasks")

      #== hand out quartets to slaves dynamically ==#
      #println("Start sending out rest of tasks")
      while task[1] > 0
        status = MPI.Probe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG,
          comm)
        #rreq = MPI.Recv!(recv_mesg_master, status.source, status.tag,
        #  comm)
        #println("Sending task $task to rank ", status.source)
        sreq = MPI.Isend(task, status.source, status.tag, comm)
        #println("Task $task sent to rank ", status.source)
        task[1] -= batch_size
      end
      #println("Done sending out rest of tasks")

      #== hand out ending signals once done ==#
      #println("Start sending out enders")
      for rank in 1:(MPI.Comm_size(comm)-1)
        for thread in 1:Threads.nthreads()
          sreq = MPI.Isend([ -1 ], rank, thread, comm)
        end
      end
      #println("Done sending out enders")
    #== worker ranks perform actual computations on quartets ==#
    elseif MPI.Comm_rank(comm) > 0
      #== create needed mutices ==#
      mutex_mpi_worker = Base.Threads.ReentrantLock()

      #== execute kernel ==#
      @sync for thread in 1:Threads.nthreads()
        Threads.@spawn begin
          #== initial set up ==#
          recv_mesg = [ 0 ]
          send_mesg = [ MPI.Comm_rank(comm) ]

          eri_quartet_batch_priv = eri_quartet_batch_thread[thread]
          jeri_tei_engine_priv = jeri_engine_thread[thread]
          F_priv = F_thread[thread]

          #== complete first task ==#
          lock(mutex_mpi_worker)
            status = MPI.Probe(0, $thread, comm)
            rreq = MPI.Recv!(recv_mesg, status.source, status.tag, comm)
            ijkl_index = recv_mesg[1]
          unlock(mutex_mpi_worker)

          for ijkl in ijkl_index:-1:(max(1,ijkl_index-batch_size+1))
            fock_build_thread_kernel(F_priv, D,
              H, basis, eri_quartet_batch_priv, #mutex,
              ijkl, jeri_tei_engine_priv,
              schwarz_bounds, Dsh,
              cutoff, debug)
          end

          #== complete rest of tasks ==#
          while ijkl_index >= 1
            lock(mutex_mpi_worker)
              status = MPI.Sendrecv!(send_mesg, 0, $thread, recv_mesg, 0,
                $thread, comm)
              ijkl_index = recv_mesg[1]
            unlock(mutex_mpi_worker)
            #println("Thread $thread ecieved task $ijkl_index from master")

            for ijkl in ijkl_index:-1:(max(1,ijkl_index-batch_size+1))
              fock_build_thread_kernel(F_priv, D,
                H, basis, eri_quartet_batch_priv, #mutex,
               ijkl, jeri_tei_engine_priv,
               schwarz_bounds, Dsh,
               cutoff, debug)
            end
          end
        end
      end

      #== reduce into Fock matrix ==#
      for ithread_fock in F_thread
        axpy!(1.0, ithread_fock, F)
      end
    end
  end
  MPI.Barrier(comm)

  for iorb in 1:basis.norb, jorb in 1:iorb
    if iorb != jorb
      F[iorb,jorb] /= 2.0
      F[jorb,iorb] = F[iorb,jorb]
    end
  end

  return F
end

@inline function fock_build_thread_kernel(F::Matrix{Float64}, D::Matrix{Float64},
  H::Matrix{Float64}, basis::Basis,
  eri_quartet_batch::Vector{Float64},
  ijkl_index::Int64,
  jeri_tei_engine, schwarz_bounds::Matrix{Float64},
  Dsh::Matrix{Float64}, cutoff::Float64, debug::Bool)

  comm=MPI.COMM_WORLD

  #== determine shells==#
  bra_pair = decompose(ijkl_index)
  ket_pair = ijkl_index - triangular_index(bra_pair)

  #quartet.bra = basis.shpair_ordering[bra_pair]
  #quartet.ket = basis.shpair_ordering[ket_pair]

  #ish = μsh.shell_id
  #jsh = νsh.shell_id
  #ksh = λsh.shell_id
  #lsh = σsh.shell_id

  ish = decompose(bra_pair)
  jsh = bra_pair - triangular_index(ish)

  ksh = decompose(ket_pair)
  lsh = ket_pair - triangular_index(ksh)

  μsh = basis[ish]
  νsh = basis[jsh]
  λsh = basis[ksh]
  σsh = basis[lsh]

  #icls = unsafe_string(μsh.class)
  #jcls = unsafe_string(νsh.class)
  #kcls = unsafe_string(λsh.class)
  #lcls = unsafe_string(σsh.class)

  #println("QUARTET($ish, $jsh, $ksh, $lsh) -> ($icls $jcls | $kcls $lcls)")

  #== Cauchy-Schwarz screening ==#
  bound = schwarz_bounds[ish, jsh]*schwarz_bounds[ksh, lsh]

  #println("SCHWARZ($ish, $jsh): $(schwarz_bounds[ish, jsh])")

  dijmax = 4.0*Dsh[ish, jsh]
  dklmax = 4.0*Dsh[ksh, lsh]

  dikmax = Dsh[ish, ksh]
  dilmax = Dsh[ish, lsh]
  djkmax = Dsh[jsh, ksh]
  djlmax = Dsh[jsh, lsh]

  maxden = max(dijmax, dklmax, dikmax, dilmax, djkmax, djlmax)
  bound *= maxden

  #== fock build for significant shell quartets ==#
  if Base.abs_float(bound) >= cutoff
    #= set up some variables =#
    nμ = μsh.nbas
    nν = νsh.nbas
    nλ = λsh.nbas
    nσ = σsh.nbas

    #== compute electron repulsion integrals ==#
    screened = compute_eris(ish, jsh, ksh, lsh, bra_pair, ket_pair, nμ, nν,
      nλ, nσ, eri_quartet_batch, jeri_tei_engine)

    if !screened
      #= axial normalization =#
      axial_normalization_factor(eri_quartet_batch, μsh, νsh, λsh, σsh,
        nμ, nν, nλ, nσ)

      #== contract ERIs into Fock matrix ==#
      contract_eris(F, D, eri_quartet_batch, ish, jsh, ksh, lsh,
        μsh, νsh, λsh, σsh, cutoff, debug)
    end
  end
  #if debug println("END TWO-ELECTRON INTEGRALS") end
end

@inline function compute_eris(ish::Int64, jsh::Int64, ksh::Int64, lsh::Int64,
  bra_pair::Int64, ket_pair::Int64,
  nμ::Int64, nν::Int64,
  nλ::Int64, nσ::Int64,
  eri_quartet_batch::Vector{Float64},
  jeri_tei_engine)

  #= actually compute integrals =#
  return JERI.compute_eri_block(jeri_tei_engine, eri_quartet_batch,
    ish, jsh, ksh, lsh, bra_pair, ket_pair, nμ*nν, nλ*nσ)

  #=
  if am[1] == 3 || am[2] == 3 || am[3] == 3 || am[4] == 3
    for idx in 1:nμ*nν*nλ*nσ
    #for idx in 1:1296
      eri = eri_quartet_batch[idx]
      println("QUARTET($ish, $jsh, $ksh, $lsh): $eri")
    end
  end
  =#
end

@inline function contract_eris(F_priv::Matrix{Float64}, D::Matrix{Float64},
  eri_batch::Vector{Float64}, ish::Int64, jsh::Int64,
  ksh::Int64, lsh::Int64,
  μsh::JCModules.Shell, νsh::JCModules.Shell,
  λsh::JCModules.Shell, σsh::JCModules.Shell,
  cutoff::Float64, debug::Bool)

  #norb = size(D,1)

  #ish = μsh.shell_id
  #jsh = νsh.shell_id
  #ksh = λsh.shell_id
  #lsh = σsh.shell_id

  pμ = μsh.pos
  nμ = μsh.nbas

  pν = νsh.pos
  nν = νsh.nbas

  pλ = λsh.pos
  nλ = λsh.nbas

  pσ = σsh.pos
  nσ = σsh.nbas

  #amμ = μsh.am
  #amν = νsh.am
  #amλ = λsh.am
  #amσ = σsh.am
  #am = [ amμ, amν, amλ, amσ ]

  for μμ::Int64 in pμ:(pμ+nμ-1)
    νmax = ish == jsh ? μμ : (pν+nν-1)
    for νν::Int64 in pν:νmax
      for λλ::Int64 in pλ:(pλ+nλ-1)
        σmax = ksh == lsh ? λλ : (pσ+nσ-1)
        for σσ::Int64 in pσ:σmax
          #if debug
            #if do_continue_print print("$μμ, $νν, $λλ, $σσ => ") end
          #end
          μνλσ = 1 + (σσ-pσ) + nσ*(λλ-pλ) + nσ*nλ*(νν-pν) + nσ*nλ*nν*(μμ-pμ)
          jlswap = (μμ < λλ) || (μμ == λλ && νν < σσ)

          eri = eri_batch[μνλσ]

          if Base.abs_float(eri) < cutoff
            #if do_continue_print println("CONTINUE SCREEN") end
            continue
          elseif jlswap && ish == ksh && jsh == lsh
            continue
          end

          μ = max(μμ, λλ)
          ν = !jlswap ? νν : σσ
          λ = min(μμ, λλ)
          σ = !jlswap ? σσ : νν

          #println("QUARTET($ish, $jsh, $ksh, $lsh): $eri")
          #println("ERI($μ, $ν, $λ, $σ) = $eri")

          eri *= (μ == ν) ? 0.5 : 1.0
          eri *= (λ == σ) ? 0.5 : 1.0
          eri *= ((μ == λ) && (ν == σ)) ? 0.5 : 1.0

          F_priv[λ,σ] += 4.0 * D[μ,ν] * eri
          F_priv[μ,ν] += 4.0 * D[λ,σ] * eri
          F_priv[μ,λ] -= D[ν,σ] * eri
          F_priv[μ,σ] -= D[ν,λ] * eri
          F_priv[max(ν,λ), min(ν,λ)] -= D[μ,σ] * eri
          F_priv[max(ν,σ), min(ν,σ)] -= D[μ,λ] * eri
        end
      end
    end
  end
end

#=
"""
	 iteration(F::Matrix{Float64}, D::Matrix{Float64}, H::Matrix{Float64}, ortho::Matrix{Float64})
Summary
======
Perform single SCF iteration.

Arguments
======
D = Current iteration's Density Matrix

H = One-electron Hamiltonian Matrix

ortho = Symmetric Orthogonalization Matrix
"""
=#
function iteration(F_μν::Matrix{Float64}, D::Matrix{Float64},
  C::Matrix{Float64}, H::Matrix{Float64}, F_eval::Vector{Float64},
  F_evec::Matrix{Float64}, workspace_a::Matrix{Float64},
  workspace_b::Matrix{Float64}, ortho::Matrix{Float64},
  basis::Basis, iter::Int, debug::Bool)

  comm=MPI.COMM_WORLD

  transpose!(workspace_b, LinearAlgebra.Hermitian(ortho))

  #== obtain new orbital coefficients ==#
  BLAS.symm!('L', 'U', 1.0, workspace_b, F_μν, 0.0, workspace_a)
  BLAS.gemm!('N', 'N', 1.0, workspace_a, ortho, 0.0, workspace_b)

  F_eval[:], F_evec[:,:] = eigen!(LinearAlgebra.Hermitian(workspace_b))

  #@views F_evec .= F_evec[:,sortperm(F_eval)] #sort evecs according to sorted evals

  if debug && MPI.Comm_rank(comm) == 0
    h5write("debug.h5","RHF/Iteration-$iter/F_evec", F_evec)
  end

  #C .= ortho*F_evec
  BLAS.symm!('L', 'U', 1.0, ortho, F_evec, 0.0, C)

  if debug && MPI.Comm_rank(comm) == 0
    h5write("debug.h5","RHF/Iteration-$iter/C", C)
  end

  #== build new density matrix ==#
  nocc = basis.nels >> 1
  norb = basis.norb

  #fill!(D, 0.0)
  for i in 1:basis.norb, j in 1:basis.norb
    D[i,j] = 2.0*BLAS.dot(nocc,pointer(C,i),norb,pointer(C,j),norb)
  end

  #== compute new SCF energy ==#
  #EHF1 = LinearAlgebra.dot(D, F_μν)
  #EHF2 = LinearAlgebra.dot(D, H)
  EHF1 = LinearAlgebra.BLAS.dot(length(D), D, 1, F_μν, 1)
  EHF2 = LinearAlgebra.BLAS.dot(length(D), D, 1, H, 1)
  E_elec = (EHF1 + EHF2)/2.0

  if debug && MPI.Comm_rank(comm) == 0
    h5write("debug.h5","RHF/Iteration-$iter/D", D)
    h5write("debug.h5","RHF/Iteration-$iter/E/EHF1", EHF1)
    h5write("debug.h5","RHF/Iteration-$iter/E/EHF2", EHF2)
    h5write("debug.h5","RHF/Iteration-$iter/E/EHF", E_elec)
  end

  return E_elec, F_eval
end


text:
module FunctionIndices

using MappedArrays

export FI, NotIndex, not, notin

"""
    AbstractFunctionIndex

Supertype of all function index types.
"""
abstract type AbstractFunctionIndex end

"""
    FunctionIndices.to_index(::Type{T<:AbstractArray}, A, ind, i)

Convert a `AbstractFunctionIndex` `i` to a array index of type `T` for `A` with `ind`.
By default, `to_index(::AbstractArray, ind, i)` will return a
`Base.LogicalIndex{Bool, ReadonlyMappedArray{Bool...}}`.
"""
@inline to_index(::Type{T}, A, ind, i) where {T<:AbstractArray} = _to_logic_index(
    # use mappedarray instead of map for less allocations and more information
    IndexStyle(A),
    A,
    mappedarray(to_function(i), ind)::ReadonlyMappedArray{Bool},
) # no type assert here, because this methods will accept any T even T is not a LogicalIndex

_to_logic_index(::IndexLinear, A, i) = Base.LogicalIndex{Int}(i)
_to_logic_index(::IndexStyle, A, i) = Base.to_index(A, i)

@inline Base.to_indices(A, inds, I::Tuple{AbstractFunctionIndex,Vararg{Any}}) = (
    to_index(indextype(A, I[1]), A, _maybefirst(inds), I[1]),
    to_indices(A, Base._maybetail(inds), Base.tail(I))...,
)

_maybefirst(::Tuple{}) = Base.OneTo(1)
_maybefirst(inds::Tuple) = inds[1]

"""
    indextype([A::AbstractArray,] I::AbstractFunctionIndex)
    indextype([::Type{TA},] ::Type{TI})

Determine the index type which the function index type `TI`
will be converted to by `Base.to_indices` for array type `TA`.
By default, it's `AbstractArray`.

!!! note

    If you define a methods `indextype(::Type{TA}, ::Type{TI}) = T`,
    while `to_index(::Type{T}, inds, I::TI)` is not defined,
    the `to_indices` will not convert a index of type `TI` to a `T`,
    but a `Base.LogicalIndex`, the default return type of `to_index`.
"""
indextype(A::AbstractArray, I::AbstractFunctionIndex) = indextype(typeof(A), typeof(I))
indextype(::Type{<:AbstractArray}, ::Type{<:AbstractFunctionIndex}) = AbstractArray

"""
    to_function(I::AbstractFunctionIndex)

Converts a `AbstractFunctionIndex` to a `Function`.
By default, `to_function(I::AbstractNotIndex)` returns `notin(parent(I))`.
"""
to_function

"""
    FunctionIndex{F} <: AbstractFunctionIndex
    FI

A implementation of function index, which make indexing with function possible.
1-d indexing `A[FI(f)]` is equivalent to `A[map(f, begin:end)]`,
multi-dimensional indexing `A[FI(f1), ..., FI(fn)]` is equivalent to
`A[map(FI(f1), axes(A, 1)), ..., map(FI(fn), axes(A, n))]`.
"""
struct FunctionIndex{F} <: AbstractFunctionIndex
    f::F
end
const FI = FunctionIndex

to_function(ind::FunctionIndex) = ind.f

"""
    AbstractNotIndex{T}

Supertype of all not types which create by [`not`](@ref).
"""
abstract type AbstractNotIndex{T} <: AbstractFunctionIndex end

"""
    notin(item, collection)
    notin(collection)

The same as `!in`, used in the default `to_function` for `AbstractNotIndex`.
"""
notin(item, collection) = !(item in collection)
notin(collect) = Base.Fix2(notin, collect)

to_function(I::AbstractNotIndex) = notin(parent(I))

"""
    NotIndex{T} <: AbstractNotIndex{T}

The default implementation of [`not`](@ref).
There are some optimization for `NotIndex(x)` comparing to `FI(!in(x))`.
For large arrays, this implementation may be faster.
but for small arrays this implementation may be slower.
See [performance comparing](@ref performance) for more details.
"""
struct NotIndex{T} <: AbstractNotIndex{T}
    parent::T
end
Base.parent(I::NotIndex) = I.parent

# convert to Vector{Int} by default for NotIndex
indextype(::Type{<:AbstractArray}, ::Type{<:NotIndex}) = Vector{Int}

"""
    not(x)

Create a `NotIndex` with given `x`, which is similar to `Not` of
[`InvertedIndices`](https://github.com/JuliaData/InvertedIndices.jl).
In most cases, `not` is much faster than `Not`.
See [performance comparing](@ref performance) for more details.

# Main differences between `not` and `Not`

For `CartesianIndex`, `A[not(CartesianIndex(i, j,..., n))]` is equivalent to
`A[not(i), not(j), ..., not(n)]` and will return a array with same dimension,
but `A[Not(CartesianIndex(i, j,..., n))]` will convert the `CartesianIndex` to
a linear index and return a vector:

```jldoctest
julia> A = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 0  3  6   9
 1  4  7  10
 2  5  8  11

julia> A[not(CartesianIndex(1, 2))] # equivalent to A[not(1), not(2)]
2×3 Matrix{Int64}:
 1  7  10
 2  8  11

julia> A[Not(CartesianIndex(1, 2))] # equivalent to A[Not(3)]
11-element Vector{Int64}:
  0
  1
  2
  4
  5
  ⋮
  8
  9
 10
 11
```

Besides, for index like `A[4, 5]` which is out of bounds,
`A[not(4), not(5)]` is equivalent to `A[:, :]`,
because inbounds indices are always not equal to the given value,
while `A[Not[4], Not(5)]` causes an error.
"""
not(x) = NotIndex(x)
not(x::Integer) = NotIndex(Int(x))
not(x::T, xs::T...) where {T} = not((x, xs...))

# CartesianIndices and CartesianIndex should be difference
const NotCartesian{N} =
    Union{AbstractNotIndex{<:CartesianIndex{N}},AbstractNotIndex{<:CartesianIndices{N}}}

# NotCartesian{0} will be ignored by to_indices
# A[not(CartesianIndex()), 1] == A[1] like A[CartesianIndex(), 1] == A[1],
@inline Base.to_indices(A, inds, I::Tuple{<:NotCartesian{0},Vararg{Any}}) =
    to_indices(A, inds, Base.tail(I))
# NotCartesian{N} will be convert to as Vararg{AbstractNotIndex{Int}, N}
@inline Base.to_indices(A, inds, I::Tuple{<:NotCartesian,Vararg{Any}}) =
    to_indices(A, inds, (_to_linear_nots(I[1])..., Base.tail(I)...))

# methods override Base.to_indices(A, ::Type{Any})
# to avoid regard NotCartesian as a LinearIndex
@inline Base.to_indices(A, I::Tuple{<:NotCartesian}) = to_indices(A, _to_linear_nots(I[1]))

# convert the NotCartesian to normal nots
_to_linear_nots(I::NotCartesian) = map(_type_unionall(I), _getIS(parent(I)))
# to strip parameter from AbstractNotIndex, T.name.wrapper is hacky (T is always a DataType)
_type_unionall(::T) where {T} = T.name.wrapper

_getIS(I::CartesianIndex) = I.I
_getIS(I::CartesianIndices) = I.indices

# utils
## SetArray as a replace of Set for more information about parent
struct SetArray{T,N,TA<:AbstractArray{T,N}} <: AbstractArray{T,N}
    A::TA
    s::Set{T}
    SetArray(A::AbstractArray{T,N}) where {T,N} = new{T,N,typeof(A)}(A, Set{T}(A))
end
# no more methods needed for SetArray
Base.parent(A::SetArray) = A.A
Base.in(x, A::SetArray) = in(x, A.s)
## Tuple as a Vector to indexing
struct TupleVector{T,P<:Tuple} <: AbstractVector{T}
    parent::P
    TupleVector(tp::Tuple{Vararg{T}}) where {T} = new{T,typeof(tp)}(tp)
end
Base.parent(V::TupleVector) = V.parent
Base.size(V::TupleVector) = (length(parent(V)),)
Base.@propagate_inbounds Base.getindex(V::TupleVector{T,N}, i::Integer) where {T,N} =
    parent(V)[i]

# Optimization for some special cases
## Optimization for any indextype(A)
### For not(::Colon), (not(::Slice) is only optimized for indextype{A} == Vector{Int})
to_index(::Type{<:AbstractArray}, A, ind, ::AbstractNotIndex{<:Colon}) = []
### For not(::AbstractArray{Bool})
to_index(::Type{<:AbstractArray}, A, ind, I::AbstractNotIndex{<:AbstractArray{Bool}}) =
    _to_logic_index(IndexStyle(A), A, mappedarray(!, parent(I)))
### For not(::AbstractArray), only test if in like a Set for AbstractArray
to_function(I::NotIndex{<:AbstractArray}) = notin(SetArray(parent(I)))
### For converted function index
to_index(
    ::Type{<:AbstractArray},
    A,
    ind,
    I::AbstractNotIndex{<:Base.LogicalIndex{<:Any,<:ReadonlyMappedArray{Bool}}},
) = (rma = parent(I).mask; _not_mapped(A, rma.f, ind))
_not_mapped(A, f, ind) = _to_logic_index(IndexStyle(A), A, mappedarray(!f, ind))
_not_mapped(A, f::Base.Fix2{typeof(notin),<:Tuple}, ::Any) = TupleVector(f.x)
_not_mapped(A, f::Base.Fix2{typeof(notin),<:SetArray}, ::Any) = parent(f.x)

## Optimize only for indextype(I) == Vector{Int}
const TVInt = Type{Vector{Int}}
function to_index(::TVInt, A, ind::AbstractUnitRange{<:Integer}, I::AbstractNotIndex{<:Integer})
    n = parent(I)
    return n in ind ? [first(ind):(n-1); (n+1):last(ind)] : collect(ind)
end

function to_index(
    ::TVInt,
    A,
    ind::AbstractUnitRange{<:Integer},
    I::AbstractNotIndex{<:AbstractUnitRange{<:Integer}},
)
    r = _check_index(ind, parent(I))
    return isempty(r) ? collect(ind) : [first(ind):(first(r)-1); (last(r)+1):last(ind)]
end
# to_index for NotIndex{<:StepRange} is suboptimal for small size arrays
function to_index(
    ::TVInt,
    A,
    ind::AbstractUnitRange{<:Integer},
    I::AbstractNotIndex{<:StepRange{<:Integer}},
)
    r = _check_index(ind, parent(I))
    isempty(r) && return collect(ind)
    # collect ind for type stable
    step(r) == 1 && return [first(ind):first(r)-1; last(r)+1:last(ind)]
    ret = Vector{Int}(undef, length(ind) - length(r))
    copyto!(ret, first(ind):first(r)-1)
    offset = first(r) - first(ind) + 1
    @inbounds for i in 1:(last(r)-first(r))
        if i % step(r) != 0
            ret[offset] = i + first(r)
            offset += 1
        end
    end
    copyto!(ret, offset, last(r)+1:last(ind))
    return ret
end
function to_index(
    ::TVInt,
    A,
    ind::AbstractUnitRange{<:Integer},
    I::AbstractNotIndex{<:Base.Slice},
)
    if parent(I) == ind # if slice is the same as the range, return a empty array
        return Int[]::Vector{Int}
    else # if slice is not equal to ind, invoke the method for unit range
        return invoke(
            to_index,
            Tuple{TVInt,typeof(A),typeof(ind),AbstractNotIndex{<:AbstractUnitRange{Int}}},
            Vector{Int},
            A,
            ind,
            I,
        )::Vector{Int}
    end
end

# process index like (1:2)[not(-1:3)] to (1:2)[not(1:2)]
# and (1:2)[not(2:-1:1)] to (1:2)[not(1:2)]
function _check_index(ind::AbstractUnitRange{<:Integer}, r::StepRange{<:Integer})
    if step(r) > 0
        return max(first(ind), first(r)):step(r):min(last(ind), last(r))
    else
        return max(first(ind), last(r)):-step(r):min(last(ind), first(r))
    end
end
_check_index(ind::AbstractUnitRange{<:Integer}, r::AbstractUnitRange{<:Integer}) =
    max(first(ind), first(r)):min(last(ind), last(r))

end # module

ext: jl


text:
module LibCURL

using LibCURL_jll
using MozillaCACerts_jll

const time_t = Int
const size_t = Csize_t
const curl_off_t = Int64
const fd_set = Union{}
const socklen_t = Int32

cacert = ""

export Mime_ext

function __init__()
    # Note: `MozillaCACerts_jll.cacert` is filled by `__init__` which requires LibCURL's
    # copy to also be filled in during initialization. Doing this ensures compatibility
    # with building system images.
    global cacert = MozillaCACerts_jll.cacert
end

include("lC_exports_h.jl")
include("lC_common_h.jl")
include("lC_curl_h.jl")

curl_easy_setopt(handle, opt, ptrval::Array) = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{Cvoid}), handle, opt, pointer(ptrval))
curl_easy_setopt(handle, opt, ptrval::Integer) = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Clong), handle, opt, ptrval)
curl_easy_setopt(handle, opt, ptrval::Ptr{T}) where {T} = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{T}), handle, opt, ptrval)
curl_easy_setopt(handle, opt, ptrval::AbstractString) = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{UInt8}), handle, opt, ptrval)

curl_multi_setopt(handle, opt, ptrval::Array) = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Ptr{Cvoid}), handle, opt, pointer(ptrval))
curl_multi_setopt(handle, opt, ptrval::Integer) = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Clong), handle, opt, ptrval)
curl_multi_setopt(handle, opt, ptrval::Ptr{T}) where {T} = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Ptr{T}), handle, opt, ptrval)
curl_multi_setopt(handle, opt, ptrval::AbstractString) = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Ptr{UInt8}), handle, opt, ptrval)

curl_easy_getinfo(handle, opt, ptrval::Array) = ccall((:curl_easy_getinfo, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{Cvoid}), handle, opt, pointer(ptrval))
curl_easy_getinfo(handle, opt, ptrval::AbstractString) = ccall((:curl_easy_getinfo, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{UInt8}), handle, opt, ptrval)

include("lC_defines_h.jl")

include("Mime_ext.jl")

end


text:
@testset "PoissonLikelihood" begin
    for args in ((), (exp,), (ExpLink(),))
        @test PoissonLikelihood(args...) isa PoissonLikelihood{ExpLink}
    end

    for lik in (PoissonLikelihood(), PoissonLikelihood(log1pexp))
        test_interface(lik, SqExponentialKernel(), rand(10))
    end
end


text:
const i10 = parse(Int, ARGS[1])

using SmartTransitionSim
using StatsBase
using Random

#Creating MapData object
const mapfile = "SanFranciscoFiltered.osm"
const datapath = "./example/";
const RoadSet = 5
const map_data = get_map_data(datapath, mapfile, use_cache=false; road_levels = Set(1:RoadSet));

#Dictionary for memorizing used k-shortest paths
const KDict = Dict{Tuple{Int,Int},Array{Vector{Int}}}()

#Create parameters grid
const density_factor = 5.0
const Params = (reps = 1:3, α = 0.05:0.05:1.0, N = 3000:500:7500, U = 50:50:300, T = [0.1,1.0,10.0], k = [1,2,3,4,5], AvgDelay = 0.0)
const ParametersGridProd = collect(Base.Iterators.product(Params.reps,Params.α, Params.N, Params.U, Params.T, Params.k, Params.AvgDelay))
const ParametersGrid = [ParametersGridProd[i] for i in 1:length(ParametersGridProd)]

#Defining starting and ending area
const Start = [Rect((37.795,-122.42),(37.811,-122.3942))]
const End = [Rect((37.7805,-122.485),(37.79,-122.45))]

#Create Agents Pool
T = 0.5
kmode = mode([ParametersGrid[i][6] for i in (i10*40-39):(i10*40)])
Delay = ParametersGrid[1][7]

Random.seed!(1);

Agents1 = generate_agents(map_data, 5000, Start, End, 0.0, kmode, T, Delay, KDict)
Agents2 = generate_agents(map_data, 5000, Start, End, 0.0, kmode, T, Delay, KDict)
const Agents = [Agents1; Agents2]
#Output file header
println("ID,rep,α,N,U,T,k,AvgDelay,runtime,overall_time_red,smart_time_red,regular_time_red,avg_time_base,avg_time_all_V2I,avg_time_smart_V2I,avg_time_regular_V2I")

for i in (i10*40-39):(i10*40)
    Random.seed!(1000);
    run_parameter_analysis(i,ParametersGrid,Agents,density_factor,KDict,map_data)
    flush(stdout)
end


text:
# prompt.jl

import MIRT # prompt
using Test: @test

@test MIRT.prompt(:state) isa Symbol


text:
export image_to_photons
export precompute_distances_and_angles, calculate_correlations_in_image_integral, symmetrize_image, calculate_correlations_in_image_using_single_photons

"""Transforms a sparse block-image into a list of photons that can be used to calculate correlations"""
function image_to_photons(img::Matrix{UInt32}, qmax::Float64=1.0)
    photons = Vector{Float64}[]
    xsize, ysize = Base.size(img)
    xhalf = floor(xsize/2.0)
    yhalf = floor(ysize/2.0)
    dq = qmax / xhalf
    #     println("Sorting $(sum(abs, img)) photons.")
    for x=1:xsize
        for y=1:ysize
            if img[x,y] > 0 && (x-xhalf,y-yhalf) != (0,0)
                for p in Base.Iterators.repeated(dq*[x-xhalf,y-yhalf,0], 1)#Int64(img[x,y])
                    push!(photons, p)
                end
            end
        end
    end
    return photons
end

############################################################################
## EXPERIMENTAL DATA
############################################################################

function precompute_distances_and_angles(image_width::Int64, N::Int64)
    distances = zeros(Int64, image_width, image_width)
    angles = zeros(Int64, image_width, image_width, image_width, image_width)

    range = 1:image_width
    center = ceil(Float64[image_width/2.0,image_width/2.0, 0.0])
    da = pi/N
    for x1 = range
        for y1 = range
            p1 = Float64[x1,y1, 0.0] - center
            @inbounds distances[x1,y1] = round(Int64,norm(p1))
            for x2 = range
                for y2 = range
                    p2 = Float64[x2,y2, 0.0] - center

                    #calculate angles
                    alpha = angle_between(p1,p2)
                    @inbounds angles[x1, y1, x2, y2] = clamp(floor(Int64, alpha/da) + 1, 1, 2*N)
                end
            end
        end
    end
    return distances, angles
end

"""Symmetrizes any image"""
function symmetrize_image(img::Matrix{Float64})
    # img_sym = deepcopy(img)
    # center = [ceil(Base.size(img)[1]/2.0), ceil(Base.size(img)[2]/2.0)]
    # for x=1:Base.size(img)[1]
    #     for y=1:x
    #         pos = [x,y]
    #         mpos = round(Int64,-1.0*([x,y] - center)+center)
    #         m = 0.5*(img_sym[x,y] + img_sym[mpos[1],mpos[2]])
    #         img_sym[x,y] = m
    #         img_sym[mpos[1],mpos[2]] = m
    #     end
    # end
    # return img_sym
    return 0.5*(img+rot180(img))
end

"""Calculates the two- and three-photon correlation from dense pixelized images"""
function calculate_correlations_in_image_integral(image_list::Array{Array{Float64,2},1}, K2::Int64, K3::Int64, N::Int64=32, filename::String="histo.dat", symmetrize::Bool=false)

    da = pi/N
    (image_width,sy) = Base.size(image_list[1])
    center = ceil(Float64[image_width/2.0,image_width/2.0])
    rotations = [[cos(phi) sin(phi); -sin(phi) cos(phi)] for phi = 0:da:2*pi]
    up = Float64[0.0, 1.0]

    c1_full = zeros(Float64, K2)
    c2_full = zeros(Float64, N, K2, K2)
    c3_full = zeros(Float64, N, 2*N, K3, K3, K3)
    number_analyzed_images = 0

    for j=1:ceil(Int64, length(image_list)/nworkers())
        println("Processing batch $j")
        #Processing next batch of nworkers() images
        c1_part,c2_part,c3_part = @sync @parallel ( (a,b) -> (a[1]+b[1], a[2]+b[2], a[3]+b[3])) for i=((j-1)*nworkers()+1):clamp(j*nworkers()+1, 1, length(image_list))
            image = image_list[i]
            if symmetrize == true
                image = symmetrize_image(image)
            end
            println("Processing image #$(i)")
            c1_local = zeros(Float64, K2)
            c2_local = zeros(Float64, N, K2, K2)
            c3_local = zeros(Float64, N, 2*N, K3, K3, K3)

            for k1 = 1:K2
                write(STDOUT,"$k1 ")
                p1 = k1*up

                sum1 = 0.0
                for rot in rotations
                    rp1 = round(Int64,rot*p1 + center)
                    sum1 += image[rp1[1], rp1[2]]
                end
                c1_local[k1] = sum1 * 1/(k1)

                for k2 = 1:k1
                    for alpha = 1:2*N
                        a = Float64(da*alpha-da/2)
                        p2 = k2*[cos(a) sin(a); -sin(a) cos(a)]*up
                        ais = alpha
                        if alpha > N ais = 2*N-ais+1 end

                        sum2 = 0.0
                        for rot in rotations
                            rp1 = round(Int64,rot*p1 + center)
                            rp2 = round(Int64,rot*p2 + center)
                            sum2 += image[rp1[1], rp1[2]]*image[rp2[1], rp2[2]]
                        end
                        c2_local[ais,k2,k1] = sum2 * 1/(k1*k2) * doubletFactor(k1,k2)

                        for k3=1:k2
                            for beta=1:2*N

                                if k1 <= K3 && k2 <= K3 && k3 <= K3
                                    b = Float64(da*beta-da/2)
                                    p3 = k3*[cos(b) sin(b); -sin(b) cos(b)]*up
                                    bis = beta
                                    if alpha > N bis = 2*N-bis+1 end

                                    sum3 = 0.0
                                    for rot in rotations
                                        rp1 = round(Int64,rot*p1 + center)
                                        rp2 = round(Int64,rot*p2 + center)
                                        rp3 = round(Int64,rot*p3 + center)
                                        sum3 += image[rp1[1], rp1[2]]*image[rp2[1], rp2[2]]*image[rp3[1], rp3[2]]
                                    end
                                    c3_local[ais,bis,k3,k2,k1] = sum3 * 1/(k1*k2*k3) * tripletFactor(k1,k2,k3)
                                end
                            end
                        end
                    end
                end
            end
            (c1_local, c2_local, c3_local)
        end
        c1_full += c1_part
        c2_full += c2_part
        c3_full += c3_part
        number_analyzed_images += nworkers()

        serializeToFile(filename, ( Dict("num_pictures"=>number_analyzed_images, "num_incident_photons"=>0, "qcut"=>1.0, "K2"=>K2, "K3"=>K3, "N"=>N, "dq"=>0.1), c2_full, c3_full, c1_full))
        flush(STDOUT)
    end
end

"""Calculates the two- and three-photon correlation from dense pixelized images"""
function calculate_correlations_in_image_using_single_photons(image_list::Array{Array{Float64,2},1}, K2::Int64, K3::Int64, N::Int64, filename::String, overall_maximum::Float64, photons_per_image::Int64, symmetrize::Bool=false)

    da = pi/N
    (image_width,sy) = Base.size(image_list[1])
    center = ceil.(Float64[image_width/2.0,image_width/2.0])
    rotations = [[cos(phi) sin(phi); -sin(phi) cos(phi)] for phi = 0:da:2*pi]
    up = Float64[0.0, 1.0]

    c1_full = zeros(Float64, K2)
    c2_full = zeros(Float64, N, K2, K2)
    c3_full = zeros(Float64, N, 2*N, K3, K3, K3)
    number_analyzed_images = 0

    for j=1:ceil(Int64, length(image_list)/nworkers())
        println("Processing batch $j")

        #Processing next batch of nworkers() images
        c1_part,c2_part,c3_part = @sync @parallel ( (a,b) -> (a[1]+b[1], a[2]+b[2], a[3]+b[3])) for i=((j-1)*nworkers()+1):j*nworkers()
            image = image_list[i]
            if symmetrize == true
                image = symmetrize_image(image)
            end
            println("Processing image #$(i)")
            c1_local = zeros(Float64, K2)
            c2_local = zeros(Float64, N, K2, K2)
            c3_local = zeros(Float64, N, 2*N, K3, K3, K3)

            #Calculating maximum intensity and expected photon counts
            integral = sum(abs,image)
            ppi = integral/overall_maximum * photons_per_image

            image_volume = ImageVolume(image, float(K2))

            photon_list,_ = pointsPerOrientation(image_volume, float(K2), float(K2), photons_per_image, rot=eye(3), lambda=0.0, beamstop_width=0.0, print_warning=false)

            histogramCorrelationsInPicture_alltoall(photon_list, c1_local, c2_local, c3_local, 1.0, N, K2, K3, 0.0)

            (c1_local, c2_local, c3_local)
        end
        c1_full += c1_part
        c2_full += c2_part
        c3_full += c3_part
        number_analyzed_images += nworkers()

        serializeToFile(filename, ( Dict("num_pictures"=>number_analyzed_images, "num_incident_photons"=>0, "qcut"=>1.0, "K2"=>K2, "K3"=>K3, "N"=>N, "dq"=>0.1), c2_full, c3_full, c1_full))
        flush(STDOUT)
    end
end


text:
module Gradientflow_module
    import ..AbstractGaugefields_module:AbstractGaugefields,Abstractfields,initialize_TA_Gaugefields,add_force!,
    exp_aF_U!,clear_U!,add_U!,evaluate_wilson_loops!,exptU!,
    Traceless_antihermitian_add!,set_wing_U!,Traceless_antihermitian,evaluate_gaugelinks!,
    construct_Λmatrix_forSTOUT!,Traceless_antihermitian!,shift_U
    import ..GaugeAction_module:GaugeAction,get_temporary_gaugefields,calc_dSdUμ!
    import ..Abstractsmearing_module:Abstractsmearing
    import Wilsonloop:make_loops_fromname,Wilsonline
    using LinearAlgebra
    import Wilsonloop:LinearAlgebra.adjoint

struct Gradientflow_general{Dim,TA,T} <: Abstractsmearing
    Nflow::Int64
    eps::Float64
    gaugeaction::GaugeAction{Dim,T}
    _temporal_TA_field::Array{TA,1}
    _temporal_G_field::Array{T,1}
    _temporal_U_field::Array{Array{T,1},1}

    function Gradientflow_general(U::Array{<:AbstractGaugefields{NC,Dim},1},linknames,linkvalues;Nflow = 1,eps = 0.01) where {NC,Dim }
        @assert length(linknames) == length(linkvalues)
        numlinks = length(linknames)
        links = Vector{Vector{Wilsonline{Dim}}}(undef,numlinks)
        for i=1:numlinks
            links[i] = make_loops_fromname(linknames[i],Dim=Dim)
        end

        return Gradientflow_general(U,links,linkvalues,Nflow = Nflow,eps = eps)
    end

    function Gradientflow_general(U::Array{<:AbstractGaugefields{NC,Dim},1},links::Vector{Vector{Wilsonline{Dim}}},linkvalues;Nflow = 1,eps = 0.01) where {NC,Dim }
        F0 = initialize_TA_Gaugefields(U)
        Ftemps = Array{typeof(F0),1}(undef,4)
        Ftemps[1] = F0
        for i=2:4
            Ftemps[i] = initialize_TA_Gaugefields(U)
        end
        T = eltype(U)
        Utemps = Array{Array{T,1},1}(undef,2)
        for i=1:2
            Utemps[i] = similar(U)
        end

        tempG = Array{T,1}(undef,3)
        for i=1:3
            tempG[i] = similar(U[1])
        end

        gaugeaction = GaugeAction(U)
        @assert length(links) == length(linkvalues)
        numlinks = length(links)
        for i=1:numlinks
            loop = links[i]#make_loops_fromname(linknames[i],Dim=Dim)
            factor = linkvalues[i]
            if typeof(factor) <: Real
                append!(loop,loop')
                push!(gaugeaction,factor,loop)
            elseif typeof(factor) <: Number
                push!(gaugeaction,factor,loop)
                push!(gaugeaction,factor',loop')
            else
                error("type of factor $(typeof(factor)) is not supported")
            end
            #append!(loop,loop')
            #factor = linkvalues[i]
            #push!(gaugeaction,factor,loop)
        end

        return new{Dim,typeof(F0),T}(Nflow,eps,gaugeaction,Ftemps,tempG,Utemps)
    end
end

mutable struct Gradientflow{TA,T} <: Abstractsmearing
    Nflow::Int64
    eps::Float64
    _temporal_TA_field::Array{TA,1}
    _temporal_G_field::Array{T,1}
    _temporal_U_field::Array{Array{T,1},1}

    function Gradientflow(U::Array{T,1};Nflow = 1,eps = 0.01,mpi = false) where T <: AbstractGaugefields
        F0 = initialize_TA_Gaugefields(U)
        Ftemps = Array{typeof(F0),1}(undef,4)
        Ftemps[1] = F0
        for i=2:4
            Ftemps[i] = initialize_TA_Gaugefields(U)
        end

        Utemps = Array{Array{T,1},1}(undef,2)
        for i=1:2
            Utemps[i] = similar(U)
        end

        tempG = Array{T,1}(undef,3)
        for i=1:3
            tempG[i] = similar(U[1])
        end

        return new{typeof(F0),T}(Nflow,eps,Ftemps,tempG,Utemps)
    end

end

function get_tempG(x::T) where T <: Gradientflow
    return x._temporal_G_field
end

function get_eps(x::T) where T <: Gradientflow
    return x.eps
end

function flow!(U,g::T) where T<: Gradientflow
    #@assert Dim == 4 "Dimension should be 4. But Dim = $Dim "
    Ftemps = g._temporal_TA_field
    Utemps = g._temporal_U_field
    temps = g._temporal_G_field

    F0 = Ftemps[1]
    F1 = Ftemps[2]
    F2 = Ftemps[3]
    Ftmp = Ftemps[4]

    #Ftmp = similar(U)
    W1 = Utemps[1]
    W2 = Utemps[2]
    temp1 = temps[1]
    temp2 = temps[2]
    temp3 = temps[3]
    eps = g.eps

    for istep=1:g.Nflow #RK4 integrator -> RK3?
        clear_U!(F0)
        add_force!(F0,U,temps,plaqonly = true)

        #add_force!(F0,U,[temp1,temp2,temp3],gparam)

        exp_aF_U!(W1,-eps*(1/4),F0,U,[temp1,temp2,temp3]) #exp(a*F)*U

        #println("W1 ",W1[1][1,1,1,1,1,1])
        #
        clear_U!(F1)
        add_force!(F1,W1,[temp1,temp2,temp3],plaqonly = true)
        #add_force!(F1,W1,[temp1,temp2,temp3],gparam) #F
        #println("F1 ",F1[1][1,1,1,1,1,1])
        clear_U!(Ftmp)
        add_U!(Ftmp,-(8/9*eps),F1)
        #println("Ftmp ",Ftmp[1][1,1,1,1,1,1])
        add_U!(Ftmp,(17/36*eps),F0)
        #println("Ftmp1 ",Ftmp[1][1,1,1,1,1,1])
        exp_aF_U!(W2,1,Ftmp,W1,[temp1,temp2,temp3]) #exp(a*F)*U
        #exp_aF_U!(W2,1,Ftmp,U,[temp1,temp2,temp3]) #exp(a*F)*U
        #println("W2 ",W2[1][1,1,1,1,1,1])

        #
        clear_U!(F2)
        add_force!(F2,W2,[temp1,temp2,temp3],plaqonly = true)
        #add_force!(F2,W2,[temp1,temp2,temp3],gparam) #F
        #calc_gaugeforce!(F2,W2,univ) #F
        clear_U!(Ftmp)
        add_U!(Ftmp,-(3/4*eps),F2)
        add_U!(Ftmp,(8/9*eps),F1)
        add_U!(Ftmp,-(17/36*eps),F0)
        #exp_aF_U!(W1,1,Ftmp,U,[temp1,temp2,temp3]) #exp(a*F)*U
        exp_aF_U!(U,1,Ftmp,W2,[temp1,temp2,temp3]) #exp(a*F)*U

        #println(typeof(U[1]))
        #println(U[1][1,1,1,1,1,1])

        #error("U")
    end

end

function flow!(U,g::Gradientflow_general{Dim,TA,T}) where {Dim,TA,T}
    #@assert Dim == 4 "Dimension should be 4. But Dim = $Dim "
    Ftemps = g._temporal_TA_field
    Utemps = g._temporal_U_field
    temps = g._temporal_G_field

    F0 = Ftemps[1]
    F1 = Ftemps[2]
    F2 = Ftemps[3]
    Ftmp = Ftemps[4]

    #Ftmp = similar(U)
    W1 = Utemps[1]
    W2 = Utemps[2]
    temp1 = temps[1]
    temp2 = temps[2]
    temp3 = temps[3]
    eps = g.eps

    for istep=1:g.Nflow #RK4 integrator -> RK3?
        clear_U!(F0)

        F_update!(F0,U,1,Dim,g.gaugeaction)

        #add_force!(F0,U,temps,plaqonly = true)

        #add_force!(F0,U,[temp1,temp2,temp3],gparam)

        exp_aF_U!(W1,-eps*(1/4),F0,U,[temp1,temp2,temp3]) #exp(a*F)*U

        #println("W1 ",W1[1][1,1,1,1,1,1])
        #
        clear_U!(F1)
        F_update!(F1,W1,1,Dim,g.gaugeaction)
        #add_force!(F1,W1,[temp1,temp2,temp3],plaqonly = true)
        #add_force!(F1,W1,[temp1,temp2,temp3],gparam) #F
        #println("F1 ",F1[1][1,1,1,1,1,1])
        clear_U!(Ftmp)
        add_U!(Ftmp,-(8/9*eps),F1)
        #println("Ftmp ",Ftmp[1][1,1,1,1,1,1])
        add_U!(Ftmp,(17/36*eps),F0)
        #println("Ftmp1 ",Ftmp[1][1,1,1,1,1,1])
        exp_aF_U!(W2,1,Ftmp,W1,[temp1,temp2,temp3]) #exp(a*F)*U
        #exp_aF_U!(W2,1,Ftmp,U,[temp1,temp2,temp3]) #exp(a*F)*U
        #println("W2 ",W2[1][1,1,1,1,1,1])

        #
        clear_U!(F2)
        F_update!(F2,W2,1,Dim,g.gaugeaction)
        #add_force!(F2,W2,[temp1,temp2,temp3],plaqonly = true)
        #add_force!(F2,W2,[temp1,temp2,temp3],gparam) #F
        #calc_gaugeforce!(F2,W2,univ) #F
        clear_U!(Ftmp)
        add_U!(Ftmp,-(3/4*eps),F2)
        add_U!(Ftmp,(8/9*eps),F1)
        add_U!(Ftmp,-(17/36*eps),F0)
        #exp_aF_U!(W1,1,Ftmp,U,[temp1,temp2,temp3]) #exp(a*F)*U
        exp_aF_U!(U,1,Ftmp,W2,[temp1,temp2,temp3]) #exp(a*F)*U

        #println(typeof(U[1]))
        #println(U[1][1,1,1,1,1,1])

        #error("U")
    end

end

function F_update!(F,U,factor,Dim,gauge_action) # F -> F +factor*U*dSdUμ
    NC = U[1].NC
    temps = get_temporary_gaugefields(gauge_action)
    dSdUμ = temps[end]

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

end

end


text:
using Molly
using Test

timestep = 0.0002
temperature = 298
n_steps = 50

forcefield, molecule, coords, box_size = readinputs(
            normpath(@__DIR__, "..", "data", "5XER", "gmx_top_ff.top"),
            normpath(@__DIR__, "..", "data", "5XER", "gmx_coords.gro"))

s = Simulation(forcefield, molecule, coords, box_size, temperature,
            timestep, n_steps)

simulate!(s)

writepdb("test.pdb", s.universe)
rm("test.pdb")


text:
module Survival

using Distributions
using LinearAlgebra
using PositiveFactorizations
using StatsBase
using StatsModels

export
    EventTime,
    isevent,
    iscensored,

    NonparametricEstimator,
    KaplanMeier,
    NelsonAalen,
    fit,
    confint,

    CoxModel,
    coxph,
    coef,
    loglikelihood,
    nullloglikelihood,
    nobs,
    dof,
    vcov,
    stderror

abstract type AbstractEstimator end
abstract type NonparametricEstimator <: AbstractEstimator end

include("eventtimes.jl")
include("estimator.jl")
include("kaplanmeier.jl")
include("nelsonaalen.jl")
include("cox.jl")
include("optimization.jl")

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Small sanity tests to ensure changing the rounding of float functions work
using Base.MathConstants

using Test

@testset "Float64 checks" begin
    # a + b returns a number exactly between prevfloat(1.) and 1., so its
    # final result depends strongly on the utilized rounding direction.
    a = prevfloat(0.5)
    b = 0.5
    c = 0x1p-54
    d = prevfloat(1.)

    @testset "Default rounding direction, RoundNearest" begin
        @test a + b === 1.
        @test - a - b === -1.
        @test a - b === -c
        @test b - a === c
    end
    @testset "RoundToZero" begin
        setrounding(Float64,RoundToZero) do
            @test a + b === d
            @test - a - b === -d
            @test a - b === -c
            @test b - a === c
        end
        # Sanity check to see if we have returned to RoundNearest
        @test a + b === 1.
        @test - a - b === -1.
        @test a - b == -c
        @test b - a == c
    end

    @testset "RoundUp" begin
        setrounding(Float64,RoundUp) do
            @test a + b === 1.
            @test - a - b === -d
            @test a - b === -c
            @test b - a === c
        end
    end
    @testset "RoundDown" begin
        setrounding(Float64,RoundDown) do
            @test a + b === d
            @test - a - b === -1.
            @test a - b === -c
            @test b - a === c
        end
    end
end

@testset "Float32 checks" begin
    a32 = prevfloat(0.5f0)
    b32 = 0.5f0
    c32 = (1.f0 - prevfloat(1.f0))/2
    d32 = prevfloat(1.0f0)

    @testset "Default rounding direction, RoundNearest" begin
        @test a32 + b32 === 1.0f0
        @test - a32 - b32 === -1.0f0
        @test a32 - b32 === -c32
        @test b32 - a32 === c32
    end
    @testset "RoundToZero" begin
        setrounding(Float32,RoundToZero) do
            @test a32 + b32 === d32
            @test - a32 - b32 === -d32
            @test a32 - b32 === -c32
            @test b32 - a32 === c32
        end

        # Sanity check to see if we have returned to RoundNearest
        @test a32 + b32 === 1.0f0
        @test - a32 - b32 === -1.0f0
        @test a32 - b32 == -c32
        @test b32 - a32 == c32
    end
    @testset "RoundUp" begin
        setrounding(Float32,RoundUp) do
            @test a32 + b32 === 1.0f0
            @test - a32 - b32 === -d32
            @test a32 - b32 === -c32
            @test b32 - a32 === c32
        end
    end
    @testset "RoundDown" begin
        setrounding(Float32,RoundDown) do
            @test a32 + b32 === d32
            @test - a32 - b32 === -1.0f0
            @test a32 - b32 === -c32
            @test b32 - a32 === c32
        end
    end
end

@testset "convert with rounding" begin
    for v = [sqrt(2),-1/3,nextfloat(1.0),prevfloat(1.0),nextfloat(-1.0),
             prevfloat(-1.0),nextfloat(0.0),prevfloat(0.0)]
        pn = Float32(v,RoundNearest)
        @test pn == convert(Float32,v)
        pz = Float32(v,RoundToZero)
        @test pz == setrounding(()->convert(Float32,v), Float64, RoundToZero)
        pd = Float32(v,RoundDown)
        @test pd == setrounding(()->convert(Float32,v), Float64, RoundDown)
        pu = Float32(v,RoundUp)
        @test pu == setrounding(()->convert(Float32,v), Float64, RoundUp)

        @test pn == pd || pn == pu
        @test v > 0 ? pz == pd : pz == pu
        @test pu - pd == eps(pz)
    end

    for T in [Float32,Float64]
        for v in [sqrt(big(2.0)),-big(1.0)/big(3.0),nextfloat(big(1.0)),
                  prevfloat(big(1.0)),nextfloat(big(0.0)),prevfloat(big(0.0)),
                  pi,ℯ,eulergamma,catalan,golden,
                  typemax(Int64),typemax(UInt64),typemax(Int128),typemax(UInt128),0xa2f30f6001bb2ec6]
            pn = T(v,RoundNearest)
            @test pn == convert(T,BigFloat(v))
            pz = T(v,RoundToZero)
            @test pz == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundToZero)
            pd = T(v,RoundDown)
            @test pd == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundDown)
            pu = T(v,RoundUp)
            @test pu == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundUp)

            @test pn == pd || pn == pu
            @test v > 0 ? pz == pd : pz == pu
            @test isinf(pu) || pu - pd == eps(pz)
        end
    end
end
@testset "fenv" begin
    @test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundNearest)) == RoundNearest
    @test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundToZero)) == RoundToZero
    @test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundUp)) == RoundUp
    @test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundDown)) == RoundDown
    @test_throws ArgumentError Base.Rounding.from_fenv(-99)
end

@testset "round error throwing" begin
    badness = 1//0
    @test_throws DivideError round(Int64,badness,RoundNearestTiesAway)
    @test_throws DivideError round(Int64,badness,RoundNearestTiesUp)
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity

using OceanTurbulenceParameterEstimation

#####
##### Parameters
#####

Nz = 16
Lz = 128
Qᵇ = 1e-8
Qᵘ = -1e-5
Δt = 20.0
f₀ = 1e-4
N² = 1e-5

stop_time = 6Δt
save_interval = 2Δt
N_ensemble = 3

experiment_name = "convective_adjustment_test"
data_path = experiment_name * ".jld2"

# "True" parameters to be estimated by calibration
convective_κz = 1.0
convective_νz = 0.9
background_κz = 1e-4
background_νz = 1e-5

@testset "OneDimensionalTimeSeries" begin
    #####
    ##### Generate synthetic observations
    #####

    grid = RegularRectilinearGrid(size=Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded))
    closure = ConvectiveAdjustmentVerticalDiffusivity(; convective_κz, background_κz, convective_νz, background_νz)
    coriolis = FPlane(f=f₀)

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        tracers = :b,
                                        buoyancy = BuoyancyTracer(),
                                        boundary_conditions = (; u=u_bcs, b=b_bcs),
                                        coriolis = coriolis,
                                        closure = closure)

    set!(model, b = (x, y, z) -> N² * z)

    simulation = Simulation(model; Δt, stop_time)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          schedule = TimeInterval(save_interval),
                                                          prefix = experiment_name,
                                                          field_slicer = nothing,
                                                          force = true)

    run!(simulation)

    #####
    ##### Load truth data as observations
    #####

    data_path = experiment_name * ".jld2"

    b_observations = OneDimensionalTimeSeries(data_path, field_names=:b)
    ub_observations = OneDimensionalTimeSeries(data_path, field_names=(:u, :b))
    uvb_observations = OneDimensionalTimeSeries(data_path, field_names=(:u, :v, :b))

    @test keys(b_observations.field_time_serieses) == tuple(:b)
    @test keys(ub_observations.field_time_serieses) == tuple(:u, :b)
    @test keys(uvb_observations.field_time_serieses) == tuple(:u, :v, :b)
end


text:
function compute()

	n=BigInt(100)
	str=string(factorial(n));
	ans=0
	for i=1:length(str)
		ans+=str[i]-'0'
	end
	return ans

end

println(compute())


text:
module mirone_

export
	mirone

using IUP
using IUP_IM
using IUP_CD
using GMT

const NEWPOINT = 0
const MOVE = 1
const CLOSE = 2
const CENTER = 3
const REPAINT = 4

type LineSegPix
	x1::Cint
	x2::Cint
	y1::Cint
	y2::Cint
	active::Bool
end
#function LineSegPix(x1, x2, y1, y2, active)
#	x1 = 0; x2 = 0; y1 = 0; y2 = 0; active = false
#end

include("draw_poly.jl")
include("mirone_uis.jl")
include("scaleto8.jl")

global poly_lastwhat = CLOSE

# ----------------------------------------------------------------------------------
function mirone()
	IupOpen()	#Initializes IUP

	imCreateButtonImages()

 	mainMenu = CreateMainMenu()

	# Associates handle "mainMenu" with mainMenu
	IupSetHandle("mainMenu", mainMenu)

	btNew = IupButton("", "imImgWinNew" )
	IupSetAttribute(btNew,"TIP","Creates a new image.");
	IupSetAttribute(btNew,"IMAGE","imImgWinNewButton");

	btOpen = IupButton("", "imImgWinOpen" )
	IupSetAttribute(btOpen,"TIP","Loads an image file from disk.")
	IupSetAttribute(btOpen,"IMAGE","imImgWinOpenButton")

	btLine = IupButton("", "imImgWinLine" )
	IupSetAttribute(btLine,"TIP","Draw line.")
	IupSetAttribute(btLine,"IMAGE","imImgWinLineButton")

	btCloseAll = IupButton("", "imImgWinCloseAll" );
	IupSetAttribute(btCloseAll,"TIP","Close all windows.");
	IupSetAttribute(btCloseAll,"IMAGE","imImgWinCloseAllButton");
	IupSetAttribute(btCloseAll, "ACTIVE", "NO");

	hbToolBar = IupHbox(
		IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO"),	# Blank space before first icon
		btNew,
		btOpen,
		btLine,
		IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO"),
		#IupSetAttributes(IupFill(),"EXPAND=YES"),
		IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO")
	)
	IupSetAttribute(hbToolBar, "btCloseAll", btCloseAll);
	IupSetAttribute(hbToolBar, "FLAT", "YES");

	cnvStatus = IupCanvas();
	IupSetAttribute(cnvStatus, "EXPAND", "HORIZONTAL");
	IupSetAttribute(cnvStatus, "RASTERSIZE", "1x25");
	IupSetAttribute(cnvStatus, "BORDER", "YES");

	iup_canvas = blabla()

	vbMainDesktop = IupVbox(
		IupSetAttributes(IupFill(),"SIZE=2, EXPAND=NO"),
		hbToolBar,
		IupSetAttributes(IupFill(),"SIZE=4, EXPAND=NO"),
		IupHbox(iup_canvas),
		IupSetAttributes(IupFill(),"SIZE=4, EXPAND=NO"),
		cnvStatus,
		IupSetAttributes(IupFill(),"SIZE=2, EXPAND=NO")
	)

	# Creates a dialog
	main_window = IupDialog(vbMainDesktop)

	handles = Handles(C_NULL, C_NULL, C_NULL)
	handles.figure1 = main_window
	handles.iup_canvas = iup_canvas
	guidata(main_window, handles)

	IupSetAttribute(main_window, "MENU","mainMenu")
	#IupSetAttribute(main_window, "CLOSE_CB","imImgWinMainExit");
	#IupSetAttribute(main_window, "DROPFILES_CB","imImgWinDropFile");
	#IupSetAttribute(main_window, "GETFOCUS_CB", "imImgWinGetFocus");

	IupSetAttribute(main_window,"MainToolBar", hbToolBar);

	#IupSetAttribute(main_window,"PLACEMENT", "MAXIMIZED");
	IupSetAttribute(main_window,"BRINGFRONT", "YES");
	IupSetAttribute(main_window, "K_cN","imImgWinWindowNext");
	IupSetAttribute(main_window, "K_cO","imImgWinOpen");

	# Sets title and size of the dialog and associates a menu to it
	IupSetAttributes(main_window, "TITLE=\"Mirone toy\", RASTERSIZE = 730x91, MENU = mainMenu");

	# associa os callbacks
	setcallbacks(handles)

	IupShowXY(main_window, IUP_CENTER, IUP_TOP) # Displays the dialog (at screen center)
	IupMainLoop()                               # Initializes IUP main loop
	IupClose()                                  # And close it when ready

end

# ---------------------------------------------------------------------------------------------
function blabla()
	iup_canvas = IupCanvas()
	IupSetAttribute(iup_canvas, "BORDER", "NO");
	IupSetCallback(iup_canvas, "ACTION", cfunction(cbCanvasRepaint, Int, (Ptr{Ihandle},)))
	IupSetCallback(iup_canvas, "MAP_CB", cfunction(cbCanvasMap, Int, (Ptr{Ihandle},)))
	return iup_canvas
end

# --------------------------------------------------------------------------------
function cbCanvasMap(iup_canvas::Ptr{Ihandle})
	handles = guidata(iup_canvas)
	cd_canvas = cdCreateCanvas(cdContextIup(), iup_canvas)
	#cd_canvas = cdCreateCanvas(cdContextDBuffer(), iup_canvas)
	IupSetAttribute(IupGetDialog(iup_canvas), "cdCanvas", cd_canvas)
	IupSetAttribute(iup_canvas, "cdCanvas", cd_canvas)			# Redundant with handles
	handles.cd_canvas = cd_canvas
	guidata(handles.figure1, handles)
	return IUP_DEFAULT
end

# ---------------------------------------------------------------------------------------------
function setcallbacks(handles::Handles)
	IupSetFunction("imImgWinOpen",  cfunction(cmOpen, Int, (Ptr{Ihandle},)))
	IupSetFunction("imImgWinLine",  cfunction(cmdDrawPoly, Int, (Ptr{Ihandle},)))

	#IupSetFunction("cmdRepaint", (Icallback) fRepaint);
	WindowButtonDownFcn(handles.iup_canvas, fButtonCB)
	WindowButtonMotionFcn(handles.iup_canvas, fMotionCB)
	#IupSetFunction("cmdResizeCB", (Icallback) fResizeCB);
end

# ---------------------------------------------------------------------------------------------
function cmOpen(iup_canvas::Ptr{Ihandle})
	# Retrieve a file name
	FileName = "*.*" * repeat(" ",253)	# File name will be returned in this variable, hence its size.
	if (IupGetFile(FileName) != 0)
		return IUP_DEFAULT
	end

	ShowImage(FileName, IupGetDialog(iup_canvas))
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function ShowGrid(FileName::String, iup_dialog::Ptr{Ihandle})
	# ...
	imgGrd = IupGetAttribute(iup_dialog, "imImage")		# typeof(imImage) => Ptr{Uint8}
	imgGrd = convert(Ptr{imImage}, imgGrd)		# If I use Ptr{imImage} it Booms???
	if (imgGrd != C_NULL)
		imImageDestroy(imgGrd)
	end
	IupSetAttribute(iup_dialog, "imImage")

	G   = load_grd_with_gmt(FileName)
	img = scaleto8(G)

	handles = guidata(iup_dialog)
	IupSetAttribute(handles.figure1, "CLIENTSIZE", @sprintf("%dx%d", max(G.n_columns,730), G.n_rows+72))
	IupRefresh(handles.figure1)

	IupSetAttribute(iup_dialog, "imImage", imgGrd);
	cbCanvasRepaint(iup_dialog)   # we can do this because canvas inherit attributes from the dialog */
	IupShow(iup_dialog)
end

# --------------------------------------------------------------------------------
function ShowImage(FileName::String, iup_dialog::Ptr{Ihandle})
	image = IupGetAttribute(iup_dialog, "imImage")		# typeof(image) => Ptr{Uint8}
	image = convert(Ptr{imImage}, image)		# If I use Ptr{imImage} it Booms???
	if (image != C_NULL)
		imImageDestroy(image)
	end
	IupSetAttribute(iup_dialog, "imImage")

#=
	error = pointer([0])
	image = imFileImageLoadBitmap(FileName, 0, error)
	error = unsafe_load(error)
	if (error != 0)	PrintError(error)	end
	if (image == C_NULL) return end
=#
	image = load_grd_with_gmt(FileName)
	img   = unsafe_load(image)		# Need to get access to the composite type, not to its pointer

	handles = guidata(iup_dialog)
	IupSetAttribute(handles.figure1, "CLIENTSIZE", @sprintf("%dx%d", max(img.width,730), img.height+72))
	IupRefresh(handles.figure1)
#@show(bytestring(IupGetAttribute(handles.figure1, "RASTERSIZE")))

	IupSetAttribute(iup_dialog, "imImage", image);
	cbCanvasRepaint(iup_dialog)   # we can do this because canvas inherit attributes from the dialog */
	IupShow(iup_dialog)
end

# --------------------------------------------------------------------------------
function load_img_with_gmt(FileName::String)
	# Load an image with GDAL but via GMT
	I = gmt("read -Ti " * FileName)
	n = I.n_columns * I.n_rows
	if (I.LayerCount >= 3)
		data = pointer([pointer(I.image[:,:,1]), pointer(I.image[:,:,2]), pointer(I.image[:,:,3])])
		pPalette = C_NULL
	else
		data = pointer([pointer(I.image)])
	end
	color_space = IM_RGB
	if (I.nColors > 1)
		color_space = IM_MAP
		cm = I.colormap
		pal = zeros(UInt32, I.nColors)
		for (k = 1:I.nColors)
			pal[k] =  (UInt32(cm[k*4]) << 24) | (UInt32(cm[k*4-3]) << 16) | (UInt32(cm[k*4-2]) << 8) | UInt32(cm[k*4-1])
		end
		pPalette = pointer(pal)
	end

	im = imImage(I.n_columns, I.n_rows, color_space, 0, 0, I.LayerCount, I.n_columns, n,
	             n*I.LayerCount, n, data, pPalette, I.nColors, C_NULL)
	image = pointer([im])
	return gc_ref(image)
end

# --------------------------------------------------------------------------------
function load_grd_with_gmt(FileName::String)
	# Load a grid with GMT (it may use GDAL)
	G = gmt("read -Tg " * FileName)
	pImg = scaleto8(G)
	return pImg
end

# --------------------------------------------------------------------------------
function cbCanvasRepaint(iup_canvas::Ptr{Ihandle})
	cd_canvas = convert(Ptr{cdCanvas}, IupGetAttribute(iup_canvas, "cdCanvas"))
	if (cd_canvas == C_NULL)
		return IUP_DEFAULT
	end

	cdCanvasActivate(cd_canvas);
	cdCanvasClear(cd_canvas);

	image = convert(Ptr{imImage}, IupGetAttribute(iup_canvas, "imImage"))
	if (image == C_NULL)
		return IUP_DEFAULT
	end

	#handles = guidata(iup_canvas)
	img = unsafe_load(image)		# Need to get access to the composite type, not to its pointer
	x_off = div(730 - img.width, 2)
	imcdCanvasPutImage(cd_canvas, img, x_off, 0, img.width, img.height, 0, 0, 0, 0)

	IupRefresh(iup_canvas)
	cdCanvasFlush(cd_canvas)
	return IUP_DEFAULT
end

end 	# module


text:
using FactCheck, Base.Test
import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    GU.read("pyrasterio/RGB.byte.tif") do raster
        @fact GU.getproj(raster) --> "PROJCS[\"UTM Zone 18, Northern Hemisphere\",GEOGCS[\"Unknown datum based upon the WGS 84 ellipsoid\",DATUM[\"Not_specified_based_on_WGS_84_spheroid\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-75],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]"
        @fact GU.getgeotransform(raster) --> roughly([101985.0,300.0379266750948,0.0,2.826915e6,0.0,-300.041782729805])
        @fact size(GU.fetch(raster, 1)) --> (791,718)
        @fact size(GU.fetch(raster, Cint[1,3])) --> (791,718,2)
        @fact size(GU.fetch(raster)) --> (791,718,3)

        band = GU.fetchband(raster, 1)
        @fact GU.getdatatype(band) --> UInt8
        @fact GU.width(band) --> 791
        @fact GU.height(band) --> 718
        band_color = GU.getcolorinterp(band)
        @fact GU.getname(band_color) --> "Red"
    end

    GU.read("pyrasterio/RGB.byte.tif") do raster
        w = GU.fetch(raster, 1, 1:100, 1:100)
        @fact size(w) --> (791, 718)
    end

    GU.create("pyrasterio/example.tif",
              "GTiff", # drivername
              width=500,
              height=300,
              nbands=1,
              dtype=UInt8 # DataType
            ) do raster
        image = fill(UInt8(127), (150, 250))
        GU.update!( raster,
                    image, # image to "burn" into the raster
                    1, # update band 1
                    30:180, # along (window) xcoords 30 to 180
                    50:300) # along (window) ycoords 30 to 180
    end

    GU.read("pyrasterio/RGB.byte.tif") do src
    GU.create("pyrasterio/example2.tif", "GTiff",
              width=500, height=300, nbands=3, dtype=UInt8) do dst
        rgb = GU.fetch(src, Cint[1,2,3]) # fetch bands 1 - 3
        # You can update all 3 bands simultaneously
        GU.update!(dst, rgb, # image to "burn" into destination dataset
                   Cint[1,2,3], # indices of the bands to be updated
                   30:269, # along (window) xcoords 30 to 269
                   50:313) # along (window) ycoords 50 to 313
    end
    end

    GU.read("pyrasterio/RGB.byte.tif") do src
        rgb = GU.fetch(src, Cint[1,2,3], 350:410, 350:450)
        GU.create("pyrasterio/example3.tif", "GTiff",
                  width=500, height=300, nbands=3, dtype=eltype(rgb)) do dst
            GU.update!(dst, rgb, Cint[1,2,3], 1:240, 1:400)
        end
    end
end

for filename in ["data1a.jl",
                 "data1b.jl",
                 "data2a.jl",
                 "data2b.jl",
                 "data3a.jl",
                 "data3b.jl",
                 "data4a.jl",
                 "data4b.jl",
                 "data5a.jl",
                 "data5b.jl",
                 "tutorial_raster.jl",
                 "tutorial_vector.jl"]
    include(filename)
end


text:
using Documenter
using CaimCore

DocMeta.setdocmeta!(CaimCore, :DocTestSetup, :(using CaimCore); recursive=true)
makedocs(
    sitename = "CaimCore",
    format = Documenter.HTML(),
    modules = [CaimCore],
    authors = "Douglas G. Moore",
    pages = Any[
        "Home" => "index.md",
        "API Reference" => "api.md",
    ]
)

deploydocs(
    repo = "github.com/dglmoore/CaimCore.jl.git"
)


text:
"""
    linearized_air_pressure(ρ, ρe_tot, ρe_pot, ρq_tot=0, ρq_liq=0, ρq_ice=0)

The air pressure, linearized around a dry rest state, from the equation of state
(ideal gas law) where:

 - `ρ` (moist-)air density
 - `ρe_tot` total energy density
 - `ρe_pot` potential energy density
and, optionally,
 - `ρq_tot` total water density
 - `ρq_liq` liquid water density
 - `ρq_ice` ice density
"""
function linearized_air_pressure(
    ρ::FT,
    ρe_tot::FT,
    ρe_pot::FT,
    ρq_tot::FT = FT(0),
    ρq_liq::FT = FT(0),
    ρq_ice::FT = FT(0),
) where {FT <: Real}
    ρ * FT(R_d) * FT(T_0) +
    FT(R_d) / FT(cv_d) * (
        ρe_tot - ρe_pot - (ρq_tot - ρq_liq) * FT(e_int_v0) +
        ρq_ice * (FT(e_int_i0) + FT(e_int_v0))
    )
end

@inline function linearized_pressure(
    ::DryModel,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    linearized_air_pressure(state.ρ, state.ρe, ρe_pot)
end
@inline function linearized_pressure(
    ::EquilMoist,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    linearized_air_pressure(state.ρ, state.ρe, ρe_pot, state.moisture.ρq_tot)
end

abstract type AtmosLinearModel <: BalanceLaw end

vars_state(lm::AtmosLinearModel, FT) = vars_state(lm.atmos, FT)
vars_gradient(lm::AtmosLinearModel, FT) = @vars()
vars_diffusive(lm::AtmosLinearModel, FT) = @vars()
vars_aux(lm::AtmosLinearModel, FT) = vars_aux(lm.atmos, FT)
vars_integrals(lm::AtmosLinearModel, FT) = @vars()
vars_reverse_integrals(lm::AtmosLinearModel, FT) = @vars()

function update_aux!(
    dg::DGModel,
    lm::AtmosLinearModel,
    Q::MPIStateArray,
    t::Real,
)
    return false
end
function flux_diffusive!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
integral_load_aux!(lm::AtmosLinearModel, integ::Vars, state::Vars, aux::Vars) =
    nothing
integral_set_aux!(lm::AtmosLinearModel, aux::Vars, integ::Vars) = nothing
reverse_integral_load_aux!(
    lm::AtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
reverse_integral_set_aux!(lm::AtmosLinearModel, aux::Vars, integ::Vars) =
    nothing
flux_diffusive!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
) = nothing
function wavespeed(lm::AtmosLinearModel, nM, state::Vars, aux::Vars, t::Real)
    ref = aux.ref_state
    return soundspeed_air(ref.T, lm.atmos.param_set)
end

function boundary_state!(
    nf::NumericalFluxNonDiffusive,
    atmoslm::AtmosLinearModel,
    args...,
)
    atmos_boundary_state!(nf, AtmosBC(), atmoslm, args...)
end
function boundary_state!(
    nf::NumericalFluxDiffusive,
    atmoslm::AtmosLinearModel,
    args...,
)
    nothing
end
init_aux!(lm::AtmosLinearModel, aux::Vars, geom::LocalGeometry) = nothing
init_state!(lm::AtmosLinearModel, state::Vars, aux::Vars, coords, t) = nothing

struct AtmosAcousticLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticLinearModel(atmos::M) where {M}
        if atmos.ref_state === NoReferenceState()
            error("AtmosAcousticLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end

function flux_nondiffusive!(
    lm::AtmosAcousticLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    ref = aux.ref_state
    e_pot = gravitational_potential(lm.atmos.orientation, aux)

    flux.ρ = state.ρu
    pL =
        linearized_pressure(lm.atmos.moisture, lm.atmos.orientation, state, aux)
    flux.ρu += pL * I
    flux.ρe = ((ref.ρe + ref.p) / ref.ρ - e_pot) * state.ρu
    nothing
end
function source!(
    lm::AtmosAcousticLinearModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

struct AtmosAcousticGravityLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticGravityLinearModel(atmos::M) where {M}
        if atmos.ref_state === NoReferenceState()
            error("AtmosAcousticGravityLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end
function flux_nondiffusive!(
    lm::AtmosAcousticGravityLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    ref = aux.ref_state
    e_pot = gravitational_potential(lm.atmos.orientation, aux)

    flux.ρ = state.ρu
    pL =
        linearized_pressure(lm.atmos.moisture, lm.atmos.orientation, state, aux)
    flux.ρu += pL * I
    flux.ρe = ((ref.ρe + ref.p) / ref.ρ) * state.ρu
    nothing
end
function source!(
    lm::AtmosAcousticGravityLinearModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    ∇Φ = ∇gravitational_potential(lm.atmos.orientation, aux)
    source.ρu -= state.ρ * ∇Φ
    nothing
end

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
using Test, Random, MPStates, LinearAlgebra

# These are just integration tests to test that the DMRG algorithms don't error.
@testset "Integration tests: variational/DMRG algorithms" begin
@testset "Variational simplification of Mps" begin
    rtest1 = MPStates.testMps("rtest1")
    simplify!(rtest1, 1)
end

@testset "DMRG algorithms" for T in [Float64, ComplexF64]
    Random.seed!(1)

    # Diagonalize a simple Hamiltonian with first neighbor interactions and test
    # that energy and variance decrease after each sweep and that the method
    # doesn't error.
    L = 15
    J = Hermitian(convert.(T, randn(L, L)))
    d = 2
    H = Mpo(T, L, d)
    add_ops!(H, "b+", "b", convert.(T, J))

    # Max allowed bond dimension per sweep.
    maxm = [10, 20, 20]

    # DMRG1.
    psi = randomMps(T, L, d, 2)
    dmrg_opts = DMRGOpts("DMRG1", maxm, 1e-10, show_trace=0)
    E1, var1 = dmrg!(psi, H, dmrg_opts)

    # DMRG2.
    psi = randomMps(T, L, d, 2)
    dmrg_opts = DMRGOpts("DMRG2", maxm, 1e-10, show_trace=0)
    E2, var2 = dmrg!(psi, H, dmrg_opts)

    # DMRG3S.
    psi = randomMps(T, L, d, 2)
    dmrg_opts = DMRGOpts("DMRG3S", maxm, 1e-10, show_trace=0)
    E3, var3 = dmrg!(psi, H, dmrg_opts)
end
end # @testset "Integration tests: variational/DMRG algorithms"


text:
using Test

@testset "TESTSETS: Test read TCX dir" begin
    @testset "CASE: Test a string that is NOT a path" begin
        err, _ = TCX.parse_tcx_dir("I am not a path")
        @test err == TCX.SERVER_ERROR
    end

    @testset "CASE: Test dir with NO TCX file" begin
	      err, _ =  TCX.parse_tcx_dir(tempdir())
        @test err == TCX.NOT_FOUND
    end

    @testset "CASE: Test dir with TCX file" begin
        err, _ =  TCX.parse_tcx_dir(@__DIR__)
        @test err == TCX.OK
    end

    @testset "CASE: Test relative dir with TCX file" begin
        err, _ =  TCX.parse_tcx_dir(".")
        @test err == TCX.OK
    end

    @testset "CASE: Test DataFrame after process dir" begin
        err, ta =  TCX.parse_tcx_dir(".")
        @test (err == TCX.OK) & (size(getDataFrame(ta), 1) == 10847)
    end
end


text:
using Flux
using Flux: flatten
using ExplainableAI
using ExplainableAI: fuse_batchnorm
using Random

pseudorand(dims...) = rand(MersenneTwister(123), Float32, dims...)
batchsize = 50

# # Test `fuse_batchnorm` on Dense layer
ins = 20
outs = 10
dense = Dense(ins, outs; init=pseudorand)
bn_dense = BatchNorm(outs, relu; initβ=pseudorand, initγ=pseudorand)
model = Chain(dense, bn_dense)

# collect statistics
x = pseudorand(ins, batchsize)
Flux.trainmode!(model)
model(x)
Flux.testmode!(model)

dense_fused = @inferred fuse_batchnorm(dense, bn_dense)
@test dense_fused(x) ≈ model(x)

# # Test `fuse_batchnorm` on Conv layer
insize = (10, 10, 3)
conv = Conv((3, 3), 3 => 4; init=pseudorand)
bn_conv = BatchNorm(4, relu; initβ=pseudorand, initγ=pseudorand)
model = Chain(conv, bn_conv)

# collect statistics
x = pseudorand(insize..., batchsize)
Flux.trainmode!(model)
model(x)
Flux.testmode!(model)

conv_fused = @inferred fuse_batchnorm(conv, bn_conv)
@test conv_fused(x) ≈ model(x)

# # Test `canonize` on models
# Sequential BatchNorm layers should be fused until they create a Dense or Conv layer
# with non-linear activation function.
model = Chain(
    Conv((3, 3), 3 => 6),
    BatchNorm(6),
    Conv((3, 3), 6 => 2, identity),
    BatchNorm(2),
    BatchNorm(2, softplus),
    BatchNorm(2),
    flatten,
    Dense(72, 10),
    BatchNorm(10),
    BatchNorm(10),
    BatchNorm(10, relu),
    BatchNorm(10),
    Dense(10, 10, gelu),
    BatchNorm(10),
    softmax,
)
Flux.trainmode!(model)
model(x)
Flux.testmode!(model)
model_canonized = canonize(model)

# 6 of the BatchNorm layers should be removed and the ouputs should match
@test length(model_canonized) == length(model) - 6
@test model(x) ≈ model_canonized(x)


text:
function optimize!(structure::DeterministicEquivalent, optimizer::MOI.AbstractOptimizer, x₀::AbstractVector)
    # Sanity check
    backend(structure.model) === optimizer || error("Stochastic program optimizer has not been connected to the deterministically equivalent problem.")
    # Crash if supported
    if MOI.supports(optimizer, MOI.VariablePrimalStart(), MOI.VariableIndex)
        for (i, dvar) in enumerate(all_decision_variables(structure.model, 1))
            MOI.set(optimizer, MOI.VariablePrimalStart(), index(dvar), x₀[i])
        end
    end
    # Run standard MOI optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function supports_structure(::MOI.AbstractOptimizer, ::DeterministicEquivalent)
    return true
end

function check_loadable(::MOI.AbstractOptimizer, ::DeterministicEquivalent)
    return nothing
end

function load_structure!(::MOI.AbstractOptimizer, ::DeterministicEquivalent, ::AbstractVector)
    return nothing
end

function set_master_optimizer!(structure::DeterministicEquivalent, optimizer)
    # Ensure decision bridges are added
    for bridge_type in structure.model.bridge_types
        JuMP._moi_add_bridge(structure.model.moi_backend, bridge_type)
    end
    return nothing
end

function set_master_optimizer_attribute!(structure::DeterministicEquivalent, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(backend(structure.model), attr, value)
    return nothing
end

function set_subproblem_optimizer!(structure::DeterministicEquivalent, optimizer)
    return nothing
end

function set_subproblem_optimizer_attribute!(::DeterministicEquivalent, ::MOI.AbstractOptimizerAttribute, value)
    return nothing
end

function optimize!(structure::DeterministicEquivalent, optimizer::AbstractStructuredOptimizer)
    # Sanity check
    supports_structure(optimizer, typeof(structure)) || throw(UnsupportedStructure{typeof(optimizer), typeof(structure)}())
    # Ensure that no decisions are fixed in the first stage
    untake_decisions!(structure.model, structure.decision_variables[1])
    # Load structure
    load_structure!(optimizer, structure)
    # Run structure-exploiting optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function cache_solution!(stochasticprogram::StochasticProgram{2}, structure::DeterministicEquivalent{2}, optimizer::MOI.AbstractOptimizer)
    cache = solutioncache(stochasticprogram)
    # Cache main solution
    variables = decision_variables_at_stage(stochasticprogram, 1)
    constraints = decision_constraints_at_stage(stochasticprogram, 1)
    cache[:solution] = SolutionCache(backend(structure.model), variables, constraints)
    # Cache first-stage solution
    cache[:node_solution_1] = SolutionCache(backend(structure.model), variables, constraints)
    try
        Q = MOIU.eval_variables(structure.sub_objectives[1][1][2]) do idx
            return MOI.get(backend(structure.model), MOI.VariablePrimal(), idx)
        end
        cache[:node_solution_1].modattr[MOI.ObjectiveValue()] = Q
        cache[:node_solution_1].modattr[MOI.DualObjectiveValue()] = Q
    catch
    end
    # Cache scenario-dependent solutions
    variables = decision_variables_at_stage(stochasticprogram, 2)
    constraints = decision_constraints_at_stage(stochasticprogram, 2)
    for i in 1:num_scenarios(stochasticprogram, 2)
        key = Symbol(:node_solution_2_, i)
        cache[key] = SolutionCache(backend(structure.model))
        cache_model_attributes!(cache[key], backend(structure.model))
        cache_variable_attributes!(cache[key], backend(structure.model), variables) do index
            structure.variable_map[(index, i)]
        end
        cache_constraint_attributes!(cache[key], backend(structure.model), constraints) do ci
            if ci isa CI{SingleDecision}
                return typeof(ci)(structure.variable_map[(MOI.VariableIndex(ci.value), i)].value)
            else
                return structure.constraint_map[(ci, i)]
            end
        end
        # Cache subobjective
        if cache[key].modattr[MOI.TerminationStatus()] == MOI.OPTIMAL
            Q = 0.0
            try
                Qᵢ = MOIU.eval_variables(structure.sub_objectives[2][i][2]) do idx
                    return MOI.get(backend(structure.model), MOI.VariablePrimal(), idx)
                end
                cache[key].modattr[MOI.ObjectiveValue()] = Qᵢ
                cache[key].modattr[MOI.DualObjectiveValue()] = Qᵢ
                Q += Qᵢ
            catch
            end
            # Correct first-stage objective
            cache[:node_solution_1].modattr[MOI.ObjectiveValue()] -= Q
        end
    end
    return nothing
end