text:
    # Functions and methods for cheme::Cascade.HollowIonScheme computations

    """
    `Cascade.computeSteps(scheme::Cascade.HollowIonScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})`
        ... computes in turn all the requested (decay) transition amplitudes as well as AutoIonization.Line's, etc. for all
            pre-specified decay steps of the cascade. When compared with standard computations of these atomic
            processes, however, the amount of output is largely reduced and often just printed into the summary file.
            A set of  data::Cascade.DrRateData  is returned.
    """
    function computeSteps(scheme::Cascade.HollowIonScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})
        linesA = AutoIonization.Line[];    linesR = PhotoEmission.Line[];    linesC = ElectronCapture.Line[];    cOrbitals = Dict{Subshell, Orbital}()
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        nt = 0;   st = 0;   previousMeanEn = 0.
        for  step  in  stepList
            st = st + 1
            nc = length(step.initialMultiplet.levels) * length(step.finalMultiplet.levels)
            sa = "\n  $st) Perform $(string(step.process)) amplitude computations for up to $nc decay lines (without selection rules): "
            println(sa);    if  printSummary   println(iostream, sa)   end

            if      step.process == Basics.ElecCapture()
                # Compute continuum orbitals independently for all transitions in the given block.
                newLines = DielectronicCapture.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                               step.settings, output=true, printout=false)
                append!(linesC, newLines);    nt = length(linesC)
            elseif  step.process == Basics.Auger()
                # Compute continuum orbitals independently for all transitions in the given block.
                newLines = AutoIonization.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                              step.settings, output=true, printout=false)
                append!(linesA, newLines);    nt = length(linesA)
            elseif  step.process == Basics.Radiative()
                newLines = PhotoEmission.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.grid,
                                                             step.settings, output=true, printout=false)
                append!(linesR, newLines);    nt = length(linesR)
            else   error("Unsupported atomic process for cascade computations.")
            end
            sa = "     Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, giving now rise " *
                 "to a total of $nt $(string(step.process)) decay lines."
            println(sa);    if  printSummary   println(iostream, sa)   end
        end
        #
        data = Cascade.DecayData(linesR, linesA)
    end

    """
    `Cascade.determineSteps(scheme::Cascade.HollowIonScheme, comp::Cascade.Computation, capturedList::Array{Cascade.Block,1})`
        ... determines all step::Cascade.Step's that need to be computed for this decay cascade. It cycles through all decay processes
            of the given scheme and selects all pairs of blocks due to the selected cascade approach. It checks that at least
            on pair of levels supports a decay within the step. A stepList::Array{Cascade.Step,1} is returned, and for which subsequently
            all required transition amplitudes and rates/cross sections are computed.
    """
    function determineSteps(scheme::Cascade.HollowIonScheme, comp::Cascade.Computation, capturedList::Array{Cascade.Block,1})
        stepList = Cascade.Step[]
        if  comp.approach  in  [Cascade.AverageSCA(), Cascade.SCA()]
            @warn("Transition within a cascade block are not yet considered.")
            for  (ia, blocka) in enumerate(capturedList)
                for  (ib, blockb) in enumerate(capturedList)
                    if   ia == ib   continue     end
                    for  process  in  comp.scheme.processes
                        if      process == Basics.Radiative()
                            if  blocka.NoElectrons == blockb.NoElectrons   &&
                                Basics.determineMeanEnergy(blocka.multiplet) - Basics.determineMeanEnergy(blockb.multiplet) > 0.
                                @show  Basics.determineMeanEnergy(blocka.multiplet) - Basics.determineMeanEnergy(blockb.multiplet)
                                settings = PhotoEmission.Settings(PhotoEmission.Settings(); multipoles=scheme.multipoles)
                                push!( stepList, Cascade.Step(process, settings, blocka.confs, blockb.confs, blocka.multiplet, blockb.multiplet) )
                            end
                        end
                    end
                end
            end
            #
        else  error("Unsupported cascade approach.")
        end
        return( stepList )
    end

    """
    `Cascade.generateBlocks(scheme::Cascade.HollowIonScheme, comp::Cascade.Computation, confs::Array{Configuration,1}; printout::Bool=true)`
        ... generate all block::Cascade.Block's, that need to be computed for this electron-capture and subsequent stabilization (DR) cascade,
            and compute also the corresponding multiplets. The different cascade approches realized different strategies how these block are
            selected and computed. A blockList::Array{Cascade.Block,1} is returned.
    """
    function generateBlocks(scheme::Cascade.HollowIonScheme, comp::Cascade.Computation, confs::Array{Configuration,1}; printout::Bool=true)
        blockList = Cascade.Block[]
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        if    comp.approach == AverageSCA()
            sa = "\n* Generate blocks for hollow-ion cascade computations: \n" *
                 "\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: " *
                 "\n    + orbitals are generated independently for each block for a Dirac-Fock-Slater potential; " *
                 "\n    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; " *
                 "\n    + only the Coulomb interaction is considered for the electron capture. " *
                 "\n    + only E1 excitations are considered for the stabilization. \n"
            if  printout       println(sa)              end
            if  printSummary   println(iostream, sa)    end
            #
            for  confa  in confs
                print("  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")
                if  printSummary   println(iostream, "\n*  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")   end
                    ##x @show confa, comp.asfSettings
                    basis     = Basics.performSCF([confa], comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                    multiplet = Basics.perform("computation: mutiplet from orbitals, no CI, CSF diagonal", [confa],  basis.orbitals,
                                               comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
            end
        elseif    comp.approach == SCA()
            sa = "\n* Generate blocks for hollow-ion cascade computations: \n" *
                 "\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: " *
                 "\n    + each single configuration forms an individual cascade block; " *
                 "\n    + orbitals are generated independently for each block for a Dirac-Fock-Slater potential; " *
                 "\n    + configuration mixing is included for each block, based on H^(DC); " *
                 "\n    + all requested multipoles are considered for the stabilization. \n"
            if  printout       println(sa)              end
            if  printSummary   println(iostream, sa)    end
            #
            for  confa  in confs
                print("  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")
                if  printSummary   println(iostream, "\n*  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")   end
                    basis     = Basics.performSCF([confa], comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                    multiplet = Basics.performCI(basis,    comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
            end
        else  error("Unsupported cascade approach.")
        end

        return( blockList )
    end

    """
    `Cascade.generateConfigurationsForDielectronicCapture(multiplets::Array{Multiplet,1},  scheme::HollowIonScheme,
                                                      nm::Nuclear.Model, grid::Radial.Grid)`
        ... generates all possible doubly-excited configurations due to (dielectronic) electron capture into the given multiplets.
            The number and type of such doubly-generated configurations depend on (1) the maximum (electron) energy for capturing an electron
            that is closely related to the (maximum) temperature of the plasma; (2) the fromShells from which (and how many displacements)
            are accepted as well as (3) the maximum principle and orbital angular quantum number of the additional (to-) shellsthe fromShells
            into which electrons excited and/or captured. A Tuple(initialConfList::Array{Configuration,1}, confList::Array{Configuration,1})
            is returned.
    """
    function generateConfigurationsForDielectronicCapture(multiplets::Array{Multiplet,1},  scheme::HollowIonScheme,
                                                      nm::Nuclear.Model, grid::Radial.Grid)
        # Determine all (reference) configurations from multiplets and generate the 'excited' configurations due to the specificed excitations
        initialConfList = Configuration[]
        for mp  in  multiplets
            confList = Basics.extractNonrelativisticConfigurations(mp.levels[1].basis)
            for  conf in confList   if  conf in initialConfList   nothing   else   push!(initialConfList, conf)      end      end
        end
        captureConfList = copy(initialConfList)
        for nc = 1:scheme.NoCapturedElectrons
            captureConfList = Basics.generateConfigurationsWithElectronCapture(captureConfList, Shell[], scheme.intoShells, 0);   @show nc, captureConfList
        end
        shellList       = Basics.extractNonrelativisticShellList(multiplets)
        shellList       = Basics.merge(scheme.intoShells, scheme.decayShells);  @show shellList
        blockConfList   = copy(initialConfList)
        for nc = 1:scheme.NoCapturedElectrons
            blockConfList   = Basics.generateConfigurationsWithElectronCapture(blockConfList, Shell[], shellList, 0);   @show nc, blockConfList
        end
        @show initialConfList
        @show blockConfList
        @warn("blockConfList does not yet support the autoionization of the high-n electrons.")
        #
        # Determine first a hydrogenic spectrum for all subshells of the initial and doubly-excited states
        allConfList   = Configuration[];      append!(allConfList, initialConfList);      append!(allConfList, blockConfList)
        allSubshells  = Basics.extractRelativisticSubshellList(allConfList)
        primitives    = Bsplines.generatePrimitives(grid)
        orbitals      = Bsplines.generateOrbitalsHydrogenic(primitives, nm, allSubshells, printout=true)
        # Exclude configurations with too high mean energies
        en            = Float64[];
        for conf in initialConfList
            wen = Basics.determineMeanEnergy(conf, orbitals, nm, grid)
            push!(en, wen)
        end
        initialMean = sum(en) / length(en)
        println(">>> initial configuration(s) have mena energies  $initialMean  [a.u.].")

        return( (initialConfList, blockConfList)  )
    end

    """
    `Cascade.perform(scheme::HollowIonScheme, comp::Cascade.Computation)`
        ... to set-up and perform a dielectronic-recombination (DR) plasma rate coefficient computation that combines the electron capture
            and the subsequent radiative stabilization steps. Such a computation starts from a given set of initial configurations xor
            initial multiplets and (1) generates all doubly-excited configurations due to the capture of an electron with a given maximum
            electron energy; (2) selects all electron capture (inverse Auger) and re-autoionization (Auger) steps and (3) selects
            all steps for radiative stabilization due to given parameters of the scheme::HollowIonScheme. The results of
            these DR plasma rate computation are comprised into (output) data::ExcitationData, while these data are only printed during
            the generation and nothing is returned.

    `Cascade.perform(scheme::HollowIonScheme, comp::Cascade.Computation; output=true, outputToFile::Bool=true)`
        ... to perform the same but to return the complete output in a dictionary that is written to disk and can be used in subsequent
            cascade simulation. The particular output depends on the specifications of the cascade.
    """
    function perform(scheme::HollowIonScheme, comp::Cascade.Computation; output::Bool=false, outputToFile::Bool=true)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        # Perform the SCF and CI computation for the intial-state multiplets if initial configurations are given
        if  comp.initialConfigs != Configuration[]
            basis      = Basics.performSCF(comp.initialConfigs, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplet  = Basics.performCI(basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplets = [Multiplet("initial states", multiplet.levels)]
        else
            multiplets = comp.initialMultiplets
        end
        # Print out initial configurations and levels
        Cascade.displayLevels(stdout, multiplets, sa="initial ")
        if  printSummary   Cascade.displayLevels(iostream, multiplets, sa="initial ")                            end
        #
        # Generate subsequent cascade configurations as well as display and group them together
        wa  = Cascade.generateConfigurationsForDielectronicCapture(multiplets, comp.scheme, comp.nuclearModel, comp.grid)
        ##x wb1 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[1], sa="(initial part of the) electron capture ")
        wb2 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[2], sa="hollow ion configurations ")
        #
        # Determine first all configuration 'blocks' and from them the individual steps of the cascade
        ##x wc1 = Cascade.generateBlocks(scheme, comp::Cascade.Computation, wb1)
        wc2 = Cascade.generateBlocks(scheme, comp::Cascade.Computation, wb2, printout=false)
        ##x Cascade.displayBlocks(stdout, wc1, sa="for the (initial part of the) excited cascade ");
        Cascade.displayBlocks(stdout, wc2, sa="for the hollow ion cascade ")
        if  printSummary   ##x Cascade.displayBlocks(iostream, wc1, sa="for the (initial part of the) excited cascade ")
                           Cascade.displayBlocks(iostream, wc2, sa="for the hollow ion cascade ")    end
        #
        # Determine, modify and compute the transition data for all steps, ie. the PhotoIonization.Line's, etc.
        gMultiplets = Multiplet[];     for block in wc2  push!(gMultiplets, block.multiplet)    end
        ##x we = Cascade.determineSteps(scheme, comp, wc1, wc2)
        we = Cascade.determineSteps(scheme, comp, wc2)
        Cascade.displaySteps(stdout, we, sa="hollow ion decay ")
        if  printSummary   Cascade.displaySteps(iostream, we, sa="hollow ion decay ")    end
        wf   = Cascade.modifySteps(we)
        #
        data = Cascade.computeSteps(scheme, comp, wf)
        if output
            results = Base.merge( results, Dict("name"                                  => comp.name) )
            results = Base.merge( results, Dict("cascade scheme"                        => comp.scheme) )
            results = Base.merge( results, Dict("initial multiplets:"                   => multiplets) )
            results = Base.merge( results, Dict("generated multiplets:"                 => gMultiplets) )
            results = Base.merge( results, Dict("hollow-ion line data:" => data) )
            #
            #  Write out the result to file to later continue with simulations on the cascade data
            filename = "zzz-cascade-hollow-ion-computations-" * string(Dates.now())[1:13] * ".jld"
            println("\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                    "\n   results = JLD.load(''$filename'')    ... to load the results back from file.")
            if  printSummary   println(iostream, "\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                                                 "\n   results = JLD.load(''$filename'')    ... to load the results back from file." )      end
            JLD2.@save filename results
        end
        ## return( results )
        return( results )
    end


text:
################################################################################
#
#        NfOrdClsZeta.jl: Zeta functions of orders in number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export zeta_log_residue

################################################################################
#
# Macros to make floating point life (hopefully) easier
#
################################################################################

macro with_round_down(expr, T)
  quote
    setrounding($(esc(T)), RoundDown) do
      $(esc(expr))
    end
  end
end

macro with_round_up(expr, T)
  quote
    setrounding($(esc(T)), RoundUp) do
      $(esc(expr))
    end
  end
end

################################################################################
#
#  Helper
#
################################################################################

# compute largest m such that a^m < b (strictly less then!)

function _max_power_in(a::fmpz, b::Int)
  m = flog(fmpz(b), a)
  if a^m == b
    m = m - 1
  end
  return m
end

_max_power_in(a::Int, b::Int) = _max_power_in(fmpz(a), b)

################################################################################
#
#  Residue approximation apres
#  Belabas, Friedmann, "Computing the residue of the Dedekind zeta function"
#
################################################################################

# bounding the error
function _approx_error_bf(O::NfOrd, Tc = BigFloat)
  return _approx_error_bf(discriminant(O), degree(O), Tc)
end

# This function gives an upper bound on the error term of the Belabas-Friedmann
# approximation.
function _approx_error_bf(disc::fmpz, degree::Int, Tc = BigFloat)

  logd_up = Tc(0)::Tc
  logd_down = Tc(0)::Tc
  sqrt_logd_up = Tc(0)::Tc

  setrounding(Tc,RoundDown) do
    logd_down = log(Tc(abs(disc)))
  end

  setrounding(Tc,RoundUp) do
    logd_up = log(Tc(abs(disc)))
    sqrt_logd_up = sqrt(logd_up)
  end

  n = BigFloat(degree)

  C1 = @with_round_down(Tc(FlintQQ(2324)//FlintQQ(1000)),Tc)
  C2 = @with_round_down(Tc(FlintQQ(388)//FlintQQ(100)),Tc)
  C3 = Tc(2)
  C4 = @with_round_down(Tc(FlintQQ(426)//FlintQQ(100)),Tc)

  function F(X)#(X::Tc)
    A1 = @with_round_down(C1*logd_down/(@with_round_up(sqrt(X)*log(3*X),Tc)),Tc)
    A2 = @with_round_down(1 + C2/@with_round_up(log(X/9),Tc),Tc)
    A3 = @with_round_down(1 + 2/sqrt_logd_up,Tc)
    A4 = @with_round_down(C4*(n-1)/@with_round_up(sqrt(X)*logd_up,Tc),Tc)

    return A1*(A2*A3^2 + A4)
  end
  return F
end

# Given f and C, this function computes a small x such that f(x) < C.
# It is assumed that f is monotone
function _find_threshold(f, C, ste, decreasing::Bool, Tc = BigFloat)
  T = Tc
  x0 =  T(70)
  x1 = x0

  y = f(x0)

  while y > C
    x0 = x1
    x1 = 2*x0
    y = f(x1)
  end

  dista = abs(x0-x1)

  while !( y < C && dista < ste)
    if y > C
      x1 = x0 + 3*dista/2
    else
      x1 = x0 - dista/2
    end

    dista = abs(x1-x0)

    x0 = x1
    y = f(x0)
  end
  return x0::Tc
end

  function _comp_summand(R, p::fmpz, m::Int, aa::arb)
    logp = log(R(p))

    pm2 = R(p)^(R(FlintZZ(m)//FlintZZ(2)))

    pm2inv = inv(pm2)

    pro = logp * pm2inv

    # now p^(m/2) m log(p)

    pro2 = logp*pm2
    pro2 = pro2*m

    # Now the inverse
    inv2 = inv(pro2)

    # Now sqrt(x)log(X)/p^(m/2)*m*p
    pro3 = aa*inv2
    pro3 = pro3 - 1

    return pro*pro3
  end

  function _comp_summand(R, p::Int, m::Int, aa::arb)
    return _comp_summand(R, fmpz(p), m, aa)
  end

# Computing the g_K(X) term of Belabas-Friedmann
function _term_bf(O::NfOrd, B::Int, R::ArbField)

  xx0 = B

  # I want xx0 to be divisible by 9!
  # Because I want xx0/9 to be an integer

  while !(mod(xx0,9) == 0)
    xx0 += 1
  end

  x = Float64(0)

  xx09 = div(xx0,9)

  p = 2

  summ = R(0)
  esumm = R(0)

  logxx0 = log(R(xx0))
  logxx09 = log(R(xx09))
  sqrtxx0 = sqrt(R(xx0))
  sqrtxx09 = sqrt(R(xx09))
  prodx = logxx0 * sqrtxx0
  prodx9 = logxx09 * sqrtxx09

  # small helper function (is this fast?)

  while p < xx0

    max_exp = _max_power_in(p, xx0)

    #println("maximal power is $max_exp")

    for m in 1:max_exp
      summand = _comp_summand(R, p, m, prodx)
      summ = summ - summand
    end

    #x += @elapsed

    lP = prime_decomposition_type(O, p)

    for P in lP
      Pnorm = fmpz(p)^P[1]
      if Pnorm < xx0
        max_exp = _max_power_in(Pnorm, xx0)

        for m in 1:max_exp
          summand = _comp_summand(R, Pnorm, m, prodx)
          summ = summ + summand
        end
      end
    end

    if p < xx09

      max_exp = _max_power_in(p, xx09)

      for m in 1:max_exp
        summand  = _comp_summand(R, p, m, prodx9)
        summ = summ + summand
      end

      for P in lP
        Pnorm = fmpz(p)^P[1]
        if (Pnorm < xx09)
          max_exp = _max_power_in(Pnorm, xx09)

          for m in 1:max_exp
            summand = _comp_summand(R, Pnorm, m, prodx9)
            summ = summ - summand
          end
        end
      end
    end
    p = next_prime(p)
  end

  log3xx0 = log(R(3*xx0))
  pr = sqrtxx0*2
  pr = pr * log3xx0
  pr = inv(pr)
  pr = pr*3
  pr = pr*summ

  return pr
end

# Approximate the residue
function _residue_approx_bf(O::NfOrd, error::Float64)
  F = _approx_error_bf(O, BigFloat)

  # magic constant
  # should be adapted to the input

  der = Int(20)

  @assert error > 0.5^der

  error_prime = @with_round_down(error - 0.5^der, BigFloat)

  error_prime = Float64(error_prime) # Do some correct rounding

  error_prime_arf = arf_struct(0, 0, 0, 0)
  ccall((:arf_init, :libarb), Nothing, (Ref{arf_struct}, ), error_prime_arf)
  ccall((:arf_set_d, :libarb), Nothing, (Ref{arf_struct}, Float64), error_prime_arf, error_prime)

  error_arf = arf_struct(0, 0, 0, 0)
  ccall((:arf_init, :libarb), Nothing, (Ref{arf_struct}, ), error_arf)
  ccall((:arf_set_d, :libarb), Nothing, (Ref{arf_struct}, Float64), error_arf, error)

  x0 = Int(ceil(_find_threshold(F, error_prime, Float64(10), true, Float64)))
  x0 = x0 + 1

  prec = 64

  val = _term_bf(O, x0, ArbField(prec, false))

  valaddederror = deepcopy(val)
  ccall((:arb_add_error_arf, :libarb), Nothing,
              (Ref{arb}, Ref{arf_struct}), valaddederror, error_prime_arf)

  while (!radiuslttwopower(val, -der)) ||
                !((radius(valaddederror)) < error)

    #println("Precision is now $prec")

    if prec > 2048
      error("Something wrong")
    end

    prec = 2*prec
    #println("increasing precision to $prec")
    val = _term_bf(O, x0, ArbField(prec, false))
    valaddederror = deepcopy(val)
    ccall((:arb_add_error_arf, :libarb), Nothing,
                (Ref{arb}, Ref{arf_struct}), valaddederror, error_prime_arf)
  end

  ccall((:arf_clear, :libarb), Nothing, (Ref{arf_struct}, ), error_prime_arf)
  ccall((:arf_clear, :libarb), Nothing, (Ref{arf_struct}, ), error_arf)

  return valaddederror
end

################################################################################
#
#  Toplevel function for users
#
################################################################################

@doc Markdown.doc"""
    zeta_log_residue(O::NfOrd, error::Float64) -> arb

Computes the residue of the zeta function of $\mathcal O$ at $1$.
The output will be an element of type `arb` with radius less then
`error`.
"""
function zeta_log_residue(O::NfOrd, abs_error::Float64)
  degree(O) == 1 && error("Number field must be of degree 1")
  return _residue_approx_bf(O, abs_error)
end

ext: jl


text:
using VertexEliminationOrder
using Test

@testset "VertexEliminationOrder.jl" begin
    include("test_utils.jl")
    include("test_heuristics.jl")
    include("test_flowcutter.jl")
    include("test_dissection.jl")

end


text:
using AtmosphericProfilesLibrary, Documenter
using DocumenterCitations
import Glob

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))

include("generate_profile_library.jl")

#! format: off
pages = Any[
    "Home" => "index.md",
    "Profile library" => "generated_profile_library.md",
    "References" => "References.md",
]

mathengine = MathJax(Dict(
    :TeX => Dict(
        :equationNumbers => Dict(:autoNumber => "AMS"),
        :Macros => Dict(),
    ),
))

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
)
#! format: on

makedocs(
    bib,
    sitename = "AtmosphericProfilesLibrary.jl",
    strict = true,
    format = format,
    checkdocs = :exports,
    clean = true,
    doctest = true,
    modules = [AtmosphericProfilesLibrary],
    pages = pages,
)

deploydocs(
    repo = "github.com/CliMA/AtmosphericProfilesLibrary.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)

# Remove temp files (e.g., .DS_Store)
for (root, _, files) in Base.Filesystem.walkdir(pkgdir(AtmosphericProfilesLibrary))
    for f in files
        file = joinpath(root, f)
        if endswith(file, ".DS_Store")
            rm(file; force = true)
        end
        occursin("build/", file) && continue
        if endswith(file, ".png")
            rm(file; force = true)
        end
    end
end


text:
#########################################################################
#    testing script for simple examples
#########################################################################

using Base.Test
using Distributions
using MCMC

# generate a random dataset
srand(1)
n = 1000
nbeta = 10
X = [ones(n) randn((n, nbeta-1))]
beta0 = randn((nbeta,))
Y = rand(n) .< ( 1 ./ (1 .+ exp(X * beta0)))

# define model
ex = quote
	vars ~ Normal(0, 1.0)
	prob = 1 ./ (1. .+ exp(X * vars))
	Y ~ Bernoulli(prob)
end

m = model(ex, vars=zeros(nbeta), gradient=true)

# different samplers
res = run(m * RWM(0.05) * SerialMC(100:1000))
res = run(m * HMC(2, 0.1) * SerialMC(100:1000))
res = run(m * NUTS() * SerialMC(100:1000))
res = run(m * MALA(0.001) * SerialMC(100:1000))
# TODO : add other samplers here

# different syntax
res = run(m, RWM(), SerialMC(steps=1000, thinning=10, burnin=0))
res = run(m, HMC(2,0.1), SerialMC(thinning=10, burnin=0))
res = run(m, HMC(2,0.1), SerialMC(burnin=20))

### README examples

mymodel1 = model(v-> -dot(v,v), init=ones(3))
mymodel2 = model(v-> -dot(v,v), grad=v->-2v, init=ones(3))

modelxpr = quote
    v ~ Normal(0, 1)
end

mymodel3 = model(modelxpr, v=ones(3))
mymodel4 = model(modelxpr, gradient=true, v=ones(3))

mychain = run(mymodel1, RWM(0.1), SerialMC(steps=1000, burnin=100))
mychain = run(mymodel1, RWM(0.1), SerialMC(steps=1000, burnin=100, thinning=5))
mychain = run(mymodel1, RWM(0.1), SerialMC(101:5:1000))
mychain1 = run(mymodel1 * RWM(0.1) * SerialMC(101:5:1000))

mychain2 = run(mymodel2, HMC(0.75), SerialMC(steps=10000, burnin=1000))

acceptance(mychain2)

# describe(mychain2)

ess(mychain2)

actime(mychain2)

# var(mychain2)
# var(mychain2, vtype=:iid)
# var(mychain2, vtype=:ipse)
# var(mychain2, vtype=:bm)

mychain1 = resume(mychain1, steps=10000)

@test_throws ErrorException run(mymodel3 * MALA(0.1) * SerialMC(1:1000))

run(mymodel4 * MALA(0.1) * SerialMC(1:1000))

mychain = run(mymodel2 * [RWM(0.1), MALA(0.1), HMC(3,0.1)] * SerialMC(steps=1000))
mychain[2].samples

mychain = run(mymodel2 * [HMC(i,0.1) for i in 1:5] * SerialMC(steps=1000))

nmod = 10
mods = Array(MCMCLikModel, nmod)
sts = logspace(1, -1, nmod)
for i in 1:nmod
  m = quote
    y = abs(x)
    y ~ Normal(1, $(sts[i]))
  end
  mods[i] = model(m, x=0)
end

targets = MCMCTask[mods[i] * RWM(sts[i]) * SeqMC(steps=10, burnin=0) for i in 1:nmod]
particles = [[randn()] for i in 1:1000]

mychain3 = run(targets, particles=particles)

# mychain4 = wsample(mychain3.samples, mychain3.diagnostics["weigths"], 1000)
# mean(mychain4)


text:
using MixedAnova, CSV, RDatasets, DataFrames, CategoricalArrays
using Test
import Base.isapprox, MixedAnova.dof_residual

test_show(x) = show(IOBuffer(), x)
macro test_error(x)
    return quote
        try
            $x
            false
        catch e
            @error e
            true
        end
    end
end

glm_init()
mm_init()
fem_init()

const anova_datadir = joinpath(dirname(@__FILE__), "..", "data")

"Data from R datasets"
iris = dataset("datasets", "iris")

"Data source: https://www.qogdata.pol.gu.se/dataarchive/qog_bas_cs_jan18.dta"
qog18 = CSV.read(joinpath(anova_datadir, "qog18.csv"), DataFrame)

"Data source: https://www.sciencedirect.com/science/article/abs/pii/0160289691900318"
iq = CSV.read(joinpath(anova_datadir, "iqsize.csv"), DataFrame)

"Data from R package MASS"
quine = dataset("MASS", "quine")

"Data from R datasets"
mtcars = dataset("datasets", "mtcars")
transform!(mtcars, [:VS, :Model] .=> categorical, renamecols = false)

"Data source: https://www.cambridge.org/core/journals/astin-bulletin-journal-of-the-iaa/article/two-studies-in-automobile-insurance-ratemaking/1A169017AE99FCFE6E5F6BF56BCCA9C2"
insurance = CSV.read(joinpath(anova_datadir, "automobile_insurance.csv"), DataFrame)
transform!(insurance, [1, 2] .=> categorical, renamecols = false)

"Data source: https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta"
nb = CSV.read(joinpath(anova_datadir, "nb_data.csv"), DataFrame)
transform!(nb, [1, 2, 5] .=> categorical, renamecols = false)

"Data source: https://stats.idre.ucla.edu/stat/stata/dae/poisson_sim"
sim = CSV.read(joinpath(anova_datadir, "poisson_sim.csv"), DataFrame)
transform!(sim, [:id, :prog] .=> categorical, renamecols = false)

"Examples from https://m-clark.github.io/mixed-models-with-R/"
gpa = CSV.read(joinpath(anova_datadir, "gpa.csv"), DataFrame)
transform!(gpa,
    7 => x->replace(x, "yes" => true, "no" => false, "NA" => missing),
    4 => x->categorical(x, levels = ["1 hour", "2 hours", "3 hours"], ordered = true),
    renamecols = false)
transform!(gpa, [1, 2, 5, 7] .=> categorical, renamecols = false)

nurses = CSV.read(joinpath(anova_datadir, "nurses.csv"), DataFrame)
transform!(nurses,
    [1, 2, 3, 4, 6, 9, 11] .=> categorical,
    10 => x->categorical(x, levels = ["small", "medium", "large"], ordered = true),
    renamecols = false)

"Data from R package datarium"
anxiety = CSV.read(joinpath(anova_datadir, "anxiety.csv"), DataFrame)
transform!(anxiety, :id => categorical, renamecols = false)

"Data from R package HSAUR2"
toenail = CSV.read(joinpath(anova_datadir, "toenail.csv"), DataFrame)
transform!(toenail, [1, 2, 3] .=> categorical, renamecols = false)

# custimized approx
isapprox(x::NTuple{N, Float64}, y::NTuple{N, Float64}, atol::NTuple{N, Float64} = x ./ 1000) where N =
    all(map((a, b, c)->isapprox(a, b, atol = c > eps(Float64) ? c : eps(Float64)), x, y, atol))

dof_residual(x::Int) = x
@testset "LinearModel" begin
    @testset "Simple linear regression" begin
        lm0, lm1, lm2, lm3, lm4 = nestedmodels(LinearModel, @formula(SepalLength ~ SepalWidth * Species), iris, dropcollinear = false)
        global aov3 = anova(lm4, type = 3)
        global aov2 = anova_lm(@formula(SepalLength ~ SepalWidth * Species), iris, type = 2)
        global aov1 = anova(lm4)
        global aovf = anova(lm0, lm1, lm2, lm3, lm4)
        global aovlr = anova(LRT, lm0, lm1, lm2, lm3, lm4)
        global aov1lr = anova(LRT, lm4)
        global aovlf = anova_lm(FTest, @formula(wdi_lifexp ~ log(gle_rgdpc) * ti_cpi), qog18, type = 2)
        ft = ftest(lm1.model, lm2.model, lm3.model, lm4.model)
        @test !(@test_error(test_show(aov1)))
        @test !(@test_error(test_show(aovf)))
        @test !(@test_error(test_show(aov1lr)))
        @test !(@test_error(test_show(aovlr)))
        @test !(@test_error(test_show(aovlf)))
        @test anova_type(aov1) == 1
        @test nobs(aov2) == ft.nobs
        @test dof(aov3) == (1, 1, 2, 2, 144)
        @test dof_residual(aov1) == 144
        @test isapprox(deviance(aovf)[2:end], ft.ssr)
        @test isapprox(deviance(aov1)[1:end - 1], MixedAnova._diffn(deviance(aovf)))
        @test isapprox(filter(!isnan, teststat(aov1)), filter(!isnan, teststat(aovf)))
        @test isapprox(filter(!isnan, teststat(aov2)), (26485.300978452644, 56.637879295914324, 188.10911669921464, 0.4064420847481629))
        @test isapprox(filter(!isnan, teststat(aov1lr)), filter(!isnan, teststat(aovlr)))
        @test isapprox(coef(lm4), coef(aov3.model))
        @test coefnames(lm4, Val(:anova)) ==  ["(Intercept)", "SepalWidth", "Species", "SepalWidth & Species"]
        @test collect(coefnames(formula(lm1), Val(:anova))) == coefnames((formula(lm1).lhs, formula(lm1).rhs), Val(:anova))

    end

    @testset "Linear regression with frequency weights" begin
        wlm1, wlm2 = nestedmodels(LinearModel, @formula(Cost / Claims ~ Insured + Merit), insurance,
                                wts = insurance.Claims ./ sum(insurance.Claims) .* length(insurance.Claims))
        global aov = anova(wlm2)
        global aovf = anova(wlm1, wlm2)
        @test !(@test_error test_show(aov))
        @test !(@test_error test_show(aovf))
        @test nobs(aov) == sum(aov.model.model.rr.wts)
        @test dof(aov) == (1, 1, 18)
        @test isapprox(filter(!isnan, teststat(aov))[2:end], filter(!isnan, teststat(aovf)))
    end
end

@testset "GeneralizedLinearModel" begin
    @testset "Gamma regression" begin
        global aov = anova_glm(FTest, @formula(Cost / Claims ~ 0 + Insured + Merit), insurance, Gamma(),
                                wts = insurance.Claims ./ sum(insurance.Claims) .* length(insurance.Claims), type = 2)
        @test !(@test_error test_show(aov))
        @test nobs(aov) == nobs(aov.model)
        @test dof(aov) == (1, 4, 15)
        @test isapprox(filter(!isnan, deviance(aov)), MixedAnova.deviances(aov.model, type = anova_type(aov)))
        @test isapprox(filter(!isnan, teststat(aov)), (6.163653078060339, 2802.3252386290533))
        @test isapprox(filter(!isnan, pval(aov)), (0.025357816283854216, 2.3626325930920802e-21))
    end

    @testset "NegativeBinomial regression" begin
        global aov = anova_glm(@formula(Days ~ Eth + Sex + Age + Lrn), quine, NegativeBinomial(2.0), LogLink(), type = 3)
        @test !(@test_error test_show(aov))
        @test nobs(aov) == nobs(aov.model)
        @test dof(aov) == (1, 1, 1, 3, 1, 139)
        @test anova_test(aov) == FTest
        @test isapprox(deviance(aov), MixedAnova.deviances(aov.model, type = 3))
        @test isapprox(filter(!isnan, teststat(aov)), (227.97422313423752, 13.180680587112887, 0.2840882132754838, 4.037856229143672, 2.6138558930314595))
        @test isapprox(filter(!isnan, pval(aov)), (4.251334236308285e-31, 0.00039667906264309605, 0.5948852219093326, 0.008659894572621351, 0.10820118567663468))
        @test coefnames(aov.model, Val(:anova)) == ["(Intercept)", "Eth", "Sex", "Age", "Lrn"]
    end

    @testset "Poisson regression" begin
        gmp = nestedmodels(GeneralizedLinearModel, @formula(num_awards ~ prog * math), sim, Poisson())
        global aov = anova(gmp...)
        lr = lrtest(gmp...)
        @test !(@test_error test_show(aov))
        @test first(nobs(aov)) == lr.nobs
        @test dof(aov) == lr.dof
        @test anova_test(aov) == LRT
        @test isapprox(deviance(aov), lr.deviance)
        @test isapprox(filter(!isnan, pval(aov)), filter(!isnan, lr.pval))
    end

    @testset "Logit regression" begin
        gml = glm(@formula(AM ~ Cyl + HP + WT), mtcars, Binomial(), LogitLink())
        global aov = anova(gml)
        lr = lrtest(nestedmodels(gml)...)
        @test !(@test_error test_show(aov))
        @test nobs(aov) == lr.nobs
        @test dof(aov) == MixedAnova._diff(lr.dof)
        @test isapprox(deviance(aov), lr.deviance[2:end])
        @test isapprox(MixedAnova.deviances(aov.model)[2:end], MixedAnova._diffn((deviance(aov)..., 0.0)))
        @test isapprox(filter(!isnan, pval(aov)), filter(!isnan, lr.pval))
    end

    @testset "Probit regression" begin
        gmp0 = glm(@formula(AM ~ 1), mtcars, Binomial(), ProbitLink())
        gmp1 = glm(@formula(AM ~ Cyl + HP + WT), mtcars, Binomial(), ProbitLink())
        gmp2 = glm(@formula(AM ~ Cyl * HP * WT), mtcars, Binomial(), ProbitLink())
        global aov = anova(gmp0, gmp1, gmp2)
        lr = lrtest(gmp0, gmp1, gmp2)
        @test !(@test_error test_show(aov))
        @test first(nobs(aov)) == lr.nobs
        @test dof(aov) == lr.dof
        @test isapprox(deviance(aov), lr.deviance)
        @test isapprox(filter(!isnan, pval(aov)), filter(!isnan, lr.pval))
    end

    @testset "InverseGaussian regression" begin
        gmi = glm(@formula(SepalLength ~ SepalWidth * Species), iris, InverseGaussian())
        global aov = anova(gmi)
        @test !(@test_error test_show(aov))
        @test nobs(aov) == nobs(gmi)
        @test dof(aov) == (1, 2, 2, 144)
        @test isapprox(filter(!isnan, teststat(aov)), (8.172100334461327, 217.34941014323272, 1.8933247444892272))
        @test isapprox(filter(!isnan, pval(aov)), (0.004885873691352542, 3.202701170052312e-44, 0.15429963193830074))
    end
end

# test contrast
@testset "LinearMixedModel" begin
    @testset "One random effect on intercept" begin
        lmm1 = lme(@formula(gpa ~ occasion + sex + job + (1|student)), gpa)
        lmm2 = lme(@formula(gpa ~ occasion * sex + job + (1|student)), gpa)
        lmm3 = lme(@formula(gpa ~ occasion * sex * job + (1|student)), gpa)
        lme3 = anova_lme(@formula(score ~ group * time + (1|id)), anxiety, type = 3)
        lme1 = anova_lme(LRT, @formula(score ~ group * time + (1|id)), anxiety)
        global aovf = anova(lmm1)
        global aovlr = anova(lmm1, lmm2, lmm3)
        global aov3 = anova(lmm1, type = 3)
        lr = MixedModels.likelihoodratiotest(lmm1, lmm2, lmm3)
        @test !(@test_error test_show(aovf))
        @test !(@test_error test_show(aovlr))
        @test !(@test_error test_show(aov3))
        @test anova_type(aov3) == 3
        @test first(nobs(aovlr)) == nobs(lmm1)
        @test dof(aovf) == (1, 1, 1, 2)
        @test dof_residual(aovf) == (997, 997, 198, 997)
        @test isapprox(teststat(aovf), (28899.81231026455, 714.3647106859062, 19.19804098862722, 45.425989531937134))
        @test isapprox(deviance(aovlr), tuple(lr.deviance...))
        @test isapprox(filter(!isnan, pval(aovlr)), tuple(lr.pvalues...))
        @test lme3.model.optsum.REML
        @test !last(lme1.model).optsum.REML
        @test coefnames(lmm2, Val(:anova)) == ["(Intercept)",  "occasion", "sex", "job", "occasion & sex"]
    end
    @testset "Random effect on slope and intercept" begin
        lmm1 = lme(@formula(gpa ~ occasion + (1|student)), gpa, REML = true)
        lmm2 = lme(@formula(gpa ~ occasion + (occasion|student)), gpa, REML = true)
        lms = nestedmodels(LinearMixedModel, @formula(gpa ~ occasion + (1|student)), gpa, REML = true)
        global aovlr = anova(lmm1, lmm2)
        global aovf = anova(lmm2)
        lr = MixedModels.likelihoodratiotest(lmm1, lmm2)
        @test !(@test_error test_show(aovlr))
        @test !(@test_error test_show(aovf))
        @test first(nobs(aovlr)) == nobs(lmm1)
        @test dof(aovlr) == tuple(lr.dof...)
        @test dof_residual(aovf) == (1000, 198)
        @test isapprox(deviance(aovlr), tuple(lr.deviance...))
        @test isapprox(filter(!isnan, pval(aovlr)), tuple(lr.pvalues...))
    end
    @testset "Multiple random effects" begin
        lmm1 = lme(@formula(stress ~ age  + sex + experience + treatment + wardtype + hospsize + (1|hospital)), nurses)
        lmm2 = lme(@formula(stress ~ age  + sex + experience + treatment + wardtype + hospsize + (1|hospital/wardid)), nurses)
        global aov = anova(lmm1, lmm2)
        lr = MixedModels.likelihoodratiotest(lmm1, lmm2)
        @test !(@test_error test_show(aov))
        @test first(nobs(aov)) == nobs(lmm1)
        @test dof(aov) == tuple(lr.dof...)
        @test isapprox(deviance(aov), tuple(lr.deviance...))
        @test isapprox(filter(!isnan, pval(aov)), tuple(lr.pvalues...))
    end
    #=
    @testset "Random effects on slope and intercept" begin
        aov1 = anova_lme(@formula(extro ~ open + agree + social + (1|school) + (1|class)), school, REML = true)
        aov2 = anova_lme(@formula(extro ~ open + agree + social + (open|school) + (open|class)), school, REML = true)
        @test aov1.stats.type == 1
        @test aov2.stats.nobs == 1200
        @test all(aov1.stats.dof .== (1, 1, 1, 1))
        @test all(aov2.stats.resdof .== (1192, 2, 1192, 1192))
        @test isapprox(aov1.stats.fstat, (208.372624411485, 1.673863453540288, 0.3223490172035519, 0.32161077338130784))
        @test isapprox(aov2.stats.pval, (1.5264108749941817e-43, 0.32443726673109696, 0.5622508640865989, 0.5681019984771887))
    end

    @testset "Nested random effects on intercept" begin
        aov = anova_lme(@formula(extro ~ open + agree + social + (1|school) + (1|school&class)), school)
        @test aov.stats.type == 1
        @test aov.stats.nobs == 1200
        @test all(aov.stats.dof .== (1, 1, 1, 1))
        @test all(aov.stats.resdof .== (1173, 1173, 1173, 1173))
        @test isapprox(aov.stats.fstat, (227.24904612916632, 1.4797389900418527, 1.8026931823189165, 0.08511022388038622))
        @test isapprox(aov.stats.pval, (4.508653887500674e-47, 0.22406007872324163, 0.1796470040402073, 0.7705396333261272))
    end

    @testset "Nested random effects on slope" begin
        aov = anova_lme(@formula(extro ~ open + agree + social + (open|school) + (open|school&class)), school)
        @test aov.stats.type == 1
        @test aov.stats.nobs == 1200
        @test all(aov.stats.dof .== (1, 1, 1, 1))
        @test all(aov.stats.resdof .== (1174, 4, 1174, 1174))
        @test isapprox(aov.stats.fstat, (250.00542522864583, 1.2322678515772565, 2.1135395635863543, 0.10258998684862923))
        @test isapprox(aov.stats.pval, (3.364459604379112e-51, 0.3292014599294774, 0.1462685388135273, 0.748800408618393))
    end
    =#
end

@testset "LinearModel and LinearMixedModel" begin
    lm1 = lm(@formula(score ~ group * time), anxiety)
    lmm1 = lme(@formula(score ~ group * time + (1|id)), anxiety)
    lmm2 = lme(@formula(score ~ group * time + (group|id)), anxiety)
    global aov = anova(lm1, lmm1, lmm2)
    lr = MixedModels.likelihoodratiotest(lm1, lmm1, lmm2)
    @test !(@test_error test_show(aov))
    @test first(nobs(aov)) == nobs(lmm1)
    @test dof(aov) == tuple(lr.dof...)
    @test isapprox(deviance(aov), tuple(lr.deviance...))
    @test isapprox(pval(aov)[2:end], tuple(lr.pvalues...))
end

@testset "GeneralizedLinearModel and GeneralizedLinearMixedModel" begin
    glm1 = glm(@formula(outcome ~ visit + treatment), toenail, Binomial(), LogitLink())
    glmm1 = glme(@formula(outcome ~ visit + treatment + (1|patientID)), toenail, Binomial(), LogitLink(), nAGQ=20, wts = ones(Float64, size(toenail, 1)))
    glmm2 = glme(@formula(outcome ~ visit * treatment + (1|patientID)), toenail, Binomial(), LogitLink(), nAGQ=20, wts = ones(Float64, size(toenail, 1)))
    global aov = anova(glm1, glmm1, glmm2)
    lr = MixedModels.likelihoodratiotest(glm1, glmm1, glmm2)
    @test !(@test_error test_show(aov))
    @test first(nobs(aov)) == nobs(glmm1)
    @test dof(aov) == tuple(lr.dof...)
    @test isapprox(deviance(aov), tuple(lr.deviance...))
    @test isapprox(pval(aov)[2:end], tuple(lr.pvalues...))
end
@testset "Miscellaneous" begin
    lm1 = lm(@formula(SepalLength ~ 0 + log(SepalWidth) * Species), iris, dropcollinear = false)
    lm2 = lm(@formula(SepalLength ~ SepalWidth + SepalWidth & Species), iris, dropcollinear = false)
    @test MixedAnova.subtablemodel(lm1, 2; reschema = true)[2].assign == [2]
    @test MixedAnova.subtablemodel(lm2, [2]; reschema = true)[2].assign == [1, 2, 2, 2]

    global ft = AnovaResult{FTest}(ntuple(identity, 7),
                            1,
                            ntuple(identity, 7),
                            ntuple(one ∘ float, 7),
                            ntuple(one ∘ float, 7),
                            ntuple(zero ∘ float, 7),
                            NamedTuple())
    global lrt = AnovaResult{LRT}(ntuple(identity, 7),
                            1,
                            ntuple(identity, 7),
                            ntuple(one ∘ float, 7),
                            ntuple(one ∘ float, 7),
                            ntuple(zero ∘ float, 7),
                            NamedTuple())
    @test @test_error test_show(ft)
    @test @test_error test_show(lrt)

    @test @test_error anova_glm(@formula(Cost / Claims ~ 0 + Class + Merit), insurance, Gamma(), type = 2)
    @test @test_error anova_glm(@formula(Cost / Claims ~ 0 + Merit), insurance, Gamma(), type = 2)
    @test @test_error anova_glm(@formula(Cost / Claims ~ 0 + Class + Merit), insurance, Gamma(), type = 3)
    @test @test_error anova_glm(@formula(Cost / Claims ~ 0 + Merit), insurance, Gamma(), type = 3)

    @test @test_error formula(1)
    @test @test_error nestedmodels(1)
    @test @test_error anova(1)
    @test @test_error anova(FTest, 1)
    @test @test_error anova(LRT, 1)
    @test @test_error MixedAnova.isnullable(1)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dash_graph3d

"""
    dash_graph3d(;kwargs...)

A Graph3D component.

Keyword arguments:
- `id` (String; required): The ID used to identify this component in Dash callbacks.
- `active` (Bool; optional): whether or not session is active. Used to enable or disable warning browser dialog when closing
- `backgroundColor` (String; optional): Getter/setter for the chart background color, default transparent
- `cameraPosition` (Array; optional): Re-position the camera, in terms of x, y, z coordinates. Each of the coordinates is optional, allowing for motion in just some dimensions. The optional second argument can be used to define the direction that the camera should aim at, in terms of an {x,y,z} point in the 3D space. The 3rd optional argument defines the duration of the transition (in ms) to animate the camera motion. A value of 0 (default) moves the camera immediately to the final position. By default the camera will face the center of the graph at a z distance proportional to the amount of nodes in the system. 3D
- `centerAt` (Array; optional): Set the coordinates of the center of the viewport. This method can be used to perform panning on the 2D canvas programmatically. Each of the x, y coordinates is optional, allowing for motion in just one dimension. An optional 3rd argument defines the duration of the transition (in ms) to animate the canvas motion.
- `centreCoordinates` (Dict with Strings as keys and values of type Real; optional): origin coordinates
- `controlType` (String; optional): Which type of control to use to control the camera on 3D mode. Choice between trackball, orbit or fly.
- `cooldownTicks` (Real; optional): How many build-in frames to render before stopping and freezing the layout engine.
- `cooldownTime` (Real; optional): How long (ms) to render for before stopping and freezing the layout engine.
- `d3AlphaDecay` (Real; optional): Sets the simulation intensity decay parameter. Only applicable if using the d3 simulation engine.
- `d3AlphaMin` (Real; optional): Sets the simulation alpha min parameter. Only applicable if using the d3 simulation engine.
- `d3ReheatSimulation` (Bool; optional): Reheats the force simulation engine, by setting the alpha value to 1. Only applicable if using the d3 simulation engine.
- `d3VelocityDecay` (Real; optional): Nodes' velocity decay that simulates the medium resistance. Only applicable if using the d3 simulation engine.
- `dagLevelDistance` (Real; optional): If dagMode is engaged, this specifies the distance between the different graph depths.
- `dagMode` (String; optional): Apply layout constraints based on the graph directionality. Only works correctly for DAG graph structures (without cycles). Choice between td (top-down), bu (bottom-up), lr (left-to-right), rl (right-to-left), zout (near-to-far), zin (far-to-near), radialout (outwards-radially) or radialin (inwards-radially).
- `dagModeOn` (Bool; optional): Apply layout constraints based on the graph directionality. Only works correctly for DAG graph structures (without cycles). Choice between td (top-down), bu (bottom-up), lr (left-to-right), rl (right-to-left), zout (near-to-far), zin (far-to-near), radialout (outwards-radially) or radialin (inwards-radially).
- `dagNodeIds` (Array of Strings; optional): array of string ids for nodes to include in DAG layout
- `emitParticle` (Dict; optional): An alternative mechanism for generating particles, this method emits a non-cyclical single particle within a specific link. The emitted particle shares the styling (speed, width, color) of the regular particle props. A valid link object that is included in graphData should be passed as a single parameter.
- `enableNodeDrag` (Bool; optional): Whether to enable the user interaction to drag nodes by click-dragging. If enabled, every time a node is dragged the simulation is re-heated so the other nodes react to the changes. Only applicable if enablePointerInteraction is true.
- `externalobject_input` (Real | String | Bool | Array | Dict; optional): externalobject_source:
- `externalobject_source` (String; optional): externalobject_source:
- `forceEngine` (String; optional): Which force-simulation engine to use (d3 or ngraph).
- `getGraphBbox` (Bool; optional): Returns the current bounding box of the nodes in the graph, formatted as { x: [<num>, <num>], y: [<num>, <num>], z: [<num>, <num>] }. If no nodes are found, returns null. Accepts an optional argument to define a custom node filter: node => <boolean>, which should return a truthy value if the node is to be included. This can be useful to calculate the bounding box of a portion of the graph.
Bounding box is saved as the graphBbox prop
- `graphData` (Dict; required): Graph data structure. Can also be used to apply incremental updates. Format {nodes:{}, links:{}}
- `gravity` (String; optional): gravity: not yet used, prop to change three gravity. Not used in 2D
- `heightRatio` (Real; optional): height of component as proportion of container
- `interactive` (Bool; optional): toggle enableZoomPanInteraction, enablePointerInteraction, enableNavigationControls with a single control
- `linkAutoColorBy` (String; optional): Link object accessor function or attribute to automatically group colors by. Only affects links without a color attribute.
- `linkColor` (String; optional): Link object accessor function or attribute for line color.
- `linkCurvature` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the curvature radius of the link line. A value of 0 renders a straight line. 1 indicates a radius equal to half of the line length, causing the curve to approximate a semi-circle. For self-referencing links (source equal to target) the curve is represented as a loop around the node, with length proportional to the curvature value.
- `linkCurveRotation` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the rotation along the line axis to apply to the curve. Has no effect on straight lines. At 0 rotation, the curve is oriented in the direction of the intersection with the XY plane. The rotation angle (in radians) will rotate the curved line clockwise around the "start-to-end" axis from this reference orientation.
- `linkDesc` (String; optional): For VR only. Link object accessor function or attribute for description (shown under label).
- `linkDirectionalArrowColor` (String; optional): Link object accessor function or attribute for the color of the arrow head.
- `linkDirectionalArrowLength` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the length of the arrow head indicating the link directionality. The arrow is displayed directly over the link line, and points in the direction of source > target. A value of 0 hides the arrow.
- `linkDirectionalArrowRelPos` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the longitudinal position of the arrow head along the link line, expressed as a ratio between 0 and 1, where 0 indicates immediately next to the source node, 1 next to the target node, and 0.5 right in the middle.
- `linkDirectionalArrowResolution` (Real; optional): Geometric resolution of the arrow head, expressed in how many slice segments to divide the cone base circumference. Higher values yield smoother arrows.
- `linkDirectionalParticleColor` (String; optional): Link object accessor function or attribute for the directional particles color.
- `linkDirectionalParticleResolution` (Real; optional): Geometric resolution of each 3D directional particle, expressed in how many slice segments to divide the circumference. Higher values yield smoother particles.
- `linkDirectionalParticleSpeed` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the directional particles speed, expressed as the ratio of the link length to travel per frame. Values above 0.5 are discouraged.
- `linkDirectionalParticleWidth` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the directional particles width. Values are rounded to the nearest decimal for indexing purposes.
- `linkDirectionalParticles` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the number of particles (small spheres) to display over the link line. The particles are distributed equi-spaced along the line, travel in the direction source > target, and can be used to indicate link directionality.
- `linkHoverPrecision` (Real; optional): Whether to display the link label when gazing the link closely (low value) or from far away (high value).
- `linkId` (String; optional): The link attribute containing the unique link id
- `linkIdsHighlight` (Array of Strings; optional): ids of highlighted links (through search)
- `linkIdsVisible` (Array of Strings; optional): ids of visible links
- `linkLabel` (String; optional): Link object accessor function or attribute for name (shown in label). Supports plain text or HTML content (except in VR).
- `linkLineDash` (Real | String; optional): Link object accessor function, attribute or number array (e.g. [5, 15]) to determine if a line dash should be applied to this rendered link. Refer to the HTML canvas setLineDash API for example values. Either a falsy value or an empty array will disable dashing.
- `linkOpacity` (Real; optional): Line opacity of links, between [0,1]. 3D, VR, AR
- `linkResolution` (Real; optional): Geometric resolution of each link 3D line, expressed in how many radial segments to divide the cylinder. Higher values yield smoother cylinders. Applicable only to links with positive width. 3D, VR, AR
- `linkSource` (String; optional): Link object accessor attribute referring to id of source node.
- `linkTarget` (String; optional): Link object accessor attribute referring to id of target node.
- `linkWidth` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the link line width.
- `maxDepth_neighbours_select` (Real; optional): max levels of neighbourhood selection around a node by repeat clicking
- `ngraphPhysics` (Dict; optional): Specify custom physics configuration for ngraph, according to its configuration object syntax. Only applicable if using the ngraph simulation engine.
- `nodeAutoColorBy` (String; optional): Node object accessor function or attribute to automatically group colors by. Only affects nodes without a color attribute.
- `nodeColor` (String; optional): The node attribute whose value should be used for coloring nodes
- `nodeDesc` (String; optional): For VR only. Node object accessor function or attribute for description (shown under label).
- `nodeIcon` (String; optional): The node attribute containing object with icon to display for each individual node.
- `nodeIcon_fontsheets` (Dict; optional): object with keys being fonts (string) and values being CSS sheets
- `nodeId` (String; optional): Node object accessor attribute for unique node id (used in link objects source/target).
- `nodeIdsHighlight` (Array of Strings; optional): ids of highlighted nodes (through searcha)
- `nodeIdsVisible` (Array of Strings; optional): ids of visible nodes
- `nodeImg` (String; optional): The node attribute containing url to image to display for each individual node
- `nodeLabel` (String; optional): Node object accessor function or attribute for name (shown in label). Supports plain text or HTML content (except in VR).
2D, 3D and VR
- `nodeOpacity` (Real; optional): Nodes sphere opacity, between [0,1]. 3D, VR, AR
- `nodeRelSize` (Real; optional): Ratio of node circle area (square px) [2D] or sphere volume (cubic px) [3D] per value unit.
- `nodeResolution` (Real; optional): Geometric resolution of each node's sphere, expressed in how many slice segments to divide the circumference. Higher values yield smoother spheres. Only applicable to 3D modes.
3D, VR, AR
- `nodeURL` (String; optional): The node attribute containing a URL
- `nodeVal` (Real | String; optional): Ratio of node circle area (square px) [2D] or sphere volume (cubic px) [3D] per value unit.
- `nodeZoomId` (String; optional): id of node to zoom to
- `nodesSelected` (Array of Dicts; optional): selected (clicked) nodes
- `numDimensions` (Real; optional): Not applicable to 2D mode. Number of dimensions to run the force simulation on. 3D, VR, AR
- `pauseAnimation` (Bool; optional): Pauses the rendering cycle of the component, effectively freezing the current view and cancelling all user interaction. This method can be used to save performance in circumstances when a static image is sufficient.
- `pixelUnitRatio` (Real; optional): pixelUnitRatio: if node attribute (in some unit of measurement) is used as coordinates, pixel:unit scale
- `refresh` (Bool; optional): Redraws all the nodes/links. 3D, VR, AR
- `rendererConfig` (Dict; optional): Configuration parameters to pass to the ThreeJS WebGLRenderer constructor. This prop only has an effect on component mount. 3D only
- `resumeAnimation` (Bool; optional): Resumes the rendering cycle of the component, and re-enables the user interaction. This method can be used together with pauseAnimation for performance optimization purposes.
- `showCoordinates` (Bool; optional): showCoordinates: whether or not to show pointer coordinates as tooltip (not yet used)
- `showNavInfo` (Bool; optional): Whether to show the navigation controls footer info.
- `size` (Dict; optional): provided react-sizeme. Contains an object with "width" and "height" attributes
- `sortRels1Descend` (Bool; optional): sort in descending order?
- `sortRels2Descend` (Bool; optional): sort in descending order?
- `sortRelsBy1` (String; optional): in zoom view, node attribute to sort relations by first
- `sortRelsBy2` (String; optional): in zoom view, node attribute to sort relations by after first sort
- `sortRoleplayers1Descend` (Bool; optional): sort in descending order?
- `sortRoleplayers2Descend` (Bool; optional): sort in descending order?
- `sortRoleplayersBy1` (String; optional): in zoom view, node attribute to sort role players by first
- `sortRoleplayersBy2` (String; optional): in zoom view, node attribute to sort role players by after first sort
- `useCoordinates` (Bool; optional): useCoordinates: whether to use node attribute to set node coordinates
- `useNodeIcon` (Bool; optional): Whether or not to use the nodeIcon
- `useNodeImg` (Bool; optional): Whether or not to use the nodeImg. Overrides nodeIcon
- `warmupTicks` (Real; optional): Number of layout engine cycles to dry-run at ignition before starting to render.
- `zoom` (Array; optional): Set the 2D canvas zoom amount. The zoom is defined in terms of the scale transform of each px. A value of 1 indicates unity, larger values zoom in and smaller values zoom out. An optional 2nd argument defines the duration of the transition (in ms) to animate the canvas motion. By default the zoom is set to a value inversely proportional to the amount of nodes in the system.
- `zoomToFit` (Array; optional): Automatically zooms/pans the canvas so that all of the nodes fit inside it. If no nodes are found no action is taken. It accepts two optional arguments: the first defines the duration of the transition (in ms) to animate the canvas motion (default: 0ms). The second argument is the amount of padding (in px) between the edge of the canvas and the outermost node (default: 10px). The third argument specifies a custom node filter: node => <boolean>, which should return a truthy value if the node is to be included. This can be useful for focusing on a portion of the graph. 2D, 3D
"""
function dash_graph3d(; kwargs...)
        available_props = Symbol[:id, :active, :backgroundColor, :cameraPosition, :centerAt, :centreCoordinates, :controlType, :cooldownTicks, :cooldownTime, :d3AlphaDecay, :d3AlphaMin, :d3ReheatSimulation, :d3VelocityDecay, :dagLevelDistance, :dagMode, :dagModeOn, :dagNodeIds, :emitParticle, :enableNodeDrag, :externalobject_input, :externalobject_source, :forceEngine, :getGraphBbox, :graphData, :gravity, :heightRatio, :interactive, :linkAutoColorBy, :linkColor, :linkCurvature, :linkCurveRotation, :linkDesc, :linkDirectionalArrowColor, :linkDirectionalArrowLength, :linkDirectionalArrowRelPos, :linkDirectionalArrowResolution, :linkDirectionalParticleColor, :linkDirectionalParticleResolution, :linkDirectionalParticleSpeed, :linkDirectionalParticleWidth, :linkDirectionalParticles, :linkHoverPrecision, :linkId, :linkIdsHighlight, :linkIdsVisible, :linkLabel, :linkLineDash, :linkOpacity, :linkResolution, :linkSource, :linkTarget, :linkWidth, :maxDepth_neighbours_select, :ngraphPhysics, :nodeAutoColorBy, :nodeColor, :nodeDesc, :nodeIcon, :nodeIcon_fontsheets, :nodeId, :nodeIdsHighlight, :nodeIdsVisible, :nodeImg, :nodeLabel, :nodeOpacity, :nodeRelSize, :nodeResolution, :nodeURL, :nodeVal, :nodeZoomId, :nodesSelected, :numDimensions, :pauseAnimation, :pixelUnitRatio, :refresh, :rendererConfig, :resumeAnimation, :showCoordinates, :showNavInfo, :size, :sortRels1Descend, :sortRels2Descend, :sortRelsBy1, :sortRelsBy2, :sortRoleplayers1Descend, :sortRoleplayers2Descend, :sortRoleplayersBy1, :sortRoleplayersBy2, :useCoordinates, :useNodeIcon, :useNodeImg, :warmupTicks, :zoom, :zoomToFit]
        wild_props = Symbol[]
        return Component("dash_graph3d", "Graph3D", "dash_react_force_graph", available_props, wild_props; kwargs...)
end

max_stars_repo_name: nodelab-aps/dash-react-force-graph

max_issues_repo_name: nodelab-aps/dash-react-force-graph

max_forks_repo_name: nodelab-aps/dash-react-force-graph


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Custom cells for the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    append_suffix_to_line!(c::CustomTextCell, l::Int, suffix::String)

Append the `suffix` to the line `l` of the custom cell text `c`.
"""
function append_suffix_to_line!(c::CustomTextCell, l::Int, suffix::String)
    error("The method `append_suffix_to_line!` is not defined for `$(typeof(c))`.")
end

"""
    apply_line_padding!(c::CustomTextCell, l::Int, left_pad::Int, right_pad::Int)

Apply to the line `l` of the custom text cell `c` the padding with `left_pad`
spaces in the left and `right_pad` spaces in the right.
"""
function apply_line_padding!(c::CustomTextCell, l::Int, left_pad::Int, right_pad::Int)
    error("The method `apply_line_padding` is not defined for `$(typeof(c))`.")
end

"""
    crop_line!(c::CustomTextCell, l::Int, num::Int)

Crop `num` characters from the line `l` of the custom text cell `c`. The number
of cropped characters must consider the left and right alignment paddings.
"""
function crop_line!(c::CustomTextCell, l::Int, num::Int)
    error("The method `crop_line!` is not defined for `$(typeof(c))`.")
end

"""
    get_printable_cell_line(c::CustomTextCell, l::Int)

Return the printable line `l` of the custom text cell `c`. The printable cell
line must consider the left and right alignment paddings.
"""
function get_printable_cell_line(c::CustomTextCell, l::Int)
    error("The method `get_printable_cell_line` is not defined for `$(typeof(c))`.")
end

"""
    get_rendered_line(c::CustomTextCell, l::Int)

Return the rendered line `l` of the custom text cell `l`.
"""
function get_rendered_line(c::CustomTextCell, l::Int)
    error("The method `get_rendered_line` is not defined for `$(typeof(c))`.")
end

"""
    parse_cell_text(c::CustomTextCell; kwargs...)

Parse the cell text and return a vector of `String` with the printable cell
text, where each element in the vector is a new line.

The returned data must contain **only** the printable characters.

The following keyword arguments are passed to this function, which is called
during the cell parsing phase. Those options are related to the input
configuration of `pretty_table`, and the user must choose whether or not support
them.

- `autowrap::Bool`: If `true`, the user wants to wrap the text in the cell. In
    this case, the option `column_width` contains the column width so that the
    text can be wrapped into multiple lines.
- `cell_first_line_only::Bool`: If `true`, the user only wants the first line.
- `column_width::Integer`: The column width.
- `compact_printing::Bool`: If `true`, the user wants compact printing (see
    `:compact` options of `IOContext`).
- `limit_printing::Bool`: If `true`, the user wants the cells to be converted
    using the option `:limit => true` in `IOContext`.
- `linebreaks::Bool`: If `true`, the user wants line breaks inside the cells.
- `renderer::Union{Val{:print}, Val{:show}}`: The render that the user wants to
    convert the cells to strings.
"""
function parse_cell_text(c::CustomTextCell; kwargs...)
    error("The method `parse_cell_text` is not defined for `$(typeof(c))`.")
end

"""
    reset!(c::CustomTextCell)

Reset all fields in the custom text cell `c`.

This function is not required for the API. It is called before parsing the
custom text cell.
"""
reset!(c::CustomTextCell) = nothing


text:
@testset "678.valid-parenthesis-string.jl" begin
    @test check_valid_string("(())((())()()(*)(*()(())())())()()((()())((()))(*") == false
    @test check_valid_string(")(") == false
    @test check_valid_string("(*))") == true
    @test check_valid_string("(*)") == true
    @test check_valid_string("()") == true
end


text:
module Advection

export
    div_ũu, div_ũv, div_ũw, div_uc,

    momentum_flux_uu,
    momentum_flux_uv,
    momentum_flux_uw,
    momentum_flux_vu,
    momentum_flux_vv,
    momentum_flux_vw,
    momentum_flux_wu,
    momentum_flux_wv,
    momentum_flux_ww,
    advective_tracer_flux_x,
    advective_tracer_flux_y,
    advective_tracer_flux_z,

    CenteredSecondOrder,
    UpwindBiasedThirdOrder,
    CenteredFourthOrder,
    WENO, WENO5

using Oceananigans.Grids
using Oceananigans.Operators

abstract type AbstractAdvectionScheme end
abstract type AbstractCenteredAdvectionScheme <: AbstractAdvectionScheme end
abstract type AbstractUpwindBiasedAdvectionScheme <: AbstractAdvectionScheme end

include("topologically_conditional_interpolation.jl")

include("centered_advective_fluxes.jl")
include("upwind_biased_advective_fluxes.jl")

include("centered_second_order.jl")
include("upwind_biased_third_order.jl")
include("centered_fourth_order.jl")

include("weno_reconstruction.jl")
include("weno_nth_order.jl")
include("weno_fifth_order.jl")

include("momentum_advection_operators.jl")
include("tracer_advection_operators.jl")

end # module


text:
module SQLLoader

using Octo: Repo, AdapterBase, DBMS, SQLElement, Structured
using .Repo: SQLKeyword, ExecuteResult

# db_dbname
function db_dbname(nt::NamedTuple)::String
    ""
end

# db_connect
function db_connect(; kwargs...)
end

# db_disconnect
function db_disconnect(conn)
end

# query
function query(conn, sql::String)
end

function query(conn, prepared::String, vals::Vector)
end

# execute
function execute(conn, sql::String)::ExecuteResult
    nothing
end

function execute(conn, prepared::String, vals::Vector)::ExecuteResult
    nothing
end

function execute(conn, prepared::String, nts::Vector{<:NamedTuple})::ExecuteResult
    nothing
end

# execute_result
function execute_result(conn, command::SQLKeyword)::NamedTuple
    NamedTuple()
end

function Base.show(io::IO, mime::MIME"text/plain", element::Union{E,Structured} where E<:SQLElement)
    dbms = DBMS.SQL()
    AdapterBase._show(io, mime, dbms, element)
end

end # module Octo.Backends.SQLLoader


text:
print("\tmodel.jl")

m = 8
n = 8
A = sprand(m, n, 0.5)
c = 1.0 - 2.0*rand(n)
b = 1.0 + rand(m)
u = 10.0 + sprand(n, 0.8)
p = nnz(u)

# Tests for PrimalDualPoint
v1 = Tulip.PrimalDualPoint(rand(n), rand(p), rand(m), rand(n), rand(p))
v2 = copy(v1)
@test n == sum(v1.x .== v2.x)
@test p == sum(v1.w .== v2.w)
@test m == sum(v1.y .== v2.y)
@test n == sum(v1.s .== v2.s)
@test p == sum(v1.z .== v2.z)

v3 = v1 + v2
@test n == sum(v3.x .== (v1.x + v2.x))
@test p == sum(v3.w .== (v1.w + v2.w))
@test m == sum(v3.y .== (v1.y + v2.y))
@test n == sum(v3.s .== (v1.s + v2.s))
@test p == sum(v3.z .== (v1.z + v2.z))

# Create random instance and run checks
model = Tulip.Model(A, b, c, u.nzind, u.nzval)

@test m == model.nconstr
@test n == model.nvars
@test model.status == :Built

sol = model.sol
@test typeof(sol) <: Tulip.PrimalDualPoint
@test n == size(sol.x, 1)
@test p == size(sol.w, 1)
@test m == size(sol.y, 1)
@test n == size(sol.s, 1)
@test p == size(sol.z, 1)

println("\tPassed.")


text:
module jcontrol
# using Term
# install_term_logger()

using Plots
using Interpolations

include("utils.jl")
include("kinematics.jl")
include("track.jl")
include("bike.jl")
include("control.jl")
include("forward_model.jl")
include("comparisons.jl")
include("run.jl")
include("trial.jl")
include("visuals.jl")
include("io.jl")

export Track, get_track_borders, FULLTRACK
export State, Bicycle
export ControlOptions, create_and_solve_control, Bounds
export Solution, run_forward_model
export summary_plot
export PATHS, load_trials, load_cached_trials
export arena
export DynamicsProblem, KinematicsProblem, realistict_control_options
export get_comparison_points
export run_mtm
export Trial, trimtrial, get_varvalue_at_frameidx

using .bicycle: State, Bicycle
using .control:
    ControlOptions,
    create_and_solve_control,
    Bounds,
    State,
    DynamicsProblem,
    KinematicsProblem,
    realistict_control_options

using .forwardmodel: Solution, run_forward_model
using .comparisons: ComparisonPoints, ComparisonPoint, track_segments, TrackSegment
using .visuals
using .Run: run_mtm
using .trial: Trial, trimtrial, get_varvalue_at_frameidx
using .io: PATHS, load_trials, load_cached_trials

end # module

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
module FaceDatasets

export load_images,
       load_shapes,
       load_labels,
       FaceDataset,
       CootesDataset,
       CKDataset,
       CKMaxDataset,
       PutDataset,
       PutFrontalDataset

include("core.jl")

end


text:
struct SpirvCrossError <: Exception
    msg::AbstractString
    return_code
end

Base.showerror(io::IO, e::SpirvCrossError) = print(io, "$(e.return_code): ", e.msg)

"""
    @check spvc_create_context(ctx)

Checks whether the expression returned SPVC_SUCCESS. Else, throw an error printing the corresponding code.
"""
macro check(expr)
    quote
        local msg = "failed to execute " * $(string(expr))
        @check $(esc(expr)) msg
    end
end

macro check(expr, msg)
    quote
        return_code = $(esc(expr))
        if Int(return_code) < 0
            throw(SpirvCrossError($msg, return_code))
        end
        return_code
    end
end

"""
Wraper around an opaque handle that refers to internal SPVC data.
"""
abstract type Handle end

Base.cconvert(T::Type{<:Ptr}, x::AbstractArray{<:Handle}) = getproperty.(x, :handle)
Base.cconvert(T::Type{<:Ptr}, x::Handle) = x

Base.unsafe_convert(T::Type{Ptr{Nothing}}, x::Handle) = x.handle

Base.broadcastable(x::Handle) = Ref(x)

null_ptr_ref() = Ref{Ptr{Cvoid}}()


text:
using XConv, CUDA, LinearAlgebra, NNlib, Flux
import Flux: mse
using BenchmarkTools
# CUDA param
CUDA.allowscalar(false)
# BLAS param
nthreads = div(Sys.CPU_THREADS, 2)
BLAS.set_num_threads(nthreads)
XConv.initXConv(0, "TrueGrad")

nx = 512
ny = 512
b = 32
n_in = 16
n_out = 16
stride = 1
n_bench = 5
nw   = 3;

@info "Running runntime comparison between CPU and CPU for N=$(nx)x$(ny), batchsize=$b, $(n_in) input channels and $(n_out) output channels"

# Flux network
C = Conv((nw, nw), n_in=>n_out, identity;pad=1, stride=stride)
X = randn(Float32, nx, ny, n_in, b)
Y = C(X)
cdims = DenseConvDims(X, C.weight; padding=1)

@btime gcpu = NNlib.∇conv_filter($X, $Y, $cdims)
@btime gcpu_ev = grad_ev($X, $Y, 8, $nw, $stride);

@info "CPU benchmark with $(nthreads) threads"

for i=1:4
  @info "test $i"
  @time gcpu = NNlib.∇conv_filter(X, Y, cdims)
  @time gcpu_ev = grad_ev(X, Y, 8, nw, stride);

  XConv.initXConv(0, "TrueGrad")
  @time gradient(C->.5f0*norm(C(X)), C)
  XConv.initXConv(2^3, "EVGrad")
  @time gradient(C->.5f0*norm(C(X)), C)
  XConv.initXConv(0, "TrueGrad")
end

XConv.initXConv(0, "TrueGrad")
@btime gradient(C->.5f0*norm(C(X)), C)
XConv.initXConv(2^3, "EVGrad")
@btime gradient(C->.5f0*norm(C(X)), C)
XConv.initXConv(0, "TrueGrad")

###### CUDA

for dtype=[Float32]
    @info "GPU benchmark in $dtype"

    C = Conv(CUDA.randn(dtype, nw, nw, n_in, n_out), CUDA.zeros(dtype, n_out), identity;pad=1)
    X = CUDA.randn(dtype, nx, ny, n_in, b)
    Y = C(X)
    w = C.weight

    cdims = DenseConvDims(X, w; padding=1)

    XConv.initXConv(0, "TrueGrad")
    @btime ggpu = NNlib.∇conv_filter($X, $Y, $cdims)
    ggpu_ev = grad_ev(X, Y, 8, nw, stride)
    @btime ggpu_ev = grad_ev($X, $Y, 8, $nw, $stride);
    @btime ggpu_ev = grad_ev($X, $Y, 8, $nw, $stride);

    for i=1:4
        @info "test $i"
        XConv.initXConv(0, "TrueGrad")
        @time ggpu = NNlib.∇conv_filter(X, Y, cdims)
        @time ggpu_ev = grad_ev(X, Y, 8, nw, stride);

        XConv.initXConv(0, "TrueGrad")
        @time gradient(C->.5f0*norm(C(X)), C)
        XConv.initXConv(2^3, "EVGrad")
        @time gradient(C->.5f0*norm(C(X)), C)
        XConv.initXConv(0, "TrueGrad")
    end

    XConv.initXConv(0, "TrueGrad")
    @time gradient(C->.5f0*norm(C(X))^2, C)
    XConv.initXConv(2^3, "EVGrad")
    @time gradient(C->.5f0*norm(C(X))^2, C)
    XConv.initXConv(0, "TrueGrad")
    @btime gradient(C->.5f0*norm(C($X))^2, $C)
    XConv.initXConv(2^3, "EVGrad")
    @btime gradient(C->.5f0*norm(C($X))^2, $C)

end

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv


text:
module CanonicalTraits
using MLStyle

export @trait, @implement, instance, @implement!

include("Utils.jl")
include("Typeclasses.jl")
include("Instances.jl")

end # module


text:
using ModuleDocstrings
using Documenter

DocMeta.setdocmeta!(ModuleDocstrings, :DocTestSetup, :(using ModuleDocstrings); recursive=true)

makedocs(;
    modules=[ModuleDocstrings],
    authors="Tim Holy <tim.holy@gmail.com> and contributors",
    repo="https://github.com/JuliaDocs/ModuleDocstrings.jl/blob/{commit}{path}#{line}",
    sitename="ModuleDocstrings.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaDocs.github.io/ModuleDocstrings.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaDocs/ModuleDocstrings.jl",
    devbranch="main",
    push_preview = true,
)

max_stars_repo_name: JuliaDocs/ModuleDocstrings

max_issues_repo_name: JuliaDocs/ModuleDocstrings

max_forks_repo_name: JuliaDocs/ModuleDocstrings


text:
__precompile__()
module REPLTetris

using Crayons, REPLGamesBase
export tetris

include("tiles.jl")
include("board.jl")
include("printboard.jl")
include("actions.jl")

function tetris(board = Board())
    rawmode() do
        clear_screen()
        update_board!(board)
        abort = [false]
        @async while !abort[1] && add_tile!(board)
            board.allowhold = true
            print_preview(board)
            print_hold(board)
            while !abort[1] && drop!(board)
                sleep((0.8 - (board.level-1) * 0.007)^(board.level-1))
            end
            delete_lines!(board)
        end

        while !abort[1]
            c = readKey()
            c in ["Up", "x"]    && rot_right!(board)
            c in ["Down"]       && drop!(board)
            c in ["Right"]      && move_right!(board)
            c in ["Left"]       && move_left!(board)
            c in [" "]          && fast_drop!(board)
            c in ["Ctrl", "z"]  && rot_left!(board)
            c in ["c"]          && hold!(board)
            c in ["Ctrl-C"]     && (abort[1]=true)
        end
    end
end

end #module

max_stars_repo_name: ChristianKurz/REPLTetris

max_issues_repo_name: ChristianKurz/REPLTetris

max_forks_repo_name: ChristianKurz/REPLTetris


text:
export TestSystem, nvariables, equations, mixed_volume, nsolutions,
    nreal_solutions, bezout_number, multi_bezout_number

struct TestSystem{T}
    equations::Vector{Polynomial{true, T}}
    bezout_number::Int
    multi_homogeneous_bezout_number::Union{Nothing, Tuple{Int, Vector{Vector{PolyVar{true}}}}}
    mixed_volume::Union{Nothing, Int}
    nsolutions::Union{Nothing, Int}
    nreal_solutions::Union{Nothing, Int}
end

"""
    TestSystem(equations;
        multi_bezout_number=nothing,
        mixed_volume=nothing,
        nsolutions=nothing,
        nreal_solutions=nothing)

Create a `TestSystem`.
"""
function TestSystem(equations;
    multi_bezout_number=nothing,
    mixed_volume=nothing,
    nsolutions=nothing,
    nreal_solutions=nothing)
    bezout_number = prod(MP.maxdegree, equations)
    TestSystem(equations, bezout_number, multi_bezout_number, mixed_volume,nsolutions, nreal_solutions)
end

"""
    equations(::TestSystem)

Get the the polynomial system.
"""
equations(S::TestSystem) = S.equations

"""
    nvariables(::TestSystem)

Obtain the number of variables of the test system.
"""
nvariables(S::TestSystem) = maximum(MP.nvariables, S.equations)

"""
    mixed_volume(system)::Union{Nothing, Int}

Returns the number mixed volume of the system if known.
"""
mixed_volume(F::TestSystem) = F.mixed_volume

"""
    nsolutions(system)::Union{Nothing, Int}

Returns the number of (complex) solutions of the system if known.
"""
nsolutions(F::TestSystem) = F.nsolutions

"""
    nreal_solutions(system)::Union{Nothing, Int}

Returns the number of real solutions of the system if known.
"""
nreal_solutions(F::TestSystem) = F.nreal_solutions

"""
    bezout_number(system)

Returns the bezout number of the system.
"""
bezout_number(F::TestSystem) = F.bezout_number

"""
    multi_bezout_number(system)::Union{Nothing, Tuple{Int, Vector{Vector{PolyVar{true}}}}}

Returns a tuple containing the multi-homogeneous bezout number
as well as the corresponding grouping of the variables if known.
"""
multi_bezout_number(F::TestSystem) = F.multi_homogeneous_bezout_number


text:
module NCBITaxonomy
using DataFrames
using Arrow
using StringDistances

if !haskey(ENV, "NCBITAXONOMY_PATH")
    @warn """
    The environmental varialbe NCBITAXONOMY_PATH is not set, so the tables will
    be stored in the package path. This is not ideal, and you really should set
    the NCBITAXONOMY_PATH.
    """
end
const taxpath = get(ENV, "NCBITAXONOMY_PATH", joinpath(@__DIR__, "..", "deps"))
ispath(taxpath) || mkpath(taxpath)

function __init__()
    name_date = mtime(joinpath(taxpath, "tables", "names.arrow"))
    return time() - name_date >= 2.6e+6 && @warn(
        "Your local taxonomy version is over 30 days old, we recommend using `] build NCBITaxonomy` to get the most recent version."
    )
end

include("types.jl")
export NCBITaxon, NCBINameClass, IDNotFoundInBackbone

include("exceptions.jl")
export NameHasNoDirectMatch, NameHasMultipleMatches

names_table = DataFrame(Arrow.Table(joinpath(taxpath, "tables", "names.arrow")))
names_table.class = NCBINameClass.(names_table.class)
names_table.lowercase = lowercase.(names_table.name)

division_table = DataFrame(Arrow.Table(joinpath(taxpath, "tables", "division.arrow")))
select!(division_table, Not(:comments))

nodes_table = DataFrame(Arrow.Table(joinpath(taxpath, "tables", "nodes.arrow")))
select!(nodes_table, Not(r"inherited_"))
select!(nodes_table, Not(r"_code_id"))
select!(nodes_table, Not(:genbank_hidden))
select!(nodes_table, Not(:hidden_subtree))
select!(nodes_table, Not(:comments))
select!(nodes_table, Not(:embl))

nodes_table = innerjoin(nodes_table, division_table; on=:division_id)
select!(nodes_table, Not(:division_id))

names_table = leftjoin(
    names_table, unique(select(nodes_table, [:tax_id, :rank])); on=:tax_id
)

include("taxon.jl")
export taxon, @ncbi_str

include("namefilters/namefilter.jl")
include("namefilters/divisions.jl")
export namefilter
export bacteriafilter,
    virusfilter,
    mammalfilter,
    vertebratefilter,
    plantfilter,
    invertebratefilter,
    rodentfilter,
    primatefilter,
    environmentalsamplesfilter,
    phagefilter

include("lineage/children.jl")
include("lineage/lineage.jl")
export children, descendants, lineage, parent, rank

include("utility/nametools.jl")
include("utility/similarnames.jl")
include("utility/taxonomicdistance.jl")
export vernacular, synonyms, authority, taxonomicdistance, taxonomicdistance!
export similarnames, alternativetaxa

end


text:
@testset "Logical" begin
    l = Logical(1, (β, xmat) -> xmat * β, false)
    @test isa(l, Logical)
    @test length(l.monitor) == 0
    l = Logical((β, xmat) -> xmat * β, false)
    @test isa(l, Logical)

    l = Logical((β, xmat) -> xmat * β, 1,  false)
    @test isa(l, Logical)

    s = Logical(Node(), (x)->x, true)
    @test isa(s, Logical)
end

@testset "Stochastic" begin
    s = Stochastic(1, () -> Normal(0, sqrt(1000)), true)
    @test isa(s, Stochastic)
    @test length(s.monitor) == 1
    s = Stochastic(() -> Normal(0, sqrt(1000)), true)
    @test isa(s, Stochastic)

    s = Stochastic(() -> Normal(0, sqrt(1000)),1, true)
    @test isa(s, Stochastic)

    s = Stochastic(Node(), ()->UniformBranchLength(), true)
    @test isa(s, Stochastic)
end

@testset "Model" begin
    model = Model(
    y = Stochastic(1, (μ, s2) -> MultivariateNormal(μ, s2), false),
    μ = Logical(1, (β, xmat) -> xmat * β, false),
    β = Stochastic(1, () -> Normal(0, sqrt(1000)), true),
    s2 = Stochastic(() -> InverseGamma(0.001, 0.001), true),
        )

    @test isa(model, Model)
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
#TODO add unit tests

@testset "test opf with variable impedance" begin
    file =  "../test/data/matpower/case5.m"
    data = PMs.parse_file(file)
    # result = PMPP.run_opf_variable_impedance_loss(data, ipopt_solver)
    # @test result["termination_status"] == PMs.LOCALLY_SOLVED
    # @test isapprox(result["objective"], 0; atol = 1e0)
end


text:
module OceanBoundaryLayerParameterizations

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using OceanTurbulenceParameterEstimation
using OceanTurbulenceParameterEstimation.Models: set!
using CairoMakie, LaTeXStrings, OrderedCollections
using Oceananigans, FileIO

export
    # lesbrary_paths.jl
    TwoDaySuite, FourDaySuite, SixDaySuite, GeneralStrat,

    # one_dimensional_ensemble_model.jl
    OneDimensionalEnsembleModel,

    # EKI_hyperparameter_search.jl
    plot_stds_within_bounds, plot_prior_variance, plot_num_ensemble_members,
    plot_observation_noise_level, plot_prior_variance_and_obs_noise_level,

    # visualize_profile_predictions.jl
    visualize!,
    visualize_and_save!

include("lesbrary_paths.jl")
include("one_dimensional_ensemble_model.jl")
include("EKI_hyperparameter_search.jl")
include("visualize_profile_predictions_utils.jl")
include("visualize_profile_predictions.jl")

end # module


text:
struct IsRadonDiag end
struct IsRadonSquare end

include("radon_diag.jl")
include("radon_square.jl")

@inline function _radon(
    ::IsRadonDiag,
    a::A,
    image::AbstractMatrix{T};
    nd::Optional{I} = nothing,
    nϕ::Optional{J} = nothing,
    α::Real = 360,
    α₀::Real = 0,
    kwargs...
) where {
    A <: AbstractProjectionAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    rows, cols = size(image)
    h::T = hypot(rows, cols)
    nd = isnothing(nd) ? round(Int, h) : nd
    nϕ = isnothing(nϕ) ? 2 * ((rows + 1) * (cols + 1) ÷ (2nd)) + 1 : nϕ
    ts = linspace(T, -1..1, nd)
    ϕ₀::T = deg2rad(α₀)
    ϕ₁::T = ϕ₀ + deg2rad(α)
    ϕs = linspace(ORI(ϕ₀..ϕ₁), nϕ)
    a(image, ts, ϕs; kwargs...)
end

@inline function _radon(
    ::IsRadonDiag,
    a::A,
    image::AbstractMatrix{T},
    tsi::Interval,
    ϕsi::Interval;
    nd::Optional{I} = nothing,
    nϕ::Optional{J} = nothing,
    kwargs...
) where {
    A <: AbstractProjectionAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    rows, cols = size(image)
    h::T = hypot(rows, cols)
    nd = isnothing(nd) ? round(Int, h) : nd
    nϕ = isnothing(nϕ) ? 2 * ((rows + 1) * (cols + 1) ÷ (2nd)) + 1 : nϕ
    ts = linspace(T, tsi, nd)
    ϕs = linspace(ϕsi, nϕ)
    a(image, ts, ϕs; kwargs...)
end

@inline function _radon(
    ::IsRadonSquare,
    a::A,
    image::AbstractMatrix{T};
    nd::Optional{I} = nothing,
    nϕ::Optional{J} = nothing,
    α::Real = 360,
    α₀::Real = 0,
    kwargs...
) where {
    A <: AbstractProjectionAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    rows, cols = size(image)
    l = min(rows, cols)
    nd = maybe(l, nd)
    nϕ = isnothing(nϕ) ? 2 * ((rows + 1) * (cols + 1) ÷ (2nd)) + 1 : nϕ
    ts = linspace(-1..1, nd)
    ϕ₀::T = deg2rad(α₀)
    ϕ₁::T = ϕ₀ + deg2rad(α)
    ϕs = linspace(ORI(ϕ₀..ϕ₁), nϕ)
    a(image, ts, ϕs; kwargs...)
end

@inline function _radon(
    ::IsRadonSquare,
    a::A,
    image::AbstractMatrix{T},
    tsi::Interval,
    ϕsi::Interval;
    nd::Optional{I} = nothing,
    nϕ::Optional{J} = nothing,
    kwargs...
) where {
    A <: AbstractProjectionAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    rows, cols = size(image)
    l = min(rows, cols)
    nd = maybe(l, nd)
    nϕ = isnothing(nϕ) ? 2 * ((rows + 1) * (cols + 1) ÷ (2nd)) + 1 : nϕ
    ts = linspace(tsi, nd)
    ϕs = linspace(ϕsi, nϕ)
    a(image, ts, ϕs; kwargs...)
end

"""
    radon(image::AbstractMatrix[, algorithm::AbstractProjectionAlgorithm]; <keyword arguments>)

Compute the Radon transform of `image` with parameters given as keyword
arguments. The parameters depend on the algorithm, please see the relative
documentation. The default algorithm is [`Radon`](@ref).

See Also: [`project_image`](@ref)
"""
@inline function radon(
    image::AbstractMatrix,
    alg::AbstractProjectionAlgorithm = Radon();
    kwargs...
)
    alg(image; kwargs...)
end

"""
    radon(image::AbstractMatrix, ts, ϕs[, algorithm::AbstractProjectionAlgorithm]; <keyword arguments>)

Compute the Radon transform of `image` with integration points given by the
vector `ts` and projection angles given by the vector `ϕs`. Additional
parameters can be passed to the algorithm through keyword arguments.Please see
the relative documentation. The default algorithm is [`Radon`](@ref).

See Also: [`project_image`](@ref)
"""
@inline function radon(
    image::AbstractMatrix,
    ts::AbstractVector,
    ϕs::AbstractVector,
    alg::AbstractProjectionAlgorithm = Radon();
    kwargs...
)
    alg(image, ts, ϕs; kwargs...)
end

@inline function radon(
    image::AbstractMatrix,
    ts::Interval,
    ϕs::Interval,
    alg::AbstractProjectionAlgorithm = Radon();
    kwargs...
)
    _radon(alg, image, ts, ϕs; kwargs...)
end

"""
    radon(image::AbstractMatrix, geometry::AbstractGeometry[, algorithm::AbstractProjectionAlgorithm]; <keyword arguments>)

Compute the Radon transform of `image` with explicit geometry. Additional
parameters can be passed to the algorithm through keyword arguments. Please see
the relative documentation. The default algorithm is [`Radon`](@ref).

See Also: [`project_image`](@ref)
"""
@inline function radon(
    image::AbstractMatrix,
    geometry::AbstractParallelBeamGeometry,
    alg::AbstractProjectionAlgorithm = Radon();
    kwargs...
)
    nd = geometry.nd
    nϕ = geometry.nϕ
    α = geometry.α
    α₀ = geometry.α₀
    radon(image, alg; nd, nϕ, α, α₀, kwargs...)
end

@inline function radon(
    image::AbstractMatrix,
    geometry::AbstractFanBeamGeometry,
    alg::AbstractProjectionAlgorithm = Radon();
    kwargs...
)
    g′ = ParallelBeamGeometry(geometry)
    sinog = radon(image, g′, alg; kwargs...)
    _, fan_sinog = para2fan(sinog, geometry)
    fan_sinog
end

@inline function radon(
    g::AbstractGeometry,
    alg::AbstractProjectionAlgorithm = Radon();
    kwargs...
)
    x -> radon(x, g, alg; kwargs...)
end