text:
# Maximum distance(difference) for a trigram in a language profile and text profile.
 const MAX_TRIGRAM_DISTANCE = 300;

# 300 trigrams where each has MAX_TOTAL_DISTANCE=300, gives us 90_000.
 const MAX_TOTAL_DISTANCE = 90_000;

# Double MAX_TRIGRAM_DISTANCE
 const TEXT_TRIGRAMS_SIZE = 600;

const RELIABLE_CONFIDENCE_THRESHOLD = 0.8;

function detect_script(text::AbstractString)
    script_counters = [
      [LatinScript()      , 0],
      [CyrillicScript()   , 0],
      [ArabicScript()     , 0],
      [MandarinScript()   , 0],
      [DevanagariScript() , 0],
      [HebrewScript()     , 0],
      [EthiopicScript()   , 0],
      [GeorgianScript()   , 0],
      [BengaliScript()    , 0],
      [HangulScript()     , 0],
      [HiraganaScript()   , 0],
      [KatakanaScript()   , 0],
      [GreekScript()      , 0],
      [KannadaScript()    , 0],
      [TamilScript()      , 0],
      [ThaiScript()       , 0],
      [GujaratiScript()   , 0],
      [GurmukhiScript()   , 0],
      [TeluguScript()     , 0],
      [MalayalamScript()  , 0],
      [OriyaScript()      , 0],
      [MyanmarScript()    , 0],
      [SinhalaScript()    , 0],
      [KhmerScript()      , 0]
    ]

    half = length(text) / 2
    local found
    for ch in text
        if is_stop_char(ch); continue; end
        for i in 1:length(script_counters)
            (script, count) = script_counters[i]
            found = false
            if is_script(script, ch)
                script_counters[i][2] += 1
                found=true
                if script_counters[i][2] > half; return script; end
            end
            # If script was found, move it closer to the front.
            # If the text contains largely 1 or 2 scripts, this will
            # cause these scripts to be eventually checked first.
            if found && i>1
                t=script_counters[i-1]
                script_counters[i-1] = script_counters[i]
                script_counters[i] = t
            end
        end
    end

    sort!(script_counters, lt=(x,y)->x[2]<y[2])
    if script_counters[1][2] > 0
        return script_counters[1][2]
    else
        return nothing
    end

end

function is_script(::CyrillicScript, ch::Char)
   if ch in '\u0400':'\u0484' || ch in '\u0487':'\u052F' || ch in '\u2DE0':'\u2DFF' || ch in '\uA640':'\uA69D' || ch == '\u1D2B' || ch == '\u1D78' || ch == '\uA69F' return true; else; return false; end
end

# https:#en.wikipedia.org/wiki/Latin_script_in_Unicode
function is_script(::LatinScript, ch::Char)
    if ch in 'a':'z' || ch in 'A':'Z' || ch in '\u0080':'\u00FF' || ch in '\u0100':'\u017F' || ch in '\u0180':'\u024F' || ch in '\u0250':'\u02AF' || ch in '\u1D00':'\u1D7F' || ch in '\u1D80':'\u1DBF' || ch in '\u1E00':'\u1EFF' || ch in '\u2100':'\u214F' || ch in '\u2C60':'\u2C7F' || ch in '\uA720':'\uA7FF' || ch in '\uAB30':'\uAB6F' ; return true; else; return false; end
end

# Based on https:#en.wikipedia.org/wiki/Arabic_script_in_Unicode
function is_script(::ArabicScript, ch::Char)
    if ch in  '\u0600':'\u06FF' || ch in '\u0750':'\u07FF' || ch in '\u08A0':'\u08FF' || ch in '\uFB50':'\uFDFF' || ch in '\uFE70':'\uFEFF'  return true; else; return false; end  #|| '\u10E60':'\u10E7F' || '\u1EE00':'\u1EEFF'
end

# Based on https:#en.wikipedia.org/wiki/Devanagari#Unicode
function is_script(::DevanagariScript, ch::Char)
    if ch in '\u0900':'\u097F' || ch in '\uA8E0':'\uA8FF' || ch in '\u1CD0':'\u1CFF' ; return true; else; return false; end
end

# Based on https:#www.key-shortcut.com/en/writing-systems/ethiopian-script/
function is_script(::EthiopicScript, ch::Char)
    if ch in '\u1200':'\u139F' || ch in '\u2D80':'\u2DDF' || ch in '\uAB00':'\uAB2F' ; return true; else; return false; end
end

# Based on https:#en.wikipedia.org/wiki/Hebrew_(Unicode_block)
function is_script(::HebrewScript, ch::Char)
    if ch in '\u0590':'\u05FF' ; return true; else; return false; end
end

function is_script(::GeorgianScript, ch::Char)
   if ch in '\u10A0':'\u10FF' ; return true; else; return false; end
end

function is_script(::MandarinScript, ch::Char)
    if ch in '\u2E80':'\u2E99' || ch in '\u2E9B':'\u2EF3' || ch in '\u2F00':'\u2FD5' || ch == '\u3005' || ch == '\u3007' || ch in '\u3021':'\u3029' || ch in '\u3038':'\u303B' || ch in '\u3400':'\u4DB5' || ch in '\u4E00':'\u9FCC' || ch in '\uF900':'\uFA6D' || ch in '\uFA70':'\uFAD9' ; return true; else; return false; end
end

function is_script(::BengaliScript, ch::Char)
   if ch in '\u0980':'\u09FF' ; return true; else; return false; end
end

function is_script(::HiraganaScript, ch::Char)
   if ch in '\u3040':'\u309F' ; return true; else; return false; end
end

function is_script(::KatakanaScript, ch::Char)
   if ch in '\u30A0':'\u30FF' ; return true; else; return false; end
end

# Hangul is Korean Alphabet. Unicode ranges are taken from: https:#en.wikipedia.org/wiki/Hangul
function is_script(::HangulScript, ch::Char)
    if ch in '\uAC00':'\uD7AF' || ch in '\u1100':'\u11FF' || ch in '\u3130':'\u318F' || ch in '\u3200':'\u32FF' || ch in '\uA960':'\uA97F' || ch in '\uD7B0':'\uD7FF' || ch in '\uFF00':'\uFFEF' ; return true; else; return false; end
end

# Taken from: https:#en.wikipedia.org/wiki/Greek_and_Coptic
function is_script(::GreekScript, ch::Char)
    if ch in '\u0370':'\u03FF' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Kannada_(Unicode_block)
function is_script(::KannadaScript, ch::Char)
    if ch in '\u0C80':'\u0CFF' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Tamil_(Unicode_block)
function is_script(::TamilScript, ch::Char)
    if ch in '\u0B80':'\u0BFF' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Thai_(Unicode_block)
function is_script(::ThaiScript, ch::Char)
    if ch in '\u0E00':'\u0E7F' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Gujarati_(Unicode_block)
function is_script(::GujaratiScript, ch::Char)
    if ch in '\u0A80':'\u0AFF' ; return true; else; return false; end
end

# Gurmukhi is the script for Punjabi language.
# Based on: https:#en.wikipedia.org/wiki/Gurmukhi_(Unicode_block)
function is_script(::GurmukhiScript, ch::Char)
    if ch in '\u0A00':'\u0A7F' ; return true; else; return false; end
end

function is_script(::TeluguScript, ch::Char)
    if ch in '\u0C00':'\u0C7F' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Malayalam_(Unicode_block)
function is_script(::MalayalamScript, ch::Char)
    if ch in '\u0D00':'\u0D7F' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Malayalam_(Unicode_block)
function is_script(::OriyaScript, ch::Char)
    if ch in '\u0B00':'\u0B7F' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Myanmar_(Unicode_block)
function is_script(::MyanmarScript, ch::Char)
    if ch in '\u1000':'\u109F' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Sinhala_(Unicode_block)
function is_script(::SinhalaScript, ch::Char)
    if ch in '\u0D80':'\u0DFF' ; return true; else; return false; end
end

# Based on: https:#en.wikipedia.org/wiki/Khmer_alphabet
function is_script(::KhmerScript, ch::Char)
    if ch in '\u1780':'\u17FF' || ch in '\u19E0':'\u19FF' ; return true; else; return false; end
end

function is_stop_char(ch::Char)
    if ch in '\u0000':'\u0040' || ch in '\u005B':'\u0060' || ch in '\u007B':'\u007E' ; return true; else; return false; end
end

### Trigram model

function to_trigram_char(ch::Char)
    if is_stop_char(ch); ' ' ; else ch; end
end

function count_trigrams(text::AbstractString)
    counter_hash  = Dict{String, Int}()
    text = string(text, ' ')
    #iterate through the string and count trigram
    chars_iter = (lowercase(to_trigram_char(x)) for x in text)
    c1 = ' ';
    #state = start(chars_iter)
    next = iterate(chars_iter)
    (c2, state) = next
    #c2, state = next(chars_iter, state)
    next = iterate(chars_iter, state)
    while next != nothing #!done(chars_iter, state)
        #c3, state  = next(chars_iter, state)
        (c3, state) = next
        if !(c2 == ' ' && (c1 == ' ' || c3 == ' '))
            trigram = string(c1, c2, c3)
            counter_hash[trigram] = get(counter_hash, trigram, 0) + 1
        end
        c1 = c2;
        c2 = c3;
        next = iterate(chars_iter, state)
    end
    return counter_hash
end

function get_trigrams_with_positions(text::String)
    count_vec = [(v,k) for (k, v) in count_trigrams(text)]
    sort!(count_vec, rev=true)
    Dict(tri => i for (i, (c, tri)) in enumerate(count_vec[1:min(length(count_vec), TEXT_TRIGRAMS_SIZE)]))
end

detect_lang_based_on_script(text::AbstractString, script::LatinScript, options) = detect_lang_trigrams(text, script, options)
detect_lang_based_on_script(text::AbstractString, script::CyrillicScript, options) = detect_lang_trigrams(text, script, options)
detect_lang_based_on_script(text::AbstractString, script::DevanagariScript, options) = detect_lang_trigrams(text, script, options)
detect_lang_based_on_script(text::AbstractString, script::HebrewScript, options) = detect_lang_trigrams(text, script, options)
detect_lang_based_on_script(text::AbstractString, script::EthiopicScript, options) = detect_lang_trigrams(text, script, options)
detect_lang_based_on_script(text::AbstractString, script::ArabicScript, options) = detect_lang_trigrams(text, script, options)
detect_lang_based_on_script(text::AbstractString, script::MandarinScript, options) = ("Cmn", 1.0)
detect_lang_based_on_script(text::AbstractString, script::BengaliScript, options) = ("Ben", 1.0)
detect_lang_based_on_script(text::AbstractString, script::HangulScript, options) = ("Kor", 1.0)
detect_lang_based_on_script(text::AbstractString, script::GeorgianScript, options) = ("Kat", 1.0)
detect_lang_based_on_script(text::AbstractString, script::GreekScript, options) = ("Ell", 1.0)
detect_lang_based_on_script(text::AbstractString, script::KannadaScript, options) = ("Kan", 1.0)
detect_lang_based_on_script(text::AbstractString, script::TamilScript, options) = ("Tam", 1.0)
detect_lang_based_on_script(text::AbstractString, script::ThaiScript, options) = ("Tha", 1.0)
detect_lang_based_on_script(text::AbstractString, script::GujaratiScript, options) = ("Guj", 1.0)
detect_lang_based_on_script(text::AbstractString, script::GurmukhiScript, options) = ("Pan", 1.0)
detect_lang_based_on_script(text::AbstractString, script::TeluguScript, options) = ("Tel", 1.0)
detect_lang_based_on_script(text::AbstractString, script::MalayalamScript, options) = ("Mal", 1.0)
detect_lang_based_on_script(text::AbstractString, script::OriyaScript, options) = ("Ori", 1.0)
detect_lang_based_on_script(text::AbstractString, script::MyanmarScript, options) = ("Mya", 1.0)
detect_lang_based_on_script(text::AbstractString, script::SinhalaScript, options) = ("Sin", 1.0)
detect_lang_based_on_script(text::AbstractString, script::KhmerScript, options) = ("Khm", 1.0)
detect_lang_based_on_script(text::AbstractString, script::Union{HiraganaScript, KatakanaScript}, options) = ("Jpn", 1.0)

@enum DetectType Trigram=1 Deep=2

struct LangDetectOptions
    model::DetectType
    whitelist::Union{Nothing, Vector{Language}}
    blacklist::Union{Nothing, Vector{Language}}
end

default_options() = LangDetectOptions(Trigram, Vector{Language}(), Vector{Language}())

detect_lang_trigrams(text::AbstractString, script::Script, options) = detect_lang_trigrams(text, trigram_models[name(script)], options)

function detect_lang_trigrams(text::AbstractString, lang_trigrams_list::Dict{String, Vector{String}}, options )
    lang_distances = Vector{Tuple{String, Int}}()
    trigrams = get_trigrams_with_positions(text)

    for (lang, lang_trigrams) in lang_trigrams_list

        if !isempty(options.whitelist) && !contains(==, options.whitelist, lang); continue; end
        if !isempty(options.blacklist) && contains(==, options.blacklist, lang); continue; end

        dist = calculate_distance(lang_trigrams, trigrams)
        push!(lang_distances, (lang, dist))
    end

    # Sort languages by distance
    sort!(lang_distances, lt=(x, y)->x[2]<y[2])

    # Return None if lang_distances is empty
    if isempty(lang_distances); return (nothing, 0.0); end
    # Return the only language with is_reliable=true if there is only 1 item
    if length(lang_distances) == 1
        return (lang_distances[1][1], 1.0)
    end

    # Calculate is_reliable based on:
    # - number of unique trigrams in the text
    # - rate (diff between score of the first and second languages)
    #
    lang_dist1 = lang_distances[1]
    lang_dist2 = lang_distances[2]
    score1 = MAX_TOTAL_DISTANCE - lang_dist1[2]
    score2 = MAX_TOTAL_DISTANCE - lang_dist2[2]

    if score1 == 0
        # If score1 is 0, score2 is 0 as well, because array is sorted.
        # Therefore there is not language to return.
        return (nothing, 0.0)
    elseif score2 == 0
        # If score2 is 0, return first language, to prevent division by zero in the rate formula.
        # In this case confidence is calculated by another formula.
        # At this point there are two options:
        # * Text contains random characters that accidentally match trigrams of one of the languages
        # * Text really matches one of the languages.
        #
        # Number 500.0 is based on experiments and common sense expectations.
        confidence = score1 / 500.0
        if confidence > 1.0
            confidence = 1.0
        end
        return (lang_dist1[1], confidence)
    end

    rate = (score1 - score2) / score2

    # Hyperbola function. Everything that is above the function has confidence = 1.0
    # If rate is below, confidence is calculated proportionally.
    # Numbers 12.0 and 0.05 are obtained experimentally, so the function represents common sense.
    #
    confident_rate = (12.0 / length(trigrams)) + 0.05
    confidence = (rate > confident_rate) ? 1.0 : rate / confident_rate

    return (lang_dist1[1], confidence)
end

function calculate_distance(lang_trigrams,  text_trigrams)
    total_dist = 0
    for (i, trigram) in enumerate(lang_trigrams)
        dist = get(text_trigrams, trigram, MAX_TRIGRAM_DISTANCE+i) - i
        total_dist += abs(dist)
    end
    total_dist
end

Base.@deprecate detect(text::AbstractString, options=default_options()) LanguageDetector()(text, options)

mutable struct LanguageDetector
end

function(m::LanguageDetector)(text::AbstractString, options=default_options())
    if text==""; throw(ArgumentError("Cannot detect language for empty text")); end
    script = detect_script(text)
    if script == nothing; return (nothing, nothing, 0); end
    lang, conf = detect_lang_based_on_script(text, script, options)
    return (from_code(lang), script, conf)
end


text:
# Artificial Viscosity
#push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))
#using PumiInterface # pumi interface

@doc """
### artificialViscosity

It is a function that adds an artificial viscosity component to the weak form
of the PDE.

**Arguments**

*  qbar: Same as q except the last dof at a node is replaced by enthalpy
*  epsilonHat: Applied artificial viscosity
*  hi : Bounding box dimensions in a particular direction
*  h  : aithmetic mean of the components of the bounding box
*  Fav : Flux corresponding to artificial viscosity

"""->
function artificialViscosity(mesh::AbstractMesh{Tmsh},
                       sbp::AbstractOperator,
                       eqn::EulerData{Tsol, Tres, Tdim}) where {Tmsh,Tsol, Tres, Tdim}

  # Create the Artificial Viscosity flux matrix
  qbar = zeros(Tsol, mesh.numDofPerNode, sbp.numnodes, mesh.numEl)
  epsilonHat = 0.01
  hi = 1
  h = 1

  # Populate qbar
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      q_vals = sview(eqn.q, :, j, i)
      calcArtViscosityFluxComp(eqn.params,q_vals,sview(qbar,:,j,i))
    end
  end

  # We now need to differentiate the above with ξ and η to get the actual
  # artificial viscosity fluxes

  # Intermediate variables which have been multiplied with mapping jacobian
  Fav = zeros(Tsol, mesh.numDofPerNode, sbp.numnodes, mesh.numEl, Tdim)
  phi = zeros(Fav) # Intermediate variable

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      phi[:,j,i,1] = qbar[:,j,i] # abs(mesh.dxidx[1,1,j,i])*qbar[:,j,i] # compute dxi/dx*qbar
      phi[:,j,i,2] = qbar[:,j,i] # abs(mesh.dxidx[2,1,j,i])*qbar[:,j,i] # compute deta/dx*qbar
    end
  end

  for k = 1:Tdim # compute dqbar/dx
    # differentiate! differentiates wrt to ξ and η
    differentiate!(sbp,k,sview(phi,:,:,:,k),sview(Fav,:,:,:,1))
  end

  phi = zeros(Fav)
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      phi[:,j,i,1] = qbar[:,j,i] # abs(mesh.dxidx[1,2,j,i])*qbar[:,j,i] # compute dxi/dy*qbar
      phi[:,j,i,2] = qbar[:,j,i] # abs(mesh.dxidx[2,2,j,i])*qbar[:,j,i] # compute deta/dy*qbar
    end
  end

  for k = 1:Tdim # compute dqbar/dy
    differentiate!(sbp,k,sview(phi,:,:,:,k),sview(Fav,:,:,:,2))
  end

  Fav = -epsilonHat*(hi/h)*Fav

  for k = 1:Tdim # Assemble into the weak form
    weakdifferentiate!(sbp,k, sview(Fav,:,:,:,k), eqn.res, trans=true)
  end

return nothing

end

@doc """
### calcArtViscosityFluxComp

Computes the Artificial Viscosity flux components

**Inputs**

*  `params` :  ParameterType{2, :conservative}
*  `q`      :  Solution vector at a node
*  `F`      :  Artificial viscosity flux

**Outputs**

*  None

"""->

function calcArtViscosityFluxComp(params::ParamType{2}, q::AbstractArray{Tsol,1},
                                  F::AbstractArray{Tsol,1}) where Tsol

  press = calcPressure(params, q)
  F[1] = q[1]
  F[2] = q[2]
  F[3] = q[3]
  F[4] = q[4] #+ press

return nothing
end

#=
function calcEpsilonHat{Tsol}(params::ParamType{2}, q::AbstractArray{Tsol,1},
                           epsilon::AbstractArray{Tsol,1})
  # This function works at the nodal level

  # Calculate the temperature
  T = (q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1])*(1/(q[1]*params.cv))
  c = sqrt(params.gamma*params.R*T)  # Speed of sound
  p = q[1]*params.R*T  # Pressure
  hBar = 1  # Arithmetic mean of element size metrics

  # Maximum Eigen value of the system
  lambdamax = sqrt((q[2]*q[2] + q[3]*q[3])/(q[1]^2)) + c

  thetah = lambdamax*hBar/p
  thetal = 0.01*thetah

  # Calculate the elsilonhat value
  if epsilon <= thetal
    epsilonHat = 0
  elseif epsilon > thetal && epsilon <= thetah
    epsilonHat = 0.5*thetah(1 + sin(pi*(epsilon - thetal)/(thetah - thetal) - 0.5))
  else
    epsilonHat = thetah
  end

return nothing
end

function EpsilonPDE{Tmsh,Tsol, Tres, Tdim}(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                                     eqn::EulerData{Tsol, Tres, Tdim})

  # epsilon = zeros(Tsol, 1, sbp.numnodes, mesh.numEl)
  epsilon = zeros(Tsol, sbp.numnodes, mesh.numEl)

return nothing
end

function shockIndicator{Tsol}(params::ParamType{2}, q::AbstractArray{Tsol,1})
# Operates at the Nodal Level

  p = calcPressure(params, q)  # Calculate Pressure

  deltaPsi = 0.5             # Empirical constants that determine when when the
  psi0 = -4 - 4.25*log10(p)  # shock indicator should take effect.

  thetaS = 1  # Maximum Value

  # Calculate F: Discontinuity sensor

  if F <= psi0 - deltaPsi
    Sk = 0
  elseif F >= psi0 + deltaPsi
    Sk = thetaS
  elseif abs(F - psi0) < deltaPsi
    Sk = 0.5*thetaS*(1 + sin(0.5*pi*(F-psi0)/deltaPsi))
  end

end
=#

function boundingBox(coord::AbstractArray{Tmsh,2}, h::AbstractArray{Tmsh,1}) where Tmsh
  # It works at the element level. Accepts the vertex coordinates of all the nodes
  # h = zeros(Tmsh,Tdim)  # Stores the dimensions of the bounding box
  for i = 1:2
    h[i] = abs(maximum(coord[i,1:3]) - minimum(coord[i,1:3]))
  end

end

#=
function AverageMeshDimen{Tmsh, Tres, Tdim}(mesh::AbstractMesh{Tmsh},
                                      hBar::AbstractArray{Tmsh,1},
                                      eqn::EulerData{Tsol, Tres, Tdim})
  # Operates at the element level
  # Reset the face iterator at the beginning of wherever this loop is called
  element = getFace()
  (vertices, nvertices) = getDownward(mesh.m_ptr, element, 0)

  #hBar = zeros(Tmsh,nvertices)

  for i = 1:nvertices
    nadjElements = countAdjacent(mesh.m_ptr, vertices, 2) # Get number of adjacent elements
    adjElements = getAdjacent(nadjElements)  # Array of adjacent element pointers
    hArray = zeros(Tmsh, Tdim,nadjElements)  # Array of bounding box of all adjacent elements
    for j = 1:nadjElements
      elemNum = getNumberJ(mesh.el_Nptr, adjElements, 0, 0) # Get the global element number
      coordinate = mesh.coords[:,:,elemNum]  # get element nodal coordinates
      boundingBox(coordinate, sview(hArray,:,j)) # Get the bounding box
    end
    hBar[i] = sum(hArray)/(Tdim*nadjacencies) # Calculate the average mesh metric for the three vertices
  end

end =#

function AvgMeshSize(mesh::AbstractMesh{Tmsh}, eqn::EulerData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}

  hArray = zeros(Tmsh, Tdim,mesh.numEl)  # Array of bounding box of all adjacent elements
  resetFaceIt() # Reset Iterator over Face
  for i = 1:mesh.numEl
    element = getFace()
    (vertices, nvertices) = getDownward(mesh.m_ptr, element, 0)
    elemNum = getNumberJ(mesh.el_Nptr, element, 0, 0) # Get the global element number
    nodal_coordinates = mesh.coords[:,:,elemNum+1]  # get element nodal coordinates
    boundingBox(nodal_coordinates, sview(hArray,:,i)) # Get the bounding box
    incrementFaceIt() # Increment face iterator
  end
  hAverage = sum(hArray)/(Tdim*mesh.numEl)

  return hAverage
end

#=
function AVSourceTerm{Tmsh, Tsol, Tres, Tdim}(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim})

epsilon
end =#

ext: jl


text:
@symbol_func function K_W(cur_reactor::AbstractReactor)
  cur_kappa = kappa_P(cur_reactor)

  cur_epsilon = epsilon_P(cur_reactor)

  cur_denominator = cur_kappa ^ 2 - 1.0

  cur_denominator *= ( 2 / cur_reactor.pi )

  cur_denominator += 1

  cur_K = K_F(cur_reactor)

  cur_K /= 5.0

  cur_K /= cur_reactor.pi ^ 2

  cur_K *= cur_kappa

  cur_K /= cur_epsilon

  cur_K /= cur_denominator

  cur_K
end


text:
export readgdv, writegdv

"""
    readgdv(fd::IO, nodes::AbstractVector)
    readgdv(file::AbstractString, nodes::AbstractVector)

Reads the .ndump2 file format that contains (static or dynamic)
graphlet counts, outputted by GraphCrunch1 (ncount program in
http://www0.cs.ucl.ac.uk/staff/natasa/graphcrunch/index.html), or
Graphcrunch2
(http://www0.cs.ucl.ac.uk/staff/natasa/graphcrunch2/index.html), or
the dynamic graphlets counting code (https://www3.nd.edu/~cone/DG/).

Graphlets are small, connected, induced sub-graphs of a network
(Przulj N, Corneil DG, Jurisica I: Modeling Interactome, Scale-Free or Geometric?, Bioinformatics 2004, 20(18):3508-3515.),
similar to network motifs.
"""
function readgdv(fd::IO, nodesidx::Dict)
    A = readdlm(fd,Any)
    X = zeros(Float64, (size(A,1),size(A,2)-1))
    X[map(i -> nodesidx[string(A[i,1])], 1:size(A,1)),:] = A[:,2:end]
    X
end
readgdv(fd::IO, nodes::AbstractVector) =
    readgdv(fd, indexmap(nodes))
readgdv(file::AbstractString,args...) =
    open(fd -> readgdv(fd,args...), file, "r")

"""
    writegdv(fd::IO, X::AbstractMatrix, nodes::AbstractVector)
    writegdv(file::AbstractString, X::AbstractMatrix, nodes::AbstractVector)

Writes to graphlets file format. See [`readgdv`](@ref).
"""
function writegdv(fd::IO, X::AbstractMatrix, nodes::AbstractVector)
    for i = 1:size(X,1)
        print(fd, nodes[i], " ")
        for j = 1:size(X,2)
            if j < size(X,2)
                print(fd, Int(X[i,j]), " ")
            else
                println(fd, Int(X[i,j]))
            end
        end
    end
end
writegdv(file::AbstractString, args...) =
    open(fd -> writegdv(fd, args...), file, "w")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule VegaFusion_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("VegaFusion")
JLLWrappers.@generate_main_file("VegaFusion", UUID("e97bca81-1c03-5bcd-840e-fab7961fc59f"))
end  # module VegaFusion_jll


text:
using Test
using LinearAlgebra
using Printf
using MFGnet

@testset "NN" begin
    include("testNN.jl")
end

@testset "param2vec" begin
    include("testParam2Vec.jl")
end

@testset "PotentialNN" begin
    include("testPotentialNN.jl")
end

@testset "PotentialResNN" begin
    include("testPotentialResNN.jl")
end

@testset "PotentialSingle" begin
    include("testPotentialSingle.jl")
end

@testset "ResNN" begin
    include("testResNN.jl")
end

@testset "singleLayer" begin
    include("testSingleLayer.jl")
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using DataStructures: DefaultDict
using Distributions
using DynamicPPL
using Turing.RandomMeasures

export conditionals, sampled_values

"""
    abstract type Cont end

`Cont` is an analytic representation of joint distribution represented in a `Graph` -- a function
from a dictionary of variable assignments to a log probability.  (For lack of a better term, I call
these functions "continuations".)

Each `Call` gets associated with a `Transformation`, and each tilde statement with a
`LogLikelihood`.  Those are both callable with a dict argument.  SSA variables get transformed to
either a `Fixed` value for constants, a `Variable` for assignments of random variables, or the
`Cont` they come from.

Something like this:

```
⟨28⟩ => logpdf(Normal(), θ[μ[2]])
⟨29⟩ => getindex(⟨8⟩(array initializer with undefined values, 2), getfield(iterate(Colon()(1, 2), getfield(iterate(Colon()(1, 2)), 2)), 1))
⟨31⟩ => logpdf(Dirichlet(2, 1.0), θ[w])
⟨32⟩ => Array{Int64,1}(array initializer with undefined values, length([0.1, -0.05, 1.0]))
⟨33⟩ => Colon()(1, length([0.1, -0.05, 1.0]))
⟨34⟩ => iterate(Colon()(1, length([0.1, -0.05, 1.0])))
⟨36⟩ => getfield(iterate(Colon()(1, length([0.1, -0.05, 1.0]))), 1)
⟨37⟩ => getfield(iterate(Colon()(1, length([0.1, -0.05, 1.0]))), 2)
⟨42⟩ => logpdf(DiscreteNonParametric(θ[w]), θ[z[1]])
```

for

```
⟨28⟩ = μ[2] ~ Normal() → -1.2107564627453093
⟨29⟩ = μ[2] = getindex(⟨9⟩, ⟨23⟩) → -1.2107564627453093
⟨31⟩ = w ~ Dirichlet(⟨4⟩, 1.0) → [0.7023731332410442, 0.2976268667589558]
⟨32⟩ = Array{Int64,1}(array initializer with undefined values, ⟨7⟩) → [139815315085536, 139815315085552, 139815315085568]
⟨33⟩ = Colon()(1, ⟨7⟩) → 1:3
⟨34⟩ = iterate(⟨33⟩) → (1, 1)
⟨36⟩ = getfield(⟨34⟩, 1) → 1
⟨37⟩ = getfield(⟨34⟩, 2) → 1
⟨42⟩ = z[1] ~ DiscreteNonParametric(⟨31⟩) → 2
```

where `θ` stands for the environment of random variable assignments.
"""
abstract type Cont end

abstract type ArgSpec <: Cont end

struct Fixed{T} <: ArgSpec
    value::T
end
Base.show(io::IO, arg::Fixed) = print(io, arg.value)
(arg::Fixed)(θ) = arg.value

struct Variable{TV<:VarName} <: ArgSpec
    vn::TV
end
Base.show(io::IO, arg::Variable) = print(io, "θ[", arg.vn, "]")
(arg::Variable)(θ) = _lookup(θ, arg.vn)

function _lookup(θ, varname)
    if haskey(θ, varname)
        result = getindex(θ, varname)
        return result
    else
        # in the case of looking up x[i] with stored x,
        # simply do it the slow way and check all elements
        for (vn, value) in θ
            if DynamicPPL.subsumes(vn, varname)
                result = foldl((x, i) -> getindex(x, i...),
                               DynamicPPL.getindexing(varname),
                               init=value)
                return result
            end
        end
        throw(BoundsError(θ, varname))
    end
end

struct Transformation{TF, N, TArgs<:NTuple{N, Cont}} <: Cont
    f::TF
    args::TArgs
end

function Base.show(io::IO, t::Transformation)
    print(io, t.f, "(")
    join(io, t.args, ", ")
    print(io, ")")
end

(t::Transformation)(θ) = t.f((deepcopy(arg(θ)) for arg in t.args)...)

function (t::Transformation{typeof(getindex)})(θ)
    # intercept this to simplify getindex(x, i) to direct lookup of x[i]
    array, indexing = first(t.args), Base.tail(t.args)
    if array isa Variable
        actual_vn = VarName(array.vn, (Tuple(ix(θ) for ix in indexing),))
        return _lookup(θ, actual_vn)
    else
        return t.f((deepcopy(arg(θ)) for arg in t.args)...)
    end
end

struct LogLikelihood{TDist, TF, TArgs, TVal} <: Cont
    dist::TDist
    f::TF # function that was used to construct the distribution
    args::TArgs
    value::TVal

    function LogLikelihood(dist::Distribution, f, args::NTuple{N, Cont}, value) where {N}
        return new{typeof(dist), typeof(f), typeof(args), typeof(value)}(
            dist, f, args, value)
    end
end

function Base.show(io::IO, ℓ::LogLikelihood{D}) where {D}
    print(io, "logpdf(", ℓ.f, "(")
    join(io, ℓ.args, ", ")
    print(io, "), ", ℓ.value, ")")
end

(ℓ::LogLikelihood{D})(θ) where {D} = begin
    logpdf(ℓ.f((arg(θ) for arg in ℓ.args)...), ℓ.value(θ))
end

_init(::Tuple{}) = ()
_init(t::Tuple{Any}) = ()
_init(t::Tuple) = (first(t), _init(Base.tail(t))...)
_last(::Tuple{}) = ()
_last((x,)::Tuple{Any}) = x
_last(t::Tuple) = _last(Base.tail(t))

"""
    continuations(graph)

Assign to each node in the graph a `Cont` object, representing a function from a variable assignment
(θ) to the value of an expression or the log-likelihood in case of a tilde statement.
"""
function continuations(graph)
    c = SortedDict{Reference, Cont}()

    function convertarg(arg)
        if arg isa Reference
            cont = c[arg]
            stmt = graph[arg]
            if cont isa LogLikelihood && !isnothing(stmt.vn)
                return Variable(stmt.vn)
            elseif cont isa Transformation{typeof(getindex)} && !isnothing(stmt.definition)
                # two cases:
                # - vn[ix] = getindex(<a>, ix) where <a> = vn ~ D: then we keep getindex(θ[vn], ix)
                #   (the whole array was sampled by the tilde)
                # - vn[ix] = getindex(<a>, ix) where <scalar> = vn[ix] ~ D: then we rewrite to
                #   `θ[vn[ix]]`  (the array element was sampled individually)
                # this truncation is realized by `_init`, which discards the last element of a tuple,
                # and maps the empty tuple to itself.
                vn, ref = stmt.definition
                ref_vn = graph[ref].vn
                array, ix = stmt.args[1], stmt.args[2:end]
                parent_vn = VarName(graph[ref].vn, _init(DynamicPPL.getindexing(graph[ref].vn)))
                parent_array = Variable(parent_vn)
                return Transformation(getindex, (parent_array, convertarg.(ix)...))
            else
                return cont
            end
        else
            return Fixed(arg)
        end
    end

    for (ref, stmt) in graph
        if stmt isa Union{Assumption, Observation}
            dist_stmt = graph[stmt.dist_ref]
            dist = getvalue(dist_stmt)
            if dist_stmt isa Call
                f, args = dist_stmt.f, convertarg.(dist_stmt.args)
            elseif f isa Constant
                # fake a constant here... `getindex(Ref(x)) == x`
                f = getindex
                args = (Fixed(Ref(dist)),)
            end
            value = Variable(stmt.vn)
            c[ref] = LogLikelihood(dist, f, args, value)
        elseif stmt isa Call
            f, args = stmt.f, convertarg.(stmt.args)
            if f isa Reference
                c[ref] = Transformation(((f, t...) -> f(t...)), (convertarg(f), args...))
            else
                c[ref] = Transformation(f, args)
            end
        elseif stmt isa Constant
            c[ref] = Fixed(deepcopy(getvalue(stmt)))
        end
    end

    return c
end

"""
    conditionals(graph, varname)

Calculate the `GibbsConditional`s of all variables in `graph` that match `varname`.  There can
be multiple observations for which this holds, so a dictionary is returned.
"""
function conditionals(graph, varname)
    dists = Dict{VarName, LogLikelihood}()
    blankets = DefaultDict{VarName, Vector{Pair{VarName, LogLikelihood}}}(
        Vector{Pair{Tuple, LogLikelihood}})
    conts = continuations(graph)

    for (ref, stmt) in graph
        if stmt isa Union{Assumption, Observation} && !isnothing(stmt.vn)
            # record distribution of this tilde if it matches the searched vn
            if DynamicPPL.subsumes(varname, stmt.vn)
                dists[stmt.vn] = conts[ref]
            end

            # add likelihood to all parents of which this RV is in the blanket
            for (pvn, p) in parent_variables(graph, stmt)
                for vn in keys(dists)
                    if DynamicPPL.subsumes(vn, p.vn)
                        push!(blankets[vn], pvn => conts[ref])
                        break
                    end
                end
            end
        end
    end

    return Dict(vn => GibbsConditional(vn, d, blankets[vn]) for (vn, d) in dists)
end

"""
    sampled_values(graph)

Extract the values of all observed and assumed random variables in the graph, including their
occuring parts.
"""
function sampled_values(graph)
    θ = Dict{VarName, Any}()
    for (ref, stmt) in graph
        if stmt isa Union{Assumption, Observation} && !isnothing(stmt.vn)
            θ[stmt.vn] = tovalue(graph, getvalue(stmt))
        elseif stmt isa Call && !isnothing(stmt.definition)
            # remember all intermediate RV values (including redundant `getindex` calls,
            # for simplicity)
            vn, _ = stmt.definition
            θ[vn] = getvalue(stmt)
        end
    end

    return θ
end

struct GibbsConditional{
    TVar<:VarName,
    TBase<:LogLikelihood,
    TBlanket}

    vn::TVar
    base::TBase
    blanket::TBlanket
end

function Base.show(io::IO, c::GibbsConditional)
    print(io, c.base)
    if !isempty(c.blanket)
        print(io, " + ")
        join(io, (β for (vn, β) in c.blanket), " + ")
    end
end

"""
    (c::GibbsConditional)(θ)

Return the conditional distribution of `vn` given the values fixed in `θ`, calculated by
normalization over the conditional distribution of `vn` itlself and its Markov blanket.

Constructed as

    P[X = x | θ] ∝ P[X = x | parents(X)] * P[children(X) | parents(children(X))]

equivalent to a distribution `D` such that

    logpdf(D, x) = ℓ(x | θ) + ∑ blanketᵢ(x | θ)

where the factors `blanketᵢ` are the log probabilities in the Markov blanket.

This only works on discrete distributions, either scalar ones (resulting in a
`DiscreteNonparametric`) or products of them (resulting in a `Product` of `DiscreteNonparametric`).
"""
function (c::GibbsConditional{V, L})(θ) where {
    V<:VarName, L<:LogLikelihood{<:DiscreteUnivariateDistribution}}

    Ω = try
        support(c.base.dist)
    catch
        throw(ArgumentError("Unable to get the support of $(c.base.dist) (probably infinite)!"))
    end

    θs_on_support = fixvalues(θ, c.vn => Ω)
    logtable = [c.base(θ′) + reduce(+, (β(θ′) for (vn, β) in c.blanket), init=0.0)
                for θ′ in θs_on_support]
    conditional = DiscreteNonParametric(Ω, softmax!(logtable))
    return conditional
end

# `Product`s can be treated as an array of iid variables
function (c::GibbsConditional{V, L})(θ) where {
    V<:VarName, L<:LogLikelihood{<:Product}}

    independent_distributions = c.base.dist.v
    Ωs = try
        support.(independent_distributions)
    catch
        throw(ArgumentError("Unable to get the support of $(c.base.dist) (probably infinite)!"))
    end

    conditionals = similar(Ωs, DiscreteNonParametric)

    for index in eachindex(Ωs, conditionals)
        sub_vn = DynamicPPL.VarName(c.vn, (DynamicPPL.getindexing(c.vn)..., (index,)))
        θs_on_support = fixvalues(θ, sub_vn => Ωs[index])
        logtable = map(θs_on_support) do θ′
            c.base(θ′) + reduce(+, (β(θ′) for (vn, β) in c.blanket if vn == sub_vn), init=0.0)
        end
        conditionals[index] = DiscreteNonParametric(Ωs[index], softmax!(vec(logtable)))
    end

    return Product(conditionals)
end

# Special treatment for CRP variables: calculate likelihoods as normal for truncated support
# (covering all existing clusters), and marginalize the creation of a new cluster
function (c::GibbsConditional{V, L})(θ) where {
    V<:VarName, L<:LogLikelihood{<:ChineseRestaurantProcess}}
    Ω = support(c.base.dist)
    Ω_init, Ω_last = Ω[1:end-1], Ω[end]
    # @show c.vn
    # @show Ω_init, Ω_last

    θs_on_init = fixvalues(θ, c.vn => Ω_init)

    # @show θs_on_init

    # t = [reduce(+, (β(θ′) for (vn, β) in c.blanket), init=0.0)
         # for θ′ in θs_on_init]
    # display(c.blanket)
    # @show [vn for (vn, β) in c.blanket for θ′ in θs_on_init]
    # println()

    logtable_init = Float64[c.base(θ′) + reduce(+, (β(θ′) for (vn, β) in c.blanket), init=0.0)
                            for θ′ in θs_on_init]

    θ_on_last = fixvalue(θ, c.vn => Ω_last)
    log_last = _estimate_last_likelihood(c, θ_on_last)
    conditional = DiscreteNonParametric(Ω, softmax!(push!(logtable_init, log_last)))
    return conditional
end

"""
Estimate the "new cluster" likelihood of a CRP mixture, given through

    𝓅(zₙ = K + 1 | z₁, ..., zₙ₋₁, μ, xₙ) ∝ (∏_{i = z ≥ n} 𝓅(zᵢ | z₁,...,zᵢ)) 𝓅(xₙ | zₙ = K + 1, μ),

by approximating

    𝓅(xₙ | zₙ = K + 1, μ) = ∫ 𝓅(xₙ, μ[zₙ] = m) dm ≈ 𝓅(xₙ | rand(m))

where Law(m) = Law(μ).
"""
function _estimate_last_likelihood(c, θ)
    l = c.base(θ)

    rpm = c.base.args[1](θ)
    if rpm isa DirichletProcess
        if length(c.base.args[1].args) == 2
            G₀ = c.base.args[1].args[2](θ)
        else
            error("Cannot handle Dirichlet process with constructor $(c.base.args[1])")
        end
    else
        error("Cannot handle CRP based on process $(c.base.args[1])")
    end

    for (vn, β) in c.blanket
        if β isa LogLikelihood{<:ChineseRestaurantProcess} && DynamicPPL.subsumes(c.vn, vn)
            # one of the CRP factors (the 𝓅(zᵢ | ...) for i > n)
            l += β(θ)
        else
            # the factor of the observed value, 𝓅(xₙ | zₙ == K + 1, μ),
            # which is what we are actually intested in in this function.
            # this is a real hack: we just assume that this vn must be the likelihood of
            # the observation, and this has distribution G₀. no checks are done.
            l += logpdf(G₀, θ[vn])
        end
    end

    return l
end

(c::GibbsConditional)(θ) =
    throw(ArgumentError("Cannot condition a non-discrete or non-univariate distribution $(c.base)."))

"""Produce `|Ω|` copies of `θ` with the `fixedvn` entries fixed to the values in the support `Ω`."""
function fixvalue(θ, (source_vn, value))
    # "source" is the value stored in Ω; "target" is the matching value in θ.

    # result = [copy(θ) for _ in Ω]
    θ′ = copy(θ)

    # for (θ′, value) in zip(result, Ω)
        for target_vn in keys(θ′)
            source_subsumes_target = DynamicPPL.subsumes(source_vn, target_vn)
            target_subsumes_source = DynamicPPL.subsumes(target_vn, source_vn)

            if source_subsumes_target && target_subsumes_source
                # target <- source, target[i] <- source[i]
                # both indices match -- we just update the complete thing
                θ′[target_vn] = value

            elseif target_subsumes_source
                # target <- source[i], target[i] <- source[i][j]
                # updating the target from a part of the source with "copy on write",
                # target[i][j] = source[i][j] <=> setindex!(copy(target[i]), source[i][j], j)

                # NB: we index into the target by the source index!
                target_indexing = DynamicPPL.getindexing(source_vn)
                ti_init, ti_last = _splitindexing(target_indexing)
                target = θ′[target_vn]
                initial_target = foldl((x, i) -> getindex(x, i...),
                                       ti_init,
                                       init=target)
                if isnothing(ti_last)
                    θ′[target_vn] = initial_target
                else
                    target = setindex!(copy(initial_target), value, ti_last...)
                    θ′[target_vn] = target
                end
            elseif source_subsumes_target
                # target[i][j] <- source[i]
                # setting the target to part of the source array
                # target[i][j] = source[i][j]

                target_indexing = DynamicPPL.getindexing(target_vn)
                target = foldl((x, i) -> getindex(x, i...),
                               target_indexing,
                               init=value)
                θ′[target_vn] = target
            end
        end
    # end

    return θ′
end

fixvalues(θ, (source_vn, Ω)) = [fixvalue(θ, source_vn => ω) for ω in Ω]

_splitindexing(::NTuple{0}) = (), nothing
_splitindexing((i,)::NTuple{1}) = (), i
_splitindexing((i, j)::NTuple{2}) = (i,), j
_splitindexing(t::Tuple) = let (init, last) = _splitindexing(Base.tail(t))
    (first(t), init...), last
end

# from https://github.com/JuliaStats/StatsFuns.jl/blob/master/src/basicfuns.jl#L259
function softmax!(x::AbstractArray{<:AbstractFloat})
    u = maximum(x)
    s = zero(u)

    @inbounds for i in eachindex(x)
        s += (x[i] = exp(x[i] - u))
    end

    s⁻¹ = inv(s)

    @inbounds for i in eachindex(x)
        x[i] *= s⁻¹
    end

    return x
end

softmax(x::AbstractArray{<:AbstractFloat}) = softmax!(copy(x))

# D_w = Dirichlet(2, 1.0)
# w = rand(D_w)
# D_p = DiscreteNonParametric([0.3, 0.6], w)
# p = rand(D_p)
# D_x = Bernoulli(p)
# x = rand(D_x)
# conditioned(D_p, [logpdf(D_x, x)])


text:
using Makie, YOLO

#First time only (downloads 5011 images & labels!)
#YOLO.download_dataset("voc2007")

settings = YOLO.pretrained.v2_tiny_voc.load(minibatch_size=1) #Process one image at a time
model = YOLO.v2_tiny.load(settings)
YOLO.loadWeights!(model, settings)

voc = YOLO.datasets.VOC.populate()
vocloaded = YOLO.load(voc, settings, indexes = [rand(1:5011)]) #Load a single random image

#Run the model
res = model(vocloaded.imstack_mat);

#Convert the output into readable predictions
predictions = YOLO.postprocess(res, settings, conf_thresh = 0.3, iou_thresh = 0.3)

scene = YOLO.renderResult(vocloaded.imstack_mat[:,:,:,1], predictions, settings, save_file = "test.png")
display(scene)


text:
using LinearAlgebra
using Combinatorics
using QuantumInformation
using SparseArrays

function superoperator(N, atol=1e-3)
    """
        Returns the superoperator of the channel:
            ϕ(Y) = ∫_U dU (U* ⊗ U^{⊗N}) Y (U* ⊗ U^{⊗N})^†
    """

    P = zeros(2^(2N+2), 2^(2N+2))
    idx = CartesianIndices(Tuple(fill(2, N+1)))

    for perm in permutations(1:N+1)
        x = zeros(fill(2, 2N+2)...)
        for k1 in idx
            k0 = CartesianIndex(Tuple(k1)[perm])
            x[k0, k1] = 1
        end
        x = vec(x)
        P += x * x'
    end

    F = eigen(P)
    Λ = Diagonal([x > atol ? 1 : 0 for x in F.values])
    P = F.vectors * Λ * F.vectors'

    sparse(permutesystems(P, [2, 2^N, 2, 2^N], [3, 2, 1, 4]))
end

function diff_q_op(N, atol=1e-03)
    """
        Returns the dynamical matrix of the difference of channels:
            ϕ_0(X) = ∫_U dU (U* ⊗ U^{⊗N}) (|0><0| ⊗ X) (U* ⊗ U^{⊗N})^†
            ϕ_1(X) = ∫_U dU (U* ⊗ U^{⊗N}) (|1><1| ⊗ X) (U* ⊗ U^{⊗N})^†
    """

    S = superoperator(N, atol)
    Φ_0 = S * (ket(1,2) ⊗ I(2^N) ⊗ ket(1,2) ⊗ I(2^N))
    Φ_1 = S * (ket(2,2) ⊗ I(2^N) ⊗ ket(2,2) ⊗ I(2^N))
    return DynamicalMatrix(reshuffle(Φ_0 - Φ_1), 2^N, 2^(N+1))
end

function calculate_p_opt(N, atol=1e-03)
    """
        Calulate the maximal value of probability (p_opt) of discrimination of
        channels:
            ϕ_0(X) = ∫_U dU (U* ⊗ U^{⊗N}) (|0><0| ⊗ X) (U* ⊗ U^{⊗N})^†
            ϕ_1(X) = ∫_U dU (U* ⊗ U^{⊗N}) (|1><1| ⊗ X) (U* ⊗ U^{⊗N})^†

        p_opt = 0.5 + 0.25 * |||ϕ_0 - ϕ_1|||,
        where ||| ⋅ ||| is the diamond norm.
    """

    1/2 + 1/4 * norm_diamond(diff_q_op(N), Val(:primal), atol)
end

function calculate_p_opt_dual(N, atol=1e-03)
    """
        Calulate the dual value for calculate_p_opt(N)
    """

    1/2 + 1/4 * norm_diamond(diff_q_op(N), Val(:dual), atol)
end

N = 1
print(calculate_p_opt_dual(N), " ≤ p_opt ≤ ", calculate_p_opt(N))

max_stars_repo_name: iitis/QML-pm-fidelity

max_issues_repo_name: iitis/QML-pm-fidelity

max_forks_repo_name: iitis/QML-pm-fidelity


text:
using Cosmology

default_cosmology() = cosmology(h=0.70, OmegaM=0.3)   #S11


text:
#==============================================================================#
# Mobile.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/mobile-2017-07-01.normal.json
#==============================================================================#

__precompile__()

module Mobile

using AWSCore

"""
    using AWSSDK.Mobile.create_project
    create_project([::AWSConfig], arguments::Dict)
    create_project([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "POST", "/projects", arguments::Dict)
    mobile([::AWSConfig], "POST", "/projects", <keyword arguments>)

# CreateProject Operation

Creates an AWS Mobile Hub project.

# Arguments

## `name = ::String`
Name of the project.

## `region = ::String`
Default region where project resources should be created.

## `contents = blob`
ZIP or YAML file which contains configuration settings to be used when creating the project. This may be the contents of the file downloaded from the URL provided in an export project operation.

## `snapshotId = ::String`
Unique identifier for an exported snapshot of project configuration. This snapshot identifier is included in the share URL when a project is exported.

# Returns

`CreateProjectResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException`, `BadRequestException`, `NotFoundException` or `LimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/CreateProject)
"""
@inline create_project(aws::AWSConfig=default_aws_config(); args...) = create_project(aws, args)

@inline create_project(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "POST", "/projects", args)

@inline create_project(args) = create_project(default_aws_config(), args)

"""
    using AWSSDK.Mobile.delete_project
    delete_project([::AWSConfig], arguments::Dict)
    delete_project([::AWSConfig]; projectId=)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "DELETE", "/projects/{projectId}", arguments::Dict)
    mobile([::AWSConfig], "DELETE", "/projects/{projectId}", projectId=)

# DeleteProject Operation

Delets a project in AWS Mobile Hub.

# Arguments

## `projectId = ::String` -- *Required*
Unique project identifier.

# Returns

`DeleteProjectResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/DeleteProject)
"""
@inline delete_project(aws::AWSConfig=default_aws_config(); args...) = delete_project(aws, args)

@inline delete_project(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "DELETE", "/projects/{projectId}", args)

@inline delete_project(args) = delete_project(default_aws_config(), args)

"""
    using AWSSDK.Mobile.describe_bundle
    describe_bundle([::AWSConfig], arguments::Dict)
    describe_bundle([::AWSConfig]; bundleId=)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "GET", "/bundles/{bundleId}", arguments::Dict)
    mobile([::AWSConfig], "GET", "/bundles/{bundleId}", bundleId=)

# DescribeBundle Operation

Get the bundle details for the requested bundle id.

# Arguments

## `bundleId = ::String` -- *Required*
Unique bundle identifier.

# Returns

`DescribeBundleResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException`, `BadRequestException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/DescribeBundle)
"""
@inline describe_bundle(aws::AWSConfig=default_aws_config(); args...) = describe_bundle(aws, args)

@inline describe_bundle(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "GET", "/bundles/{bundleId}", args)

@inline describe_bundle(args) = describe_bundle(default_aws_config(), args)

"""
    using AWSSDK.Mobile.describe_project
    describe_project([::AWSConfig], arguments::Dict)
    describe_project([::AWSConfig]; projectId=, <keyword arguments>)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "GET", "/project", arguments::Dict)
    mobile([::AWSConfig], "GET", "/project", projectId=, <keyword arguments>)

# DescribeProject Operation

Gets details about a project in AWS Mobile Hub.

# Arguments

## `projectId = ::String` -- *Required*
Unique project identifier.

## `syncFromResources = ::Bool`
If set to true, causes AWS Mobile Hub to synchronize information from other services, e.g., update state of AWS CloudFormation stacks in the AWS Mobile Hub project.

# Returns

`DescribeProjectResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException`, `BadRequestException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/DescribeProject)
"""
@inline describe_project(aws::AWSConfig=default_aws_config(); args...) = describe_project(aws, args)

@inline describe_project(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "GET", "/project", args)

@inline describe_project(args) = describe_project(default_aws_config(), args)

"""
    using AWSSDK.Mobile.export_bundle
    export_bundle([::AWSConfig], arguments::Dict)
    export_bundle([::AWSConfig]; bundleId=, <keyword arguments>)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "POST", "/bundles/{bundleId}", arguments::Dict)
    mobile([::AWSConfig], "POST", "/bundles/{bundleId}", bundleId=, <keyword arguments>)

# ExportBundle Operation

Generates customized software development kit (SDK) and or tool packages used to integrate mobile web or mobile app clients with backend AWS resources.

# Arguments

## `bundleId = ::String` -- *Required*
Unique bundle identifier.

## `projectId = ::String`
Unique project identifier.

## `platform = "OSX", "WINDOWS", "LINUX", "OBJC", "SWIFT", "ANDROID" or "JAVASCRIPT"`
Developer desktop or target application platform.

# Returns

`ExportBundleResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException`, `BadRequestException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/ExportBundle)
"""
@inline export_bundle(aws::AWSConfig=default_aws_config(); args...) = export_bundle(aws, args)

@inline export_bundle(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "POST", "/bundles/{bundleId}", args)

@inline export_bundle(args) = export_bundle(default_aws_config(), args)

"""
    using AWSSDK.Mobile.export_project
    export_project([::AWSConfig], arguments::Dict)
    export_project([::AWSConfig]; projectId=)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "POST", "/exports/{projectId}", arguments::Dict)
    mobile([::AWSConfig], "POST", "/exports/{projectId}", projectId=)

# ExportProject Operation

Exports project configuration to a snapshot which can be downloaded and shared. Note that mobile app push credentials are encrypted in exported projects, so they can only be shared successfully within the same AWS account.

# Arguments

## `projectId = ::String` -- *Required*
Unique project identifier.

# Returns

`ExportProjectResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException`, `BadRequestException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/ExportProject)
"""
@inline export_project(aws::AWSConfig=default_aws_config(); args...) = export_project(aws, args)

@inline export_project(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "POST", "/exports/{projectId}", args)

@inline export_project(args) = export_project(default_aws_config(), args)

"""
    using AWSSDK.Mobile.list_bundles
    list_bundles([::AWSConfig], arguments::Dict)
    list_bundles([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "GET", "/bundles", arguments::Dict)
    mobile([::AWSConfig], "GET", "/bundles", <keyword arguments>)

# ListBundles Operation

List all available bundles.

# Arguments

## `maxResults = ::Int`
Maximum number of records to list in a single response.

## `nextToken = ::String`
Pagination token. Set to null to start listing bundles from start. If non-null pagination token is returned in a result, then pass its value in here in another request to list more bundles.

# Returns

`ListBundlesResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException` or `BadRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/ListBundles)
"""
@inline list_bundles(aws::AWSConfig=default_aws_config(); args...) = list_bundles(aws, args)

@inline list_bundles(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "GET", "/bundles", args)

@inline list_bundles(args) = list_bundles(default_aws_config(), args)

"""
    using AWSSDK.Mobile.list_projects
    list_projects([::AWSConfig], arguments::Dict)
    list_projects([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "GET", "/projects", arguments::Dict)
    mobile([::AWSConfig], "GET", "/projects", <keyword arguments>)

# ListProjects Operation

Lists projects in AWS Mobile Hub.

# Arguments

## `maxResults = ::Int`
Maximum number of records to list in a single response.

## `nextToken = ::String`
Pagination token. Set to null to start listing projects from start. If non-null pagination token is returned in a result, then pass its value in here in another request to list more projects.

# Returns

`ListProjectsResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException` or `BadRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/ListProjects)
"""
@inline list_projects(aws::AWSConfig=default_aws_config(); args...) = list_projects(aws, args)

@inline list_projects(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "GET", "/projects", args)

@inline list_projects(args) = list_projects(default_aws_config(), args)

"""
    using AWSSDK.Mobile.update_project
    update_project([::AWSConfig], arguments::Dict)
    update_project([::AWSConfig]; projectId=, <keyword arguments>)

    using AWSCore.Services.mobile
    mobile([::AWSConfig], "POST", "/update", arguments::Dict)
    mobile([::AWSConfig], "POST", "/update", projectId=, <keyword arguments>)

# UpdateProject Operation

Update an existing project.

# Arguments

## `contents = blob`
ZIP or YAML file which contains project configuration to be updated. This should be the contents of the file downloaded from the URL provided in an export project operation.

## `projectId = ::String` -- *Required*
Unique project identifier.

# Returns

`UpdateProjectResult`

# Exceptions

`InternalFailureException`, `ServiceUnavailableException`, `UnauthorizedException`, `TooManyRequestsException`, `BadRequestException`, `NotFoundException`, `AccountActionRequiredException` or `LimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mobile-2017-07-01/UpdateProject)
"""
@inline update_project(aws::AWSConfig=default_aws_config(); args...) = update_project(aws, args)

@inline update_project(aws::AWSConfig, args) = AWSCore.Services.mobile(aws, "POST", "/update", args)

@inline update_project(args) = update_project(default_aws_config(), args)

end # module Mobile

#==============================================================================#
# End of file
#==============================================================================#


text:
# for pseudo_basis, basis_pmatrix, etc. see NfRel/NfRelOrdIdl.jl

################################################################################
#
#  Basic field access
#
################################################################################

@doc Markdown.doc"""
    order(a::NfRelOrdFracIdl) -> NfRelOrd

Returns the order of $a$.
"""
order(a::NfRelOrdFracIdl) = a.order

@doc Markdown.doc"""
    nf(a::NfRelOrdFracIdl) -> NumField

Returns the number field, of which $a$ is an fractional ideal.
"""
nf(a::NfRelOrdFracIdl) = nf(order(a))

################################################################################
#
#  Parent
#
################################################################################

parent(a::NfRelOrdFracIdl) = a.parent

################################################################################
#
#  iszero/isone
#
################################################################################

iszero(a::NfRelOrdFracIdl) = iszero(basis_matrix(a, copy = false)[1, 1])

function isone(a::NfRelOrdFracIdl)
  if denominator(a) != 1
    return false
  end
  @assert isone(basis_pmatrix(a, copy = false).matrix[1, 1])
  @assert isone(basis_pmatrix(order(a), copy = false).matrix[1, 1])

  return isone(basis_pmatrix(a, copy = false).coeffs[1])
end

################################################################################
#
#  Numerator and Denominator
#
################################################################################

function assure_has_denominator(a::NfRelOrdFracIdl)
  if isdefined(a, :den)
    return nothing
  end
  if iszero(a)
    a.den = fmpz(1)
    return nothing
  end
  O = order(a)
  n = degree(O)
  PM = basis_pmatrix(a, copy = false)
  pb = pseudo_basis(O, copy = false)
  inv_coeffs = inv_coeff_ideals(O, copy = false)
  d = fmpz(1)
  for i = 1:n
    for j = 1:i
      d = lcm(d, denominator(simplify(PM.matrix[i, j]*PM.coeffs[i]*inv_coeffs[j])))
    end
  end
  a.den = d
  return nothing
end

@doc Markdown.doc"""
    denominator(a::NfRelOrdFracIdl) -> fmpz

Returns the smallest positive integer $d$ such that $da$ is contained in
the order of $a$.
"""
function denominator(a::NfRelOrdFracIdl; copy::Bool = true)
  assure_has_denominator(a)
  if copy
    return deepcopy(a.den)
  else
    return a.den
  end
end

@doc Markdown.doc"""
    numerator(a::NfRelOrdFracIdl) -> NfRelOrdIdl

Returns the ideal $d*a$ where $d$ is the denominator of $a$.
"""
function numerator(a::NfRelOrdFracIdl; copy::Bool = true) # copy for compatibility with NfOrdFracIdl (it doesn't do anything here)
  d = denominator(a)
  PM = basis_pmatrix(a)
  if isone(d)
    return ideal_type(order(a))(order(a), PM)
  end
  for i = 1:degree(order(a))
    PM.coeffs[i] = PM.coeffs[i]*d
    PM.coeffs[i] = simplify(PM.coeffs[i])
  end
  return ideal_type(order(a))(order(a), PM)
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, s::NfRelOrdFracIdlSet)
  print(io, "Set of fractional ideals of ")
  print(io, s.order)
end

function show(io::IO, a::NfRelOrdFracIdl)
  compact = get(io, :compact, false)
  if compact
    print(io, "Fractional ideal")
    #show(IOContext(io, :compact => true), basis_pmatrix(a, copy = false))
  else
    print(io, "Fractional ideal of\n")
    show(IOContext(io, :compact => true), order(a))
    print(io, "\nwith basis pseudo-matrix\n")
    show(IOContext(io, :compact => true), basis_pmatrix(a, copy = false))
  end
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    fractional_ideal(O::NfRelOrd, M::PMat, M_in_hnf::Bool = false) -> NfRelOrdFracIdl

Creates the fractional ideal of $\mathcal O$ with basis pseudo-matrix $M$. If
M_in_hnf is set, then it is assumed that $M$ is already in lower left pseudo
HNF.
"""
function fractional_ideal(O::NfRelOrd{T, S, U}, M::PMat{T, S}, M_in_hnf::Bool = false) where {T, S, U}
  !M_in_hnf && !iszero(matrix(M)) ? M = pseudo_hnf(M, :lowerleft, true) : nothing
  return NfRelOrdFracIdl{T, S, U}(O, M)
end

@doc Markdown.doc"""
    fractional_ideal(O::NfRelOrd, M::Generic.Mat) -> NfRelOrdFracIdl

Creates the fractional ideal of $\mathcal O$ with basis matrix $M$.
"""
function fractional_ideal(O::NfRelOrd{T, S, U}, M::Generic.Mat{T}) where {T, S, U}
  coeffs = deepcopy(basis_pmatrix(O, copy = false)).coeffs
  return fractional_ideal(O, PseudoMatrix(M, coeffs))
end

function fractional_ideal(O::NfRelOrd{T, S, U}, A::Vector{U}) where {T, S, U}
  if all(iszero, A)
    M = zero_matrix(base_field(nf(O)), degree(O), degree(O))
    pb = pseudo_basis(O)
    return NfRelOrdFracIdl{T, S, U}(O, PseudoMatrix(M, [ deepcopy(pb[i][2]) for i = 1:degree(O)]))
  end

  return sum(fractional_ideal(O, a) for a in A if !iszero(a))
end

function fractional_ideal(O::NfRelOrd{T, S, U}, x::U) where {T, S, U}
  d = degree(O)
  pb = pseudo_basis(O, copy = false)
  M = zero_matrix(base_field(nf(O)), d, d)
  if iszero(x)
    return NfRelOrdFracIdl{T, S, U}(O, PseudoMatrix(M, [ deepcopy(pb[i][2]) for i = 1:d ]))
  end
  for i = 1:d
    elem_to_mat_row!(M, i, pb[i][1]*x)
  end
  M = M*basis_mat_inv(O, copy = false)
  PM = PseudoMatrix(M, [ deepcopy(pb[i][2]) for i = 1:d ])
  PM = pseudo_hnf(PM, :lowerleft)
  return NfRelOrdFracIdl{T, S, U}(O, PM)
end

*(O::NfRelOrd{T, S, U}, x::U) where {T, S, U} = fractional_ideal(O, x)

*(x::U, O::NfRelOrd{T, S, U}) where {T, S, U} = fractional_ideal(O, x)

function fractional_ideal(O::NfRelOrd{T, S, U}, a::NfRelOrdIdl{T, S, U}) where {T, S, U}
  return fractional_ideal(O, basis_pmatrix(a), true)
end

function fractional_ideal(O::NfRelOrd{T, S}, a::NfRelOrdIdl{T, S}, d::U) where { T, S, U <: Union{ fmpz, NfAbsOrdElem, NfRelOrdElem } }
  K = base_field(nf(O))
  dd = inv(K(d))
  return fractional_ideal(O, dd*basis_pmatrix(a), true)
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::NfRelOrdFracIdl{T, S, U}, dict::IdDict) where {T, S, U}
  z = NfRelOrdFracIdl{T, S, U}(a.order)
  for x in fieldnames(typeof(a))
    if x != :order && x != :parent && isdefined(a, x)
      setfield!(z, x, Base.deepcopy_internal(getfield(a, x), dict))
    end
  end
  z.order = a.order
  z.parent = a.parent
  return z
end

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(a::NfOrdRelFracIdl, b::NfRelOrdFracIdl) -> Bool

Returns whether $a$ and $b$ are equal.
"""
function ==(a::NfRelOrdFracIdl, b::NfRelOrdFracIdl)
  order(a) !== order(b) && return false
  return basis_pmatrix(a, copy = false) == basis_pmatrix(b, copy = false)
end

function ==(a::NfRelOrdIdl, b::NfRelOrdFracIdl)
  order(a) !== order(b) && return false
  return basis_pmatrix(a, copy = false) == basis_pmatrix(b, copy = false)
end
==(a::NfRelOrdFracIdl, b::NfRelOrdIdl) = b == a

################################################################################
#
#  Norm
#
################################################################################

# Assumes, that det(basis_matrix(a)) == 1
function assure_has_norm(a::NfRelOrdFracIdl)
  if a.has_norm
    return nothing
  end
  if iszero(a)
    O = order(basis_pmatrix(a, copy = false).coeffs[1])
    a.norm = nf(O)()*O
    a.has_norm = true
    return nothing
  end
  c = basis_pmatrix(a, copy = false).coeffs
  d = inv_coeff_ideals(order(a), copy = false)
  n = c[1]*d[1]
  for i = 2:degree(order(a))
    n *= c[i]*d[i]
  end
  simplify(n)
  a.norm = n
  a.has_norm = true
  return nothing
end

@doc Markdown.doc"""
    norm(a::NfRelOrdFracIdl{T, S}) -> S

Returns the norm of $a$
"""
function norm(a::NfRelOrdFracIdl{S, U, V}, copy::Type{Val{T}} = Val{true}) where {S, T, U, V}
  assure_has_norm(a)
  if copy == Val{true}
    return deepcopy(a.norm)::U
  else
    return a.norm::U
  end
end

function norm(a::NfRelOrdFracIdl, k::Union{ NfRel, AnticNumberField, NfRelNS })
  n = norm(a)
  while nf(order(n)) != k
    n = norm(n)
  end
  return n
end

function norm(a::NfRelOrdFracIdl, k::FlintRationalField)
  n = norm(a)
  while !(n isa fmpq)
    n = norm(n)
  end
  return n
end

function absolute_norm(a::NfRelOrdFracIdl)
  return norm(a, FlintQQ)
end

################################################################################
#
#  Ideal addition / GCD
#
################################################################################

@doc Markdown.doc"""
    +(a::NfRelOrdFracIdl, b::NfRelOrdFracIdl) -> NfRelOrdFracIdl

Returns $a + b$.
"""
function +(a::NfRelOrdFracIdl{T, S}, b::NfRelOrdFracIdl{T, S}) where {T, S}
  if iszero(a)
    return b
  end
  if iszero(b)
    return a
  end
  d = degree(order(a))
  H = vcat(basis_pmatrix(a), basis_pmatrix(b))
  if T != nf_elem
    H = sub(pseudo_hnf(H, :lowerleft), (d + 1):2*d, 1:d)
    return fractional_ideal(order(a), H, true)
  end
  den = lcm(denominator(a), denominator(b))
  for i = 1:d
    # We assume that the basis_pmats are lower triangular
    for j = 1:i
      H.matrix[i, j] *= den
      H.matrix[i + d, j] *= den
    end
  end
  m = simplify(den^d*(norm(a) + norm(b)))
  @assert isone(denominator(m))
  H = sub(pseudo_hnf_full_rank_with_modulus(H, numerator(m), :lowerleft), (d + 1):2*d, 1:d)
  for i = 1:d
    H.coeffs[i].den = H.coeffs[i].den*den
    H.coeffs[i] = simplify(H.coeffs[i])
  end
  return fractional_ideal(order(a), H, true)
end

+(a::NfRelOrdIdl{T, S}, b::NfRelOrdFracIdl{T, S}) where {T, S} = fractional_ideal(order(a), a) + b

+(a::NfRelOrdFracIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S} = a + fractional_ideal(order(b), b)

################################################################################
#
#  Ideal multiplication
#
################################################################################

@doc Markdown.doc"""
    *(a::NfRelOrdFracIdl, b::NfRelOrdFracIdl) -> NfRelOrdFracIdl

Returns $a \cdot b$.
"""
function *(a::NfRelOrdFracIdl{T, S, U}, b::NfRelOrdFracIdl{T, S, U}) where {T, S, U}
  if iszero(a) || iszero(b)
    return nf(order(a))()*order(a)
  end
  pba = pseudo_basis(a, copy = false)
  pbb = pseudo_basis(b, copy = false)
  ma = basis_matrix(a, copy = false)
  mb = basis_matrix(b, copy = false)
  den = denominator(a)*denominator(b)
  L = nf(order(a))
  K = base_field(L)
  d = degree(order(a))
  M = zero_matrix(K, d^2, d)
  C = Array{fractional_ideal_type(order_type(K)), 1}(undef, d^2)
  t = L()
  for i = 1:d
    for j = 1:d
      mul!(t, pba[i][1], pbb[j][1])
      T == nf_elem ? t = t*den : nothing
      elem_to_mat_row!(M, (i - 1)*d + j, t)
      C[(i - 1)*d + j] = pba[i][2]*pbb[j][2]
    end
  end
  PM = PseudoMatrix(M, C)
  PM.matrix = PM.matrix*basis_mat_inv(order(a), copy = false)
  if T != nf_elem
    H = sub(pseudo_hnf(PM, :lowerleft), (d*(d - 1) + 1):d^2, 1:d)
    return fractional_ideal(order(a), H, true)
  end
  m = simplify(den^(2*d)*norm(a)*norm(b))
  @assert isone(denominator(m))
  H = sub(pseudo_hnf_full_rank_with_modulus(PM, numerator(m), :lowerleft), (d*(d - 1) + 1):d^2, 1:d)
  for i = 1:d
    H.coeffs[i].den = H.coeffs[i].den*den
    H.coeffs[i] = simplify(H.coeffs[i])
  end
  return fractional_ideal(order(a), H, true)
end

*(a::NfRelOrdIdl{T, S}, b::NfRelOrdFracIdl{T, S}) where {T, S} = fractional_ideal(order(a), a)*b

*(a::NfRelOrdFracIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S} = a*fractional_ideal(order(b), b)

################################################################################
#
#  Division
#
################################################################################

@doc Markdown.doc"""
      divexact(a::NfRelOrdFracIdl, b::NfRelOrdFracIdl) -> NfRelOrdFracIdl
      divexact(a::NfRelOrdFracIdl, b::NfRelOrdIdl) -> NfRelOrdFracIdl
      divexact(a::NfRelOrdIdl, b::NfRelOrdFracIdl) -> NfRelOrdFracIdl

Returns $ab^{-1}$.
"""
divexact(a::NfRelOrdFracIdl{T, S}, b::NfRelOrdFracIdl{T, S}) where {T, S} = a*inv(b)

divexact(a::NfRelOrdFracIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S} = a*inv(b)

function divexact(a::NfRelOrdIdl{T, S}, b::NfRelOrdFracIdl{T, S}) where {T, S}
  O = order(a)
  return fractional_ideal(O, basis_pmatrix(a, copy = false), true)*inv(b)
end

//(a::NfRelOrdFracIdl{T, S}, b::NfRelOrdFracIdl{T, S}) where {T, S} = divexact(a, b)

//(a::NfRelOrdFracIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S} = divexact(a, b)

//(a::NfRelOrdIdl{T, S}, b::NfRelOrdFracIdl{T, S}) where {T, S} = divexact(a, b)

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

function *(a::NfRelOrdFracIdl{T, S}, b::NumFieldElem{T}) where {T, S}
  if iszero(b)
    return b*order(a)
  end
  c = b*order(a)
  return c*a
end

*(b::NumFieldElem{T}, a::NfRelOrdFracIdl{T, S}) where {T, S} = a*b

function *(a::NfRelOrdFracIdl{T, S}, b::S) where {T, S}
  pm = basis_pmatrix(a)
  pmnew = pseudo_matrix(matrix(pm), map(z -> b * z, coefficient_ideals(pm)))
  return fractional_ideal(order(a), pmnew)
end

*(a::S, b::NfRelOrdFracIdl{T, S}) where {T, S} = b * a

################################################################################
#
#  Integral ideal testing
#
################################################################################

isintegral(a::NfRelOrdFracIdl) = defines_ideal(order(a), basis_pmatrix(a, copy = false))

################################################################################
#
#  "Simplification"
#
################################################################################

# The basis_pmatrix of a NfRelOrdFracIdl should be in pseudo hnf, so it should already
# be "simple". Maybe we could simplify the coefficient ideals?
simplify(a::NfRelOrdFracIdl) = a

################################################################################
#
#  Reduction of element modulo ideal
#
################################################################################

function mod(x::S, y::T) where {S <: Union{nf_elem, NumFieldElem}, T <: Union{NfOrdFracIdl, NfRelOrdFracIdl}}
  K = parent(x)
  O = order(y)
  d = K(lcm(denominator(x, O), denominator(y)))
  dx = d*x
  dy = d*y
  if T == NfOrdFracIdl
    dy = simplify(dy)
    dynum = numerator(dy)
  else
    dynum = ideal_type(O)(O, basis_pmatrix(dy, copy = false))
  end
  dz = mod(O(dx), dynum)
  z = divexact(K(dz), d)
  return z
end

################################################################################
#
#  Inclusion of elements in ideals
#
################################################################################

@doc Markdown.doc"""
    in(x::NumFieldElem, y::NfRelOrdFracIdl)

Returns whether $x$ is contained in $y$.
"""
function in(x::NumFieldElem, y::NfRelOrdFracIdl)
  parent(x) != nf(order(y)) && error("Number field of element and ideal must be equal")
  O = order(y)
  b_pmat = basis_pmatrix(y, copy = false)
  t = zero_matrix(base_field(nf(O)), 1, degree(O))
  elem_to_mat_row!(t, 1, x)
  t = t*basis_mat_inv(O, copy = false)
  t = t*basis_mat_inv(y, copy = false)
  for i = 1:degree(O)
    if !(t[1, i] in b_pmat.coeffs[i])
      return false
    end
  end
  return true
end

################################################################################
#
#  Valuation
#
################################################################################

function valuation(A::NfRelOrdFracIdl, P::NfRelOrdIdl)
  return valuation(numerator(A), P) - valuation(denominator(A), P)
end

function valuation_naive(a::NumFieldElem, P::NfRelOrdIdl)
  @assert !iszero(a)
  #return valuation(a*order(P), P)
  d = denominator(a, order(P))
  return valuation(order(P)(d * a), P) - valuation(d, P)
end

valuation(a::NumFieldElem, P::NfRelOrdIdl) = valuation_naive(a, P)

################################################################################
#
#  Random elements
#
################################################################################

function rand(a::NfRelOrdFracIdl, B::Int)
  pb = pseudo_basis(a, copy = false)
  z = nf(order(a))()
  for i = 1:degree(order(a))
    t = rand(pb[i][2], B)
    z += t*pb[i][1]
  end
  return z
end

################################################################################
#
#  Factorization
#
################################################################################

function integral_split(a::NfRelOrdFracIdl)
  O = order(a)
  K = nf(O)
  d = inv(a + K(1)*O)
  @assert denominator(d) == 1
  n = a*d
  @assert denominator(n) == 1
  return numerator(n), numerator(d)
end

function factor(I::NfRelOrdFracIdl)
  if iszero(I)
    error("Cannot factor zero ideal")
  end
  n, d = integral_split(I)
  fn = factor(n)
  fd = factor(d)
  for (k, v) = fd
    if haskey(fn, k)
      fn[k] -= v
    else
      fn[k] = -v
    end
  end
  return fn
end

################################################################################
#
#  Hashing
#
################################################################################

function Base.hash(A::NfRelOrdFracIdl, h::UInt)
  return Base.hash(basis_pmatrix(A, copy = false), h)
end


text:
module Geodesy

using Dates
using CoordinateTransformations
using StaticArrays

import CoordinateTransformations.transform_deriv,
       CoordinateTransformations.transform_deriv_params,
       CoordinateTransformations.compose,
       CoordinateTransformations.∘

export
    # Points
    ECEF,
    ENU,
    NED,
    LLA,
    LatLon,
    UTM,
    UTMZ,

    # Other types
    Ellipsoid, ellipsoid,

    # Ellipsoids
    wgs84_ellipsoid, airy1830, clarke1866, grs80,

    # Datums
    Datum,
    WGS84, wgs84,
    OSGB36, osgb36,
    NAD27, nad27,
    NAD83, nad83,
    GRS80, grs80,
    GDA94, gda94,

    # Methods
    euclidean_distance,

    # transformation methods
    transform_deriv, transform_deriv_params, compose, ∘,
    ECEFfromLLA, LLAfromECEF, ENUfromECEF, ECEFfromENU, ENUfromLLA, LLAfromENU,
    ENUfromNED, NEDfromENU, NEDfromECEF, ECEFfromNED, NEDfromLLA, LLAfromNED,
    UTMfromLLA, LLAfromUTM, UTMfromECEF, ECEFfromUTM, ENUfromUTM, UTMfromENU, NEDfromUTM, UTMfromNED,
    UTMZfromLLA, LLAfromUTMZ, UTMZfromECEF, ECEFfromUTMZ, ENUfromUTMZ, UTMZfromENU, NEDfromUTMZ, UTMZfromNED,
    UTMZfromUTM, UTMfromUTMZ,
    WebMercatorfromLLA, LLAfromWebMercator,

    # Datum transformations
    datum_shift_ECEF,
    ITRF, GDA94,

    # UTM helpers
    utm_zone

include("ellipsoids.jl")
include("datums.jl")
include("points.jl")
include("transverse_mercator.jl")
include("polar_stereographic.jl")
include("web_mercator.jl")
include("transformations.jl")
include("conversion.jl")
include("distances.jl")
include("utm.jl")
include("datum_transformations.jl")

end # module Geodesy


text:
using HTTP
using JSON
using LightOSM
using LightXML
using Test

include("stub.jl")

@testset "LightOSM Tests" begin
    @testset "Constants" begin include("constants.jl") end
    @testset "Utilities" begin include("utilities.jl") end
    @testset "Geometry" begin include("geometry.jl") end
    @testset "Download" begin include("download.jl") end
    @testset "Nearest Node" begin include("nearest_node.jl") end
    @testset "Shortest Path" begin include("shortest_path.jl") end
end


text:
function buoyancy_flux(param_set, shf::FT, lhf, T_b, qt_b, α0_b, ts) where {FT}
    g = FT(CPP.grav(param_set))
    molmass_ratio = FT(CPP.molmass_ratio(param_set))
    lv = TD.latent_heat_vapor(param_set, T_b)
    cp_m = TD.cp_m(ts)
    return (g * α0_b / cp_m / T_b * (shf + (molmass_ratio - 1) * cp_m * T_b * lhf / lv))
end

function compute_friction_velocity(param_set, windspeed, buoy_flux, z0, z1)
    vkb = CPSGS.von_karman_const(param_set)
    logz = log(z1 / z0)
    # use neutral condition as first guess
    ustar0 = windspeed * vkb / logz
    if (abs(buoy_flux) > 1.0e-20)
        function roots(ustar)
            lmo = obukhov_length(param_set, ustar, buoy_flux)
            uf = UF.Businger(param_set, lmo)
            ζ = z1 / lmo
            ζ_0 = z0 / lmo
            Ψ_m_1 = UF.psi(uf, ζ, UF.MomentumTransport())
            Ψ_m_0 = UF.psi(uf, ζ_0, UF.MomentumTransport())
            return windspeed - ustar / vkb * (logz - Ψ_m_1 + Ψ_m_0)
        end
        sol = RS.find_zero(roots, RS.NewtonsMethodAD(ustar0), RS.CompactSolution(), RS.SolutionTolerance(1e-3), 100)
        sol.converged || error("Unconverged root solve in compute_friction_velocity")
        return sol.root
    else
        return ustar0
    end
end

"""
    exchange_coefficients_byun(param_set, Ri, zb, z0)

Ref: [Byun1990](@cite)
"""
function exchange_coefficients_byun(param_set, Ri, zb, z0)
    von_karman_const = CPSGS.von_karman_const(param_set)

    zb = zb.z
    logz = log(zb / z0)
    zfactor = zb / (zb - z0) * logz
    sb = Ri / Pr0

    if Ri > 0.0
        ζ =
            zfactor / (2.0 * beta_h * (beta_m * Ri - 1.0)) *
            ((1.0 - 2.0 * beta_h * Ri) - sqrt(1.0 + 4.0 * (beta_h - beta_m) * sb))
    else
        qb = 1.0 / 9.0 * (1.0 / (gamma_m * gamma_m) + 3.0 * gamma_h / gamma_m * sb * sb)
        pb = 1.0 / 54.0 * (-2.0 / (gamma_m * gamma_m * gamma_m) + 9.0 / gamma_m * (-gamma_h / gamma_m + 3.0) * sb * sb)
        crit = qb * qb * qb - pb * pb
        if crit < 0.0
            tb = cbrt(sqrt(-crit) + abs(pb))
            ζ = zfactor * (1.0 / (3.0 * gamma_m) - (tb + qb / tb))
        else
            angle = acos(pb / sqrt(qb * qb * qb))
            ζ = zfactor * (-2.0 * sqrt(qb) * cos(angle / 3.0) + 1.0 / (3.0 * gamma_m))
        end
    end
    lmo = zb / ζ
    ζ_0 = z0 / lmo
    uf = UF.Businger(param_set, lmo)
    Ψ_m_1 = UF.psi(uf, ζ, UF.MomentumTransport())
    Ψ_m_0 = UF.psi(uf, ζ_0, UF.MomentumTransport())
    Ψ_h_1 = UF.psi(uf, ζ, UF.HeatTransport())
    Ψ_h_0 = UF.psi(uf, ζ_0, UF.HeatTransport())

    cu = von_karman_const / (logz - Ψ_m_1 + Ψ_m_0)
    cth = von_karman_const / (logz - Ψ_h_1 + Ψ_h_0) / Pr0
    cm = cu * cu
    ch = cu * cth
    return cm, ch, lmo
end

function obukhov_length(param_set, ustar, buoy_flux)
    von_karman_const = CPSGS.von_karman_const(param_set)
    if abs(buoy_flux) > 1e-10
        lmo = -ustar^3 / buoy_flux / von_karman_const
    else
        lmo = 0.0
    end
    return lmo
end

function friction_velocity_given_windspeed(cm, windspeed)
    return sqrt(cm) * windspeed
end


text:
module DEACisf_fitness
using LinearAlgebra
const functionNames = ["mean_isf",
                       "max_isf",
                       "sum_isf",
                       "prod_isf",
                       "beta_over_two"]

function mean_isf(isfmT::Array{Float64,1},isf_model::Array{Float64,2},isf::Array{Float64,1})
    n_isf = Float64(size(isf_model,1));
    @inbounds for i = 1:size(isf_model,2)
        _s = 0.0;
        for j = 1:size(isf_model,1)
            _s += 1 + abs(1 - isf_model[j,i] / isf[j]);
        end
        isfmT[i] = _s/n_isf;
    end
    nothing
end

function max_isf(isfmT::Array{Float64,1},isf_model::Array{Float64,2},isf::Array{Float64,1})
    @inbounds for i = 1:size(isf_model,2)
        _s = 0.0;
        for j = 1:size(isf_model,1)
            _s = max(_s,abs(1 - isf_model[j,i] / isf[j]));
        end
        isfmT[i] = 1 + _s;
    end
    nothing
end

function sum_isf(isfmT::Array{Float64,1},isf_model::Array{Float64,2},isf::Array{Float64,1})
    @inbounds for i = 1:size(isf_model,2)
        _s = 1.0;
        for j = 1:size(isf_model,1)
            _s += abs(1 - isf_model[j,i] / isf[j]);
        end
        isfmT[i] = _s;
    end
    nothing
end

function prod_isf(isfmT::Array{Float64,1},isf_model::Array{Float64,2},isf::Array{Float64,1})
    @inbounds for i = 1:size(isf_model,2)
        _s = 1.0;
        for j = 1:size(isf_model,1)
            _s *= 1 + abs(1 - isf_model[j,i] / isf[j]);
        end
        isfmT[i] = _s;
    end
    nothing
end

function beta_over_two(isfmT::Array{Float64,1},isf_model::Array{Float64,2},isf::Array{Float64,1})
    @inbounds for i = 1:size(isf_model,2)
        _s = 1 + abs(1 - isf_model[end,i] / isf[end]);
        isfmT[i] = _s;
    end
    nothing
end
end


text:
#mostly legacy code used in testing

#the below funciton is now only for testing purposes
function (Θ::AbstractVolumePartsIter{TM,TV,T})(Xv::XYIterLevel,
  RegressionType::Val) where {TM<:Matrix,TV<:Vector,T}
  @unpack G,A = Θ

  #initialize A₁.==1 and all other As to the amount implied by G
  #print("prod(Θ.G): f0=$(prod(Θ.G|>Matrix)) x64=$(prod(Float64.(Θ.G |> Matrix) )) |")
  prevA₁ = A[:,1]
  factored::TM = initializeÃandfactor(Θ,Xv)

  #run the regression for A₁
  try
    A[:,1] .= reg(factored, Xv.v, RegressionType)
  catch err
    (err == InterruptException()) && throw(err)
    errmessage = "$err"
    errmessage = length(errmessage > 1000) ? errmessage[1:1000] : errmessage
    all((isfinite.(factored))) || error("Factored includes nonfinite values!
      prevA₁: $(prevA₁)
      ****************\nG: $G
      ****************\nfactored[1:1000,:]: $(factored[1:1000,:])")

    A[:,1] .= reg(factored, Xv.v, Val{:fallbackreg}())
    @warn "regression method $RegressionType failed due to $errmessage. Using fallbackreg fallback val=$(A[:,1])"
  end

  #println("A[:,1]: $(A[:,1])")

  #absμ .= vec(sum(absμₖ,dims=2))
  absμ = factored * (x->Float64(abs(x))).(A[:,1])

  return absμ
end

#primary prediction alogrithm for the demeaned version
function (Θ::AbstractVolumePartsIter{TM,TV,T})(Xv::XYIterDemean,  RegressionType::Val
    ) where {TM,TV,T}
  @unpack G,A = Θ

  #initialize A₁.==1 and all other As to the amount implied by G
  prevA₁ = A[:,1]
  factored::TM = Xv.M * initializeÃandfactor(Θ,Xv)

  #run the regression for A₁
  A[:,1] .= reg(factored, Xv.ṽ, RegressionType)
  #@eval Main demeanbeta = $(A[:,1] |> Vector)

  absμ = factored * (x->Float64(abs(x))).(A[:,1])

  return absμ
end

(Θ::VolumePartsIter)(Xv::AbstractXYIter{TM,TV,T},
    ::Val{:leveldemean}) where {TM<:CuMatrix,TV<:CuVector,T} = error(
    "Do not use CUDA for the demeaning version")

function (Θ::AbstractVolumePartsIter{TM,TV,T})(Xv::XYIterLevel, RegressionType::Val
  ) where {TM<:CuMatrix,TV<:CuVector,T}

  @unpack G,A = Θ

  #initialize A₁.==1 and all other As to the amount implied by G
  #print("prod(Θ.G): f0=$(prod(Θ.G|>Matrix)) x64=$(prod(Float64.(Θ.G |> Matrix) )) |")
  prevA₁ = A[:,1]
  factored::TM = initializeÃandfactor(Θ,Xv)

  #run the regression for A₁
  try
    A[:,1] .= reg(factored, Xv.v, RegressionType)
    #@info "culevel A[:,1]: $(A[:,1])"
    #@info "culevel factored[20,:]: $(factored[20,:])"
    #@info "culevel A[:,20]: $(A[:,20])"
  catch err #NOTE: below may be broken, could try to fix if needed
    (err == InterruptException()) && throw(err)
    errmessage = "$err"
    errmessage = length(errmessage) > 1000 ? errmessage[1:1000] : errmessage
    all((isfinite.(factored))) || @warn("Factored includes nonfinite values!
      prevA₁: $(prevA₁)
      ****************\nG: $G
      ****************\nfactored[1:1000,:]: $(factored[1:1000,:])
      attempting correction")
    factored = (x->min(x,LARGE_VAL)).(factored)

    A[:,1] .= reg(factored, Xv.v, Val{:fallbackreg}())
    @warn "regression method $RegressionType failed due to $errmessage. Using fallback val=$(A[:,1])"
  end

  #println("A[:,1]: $(A[:,1])")

  absμ = factored * (x->Float64(abs(x))).(A[:,1])

  return absμ
end

#NOTE: this is only for testing purposes, but is it even needed for that???
#handles the intercept case
function (Θ::AbstractVolumePartsIter{TM,TV,T})(Xv::XYIterLevel,
  ::Union{Val{:interceptzygote}, Val{:intercept}},

  ) where {TM<:CuMatrix,TV<:CuVector,T}

  @unpack G,A = Θ
  factored::TM = hcat(fill!(similar(Xv.v), T(1.0)), initializeÃandfactor(Θ,Xv))

  #run the regression for A₁
  A₁ = reg(factored, Xv.v, Val{:choleskyzygote}())
  absμ = factored * A₁

  return absμ
end

#not sure why I need to write this again...but it throws a method error otherwise
function (Θ::AbstractVolumePartsIter{TM,TV,T})(Xv::XYIterLevel,
  ::Union{Val{:interceptzygote},Val{:intercept}},
  ) where {TM<:Matrix,TV<:Vector,T}

  @unpack G,A = Θ
  factored::TM = hcat(fill!(similar(Xv.v), T(1.0)), initializeÃandfactor(Θ,Xv))
  #run the regression for A₁
  A₁ = reg(factored, Xv.v, Val{:choleskyzygote}())
  absμ = factored * A₁

  return absμ
end

function (Θ::VolumePartsIter)(Xv::XYIterLevel{TM,TV,T}, ::Val{:intercept}
  ) where {TM<:CuMatrix,TV<:CuVector,T}

  @unpack G,A = Θ

  #initialize A₁.==1 and all other As to the amount implied by G
  #print("prod(Θ.G): f0=$(prod(Θ.G|>Matrix)) x64=$(prod(Float64.(Θ.G |> Matrix) )) |")
  prevA₁ = A[:,1]
  factored::TM = hcat(Xv.tv1s, initializeÃandfactor(Θ,Xv))
  local b::TV
  #expectations = mean(factored, dims=1)
  #v̄ = mean(Xv.v)

  #run the regression for A₁
  try
    b = reg(factored, Xv.v, Val{:cholesky}())
    #@info "culevel A[:,1]: $(A[:,1])"
    #@info "culevel factored[20,:]: $(factored[20,:])"
    #@info "culevel A[:,20]: $(A[:,20])"
  catch err
    (err == InterruptException()) && throw(err)
    errmessage = "$err"
    errmessage = length(errmessage)  > 1000 ? errmessage[1:1000] : errmessage
    all((isfinite.(factored))) || @warn("Factored includes nonfinite values!
      prevA₁: $(prevA₁)
      ****************\nG: $G
      ****************\nfactored[1:1000,:]: $(factored[1:1000,:])
      attempting correction")
    factored .= (x->min(x,LARGE_VAL)).(factored)
    #expectations .= mean(factored, dims=1)

    b = reg(factored, Xv.v, Val{:fallbackreg}())
    @warn "regression method $RegressionType failed due to $errmessage. Using fallback val=$(A[:,1])"
  end

  #intercept = v̄ - dot(vec(expectations), A[:,1])

  #println("A[:,1]: $(A[:,1])")
  A[:,1] .= b[2:end]
  absμ = factored * b

  return absμ
end

function (Θ::AbstractVolumePartsIter{TM,TV,T})(G::Matrix{T}, Xv::XYIterLevel,
  RegressionType::Val) where {TM,TV,T}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2), Nexpanded=size(Xv.ws,1))

  #special version for zygote
  #expand a matrix for A
  (RegressionType === Val{:intercept}()) && error("Don't use the pre-compute method with an intercept regression")

  Ã = hcat(ones(T, dims.K, 1), cumprodbyrow(G))
  Ãexpanded = view(Ã, Θ.expand.index.Ã...)
  LÃexpanded = view(Ã, Θ.expand.index.LÃ...)
  factored = genabsμₖ.(Ãexpanded, LÃexpanded, Xv.ws', Xv.Rtws', Xv.RLws')' #0.4ms

  absA₁ = (x->Float64(abs(x))).(reg(factored, Xv.v, RegressionType))#0.5ms
  #println("CPU absAₜ: $absAₜ")

  #sum across signals
  absμ = factored * absA₁
  #absμ = sum(Aexpanded, dims=1) |> vec
  return absμ
end

#used in testing
function (Θ::AbstractVolumePartsIter{TM,TV,T})(G::Matrix{T}, Xv::XYIterDemean,
  RegressionType::Val) where {TM,TV,T}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2), Nexpanded=size(Xv.ws,1))

  #special version for zygote
  #expand a matrix for A
  (RegressionType === Val{:intercept}()) && error("Don't use the pre-compute method with an intercept regression")

  Ã = hcat(ones(T, dims.K, 1), cumprodbyrow(G))
  Ãexpanded = view(Ã, Θ.expand.index.Ã...)
  LÃexpanded = view(Ã, Θ.expand.index.LÃ...)
  factored = Xv.M * genabsμₖ.(Ãexpanded', LÃexpanded', Xv.ws, Xv.Rtws, Xv.RLws) #0.4ms

  absA₁ = (x->Float64(abs(x))).(reg(factored, Xv.ṽ, Val{:choleskyzygote}()))#0.5ms
  #@info ("absA₁ (:zygoteleveldemean): $(absA₁)")

  #sum across signals
  absμ = factored * absA₁
  #absμ = sum(Aexpanded, dims=1) |> vec
  return absμ
end

function (Θ::AbstractVolumePartsIter{TM,TV,T})(G::Matrix{T}, Xv::XYIterLevel,
  RegressionType::Val) where {TM<:CuMatrix,TV<:CuVector,T}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2), Nexpanded=size(Xv.ws,1))

  someones = ones(T, dims.K, 1)
  prodmatG = cumprodbyrow(G)
  Ã = hcat(someones, prodmatG)

  Ãexpanded = CuMatrix{T}(Ã[Θ.expand.index.Ã...])
  LÃexpanded = CuMatrix{T}(Ã[Θ.expand.index.LÃ...])

  factored = genabsμₖ.(Ãexpanded, LÃexpanded, Xv.ws', Xv.Rtws', Xv.RLws')' #0.4ms

  #run the regression (use a spacial non-inplace version)

  #println(sum(cholesky(factored'*factored)\(factored'*Xv.v)))
  #println("Cu absAₜT: $absAₜT")

  absA₁ = reg(factored, Xv.v, RegressionType)

  absμ = factored * absA₁

  #=return absμ=#

  return absμ
end

#verion for regressions with intercept
function (Θ::AbstractVolumePartsIter{TM,TV,T})(G::Matrix{T}, Xv::XYIterLevel,
  RegressionType::Val{:interceptzygote}) where {TM<:CuMatrix,TV<:CuVector,T}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2), Nexpanded=size(Xv.ws,1))

  someones = ones(T, dims.K, 1)
  prodmatG = cumprodbyrow(G)
  Ã = hcat(someones, prodmatG)

  Ãexpanded = CuMatrix{T}(Ã[Θ.expand.index.Ã...])
  LÃexpanded = CuMatrix{T}(Ã[Θ.expand.index.LÃ...])
  factored = CUDA.hcat(
    Xv.tv1s, genabsμₖ.(Ãexpanded, LÃexpanded, Xv.ws', Xv.Rtws', Xv.RLws')') #0.4ms

  absA₁ =reg(factored, Xv.v, Val{:choleskyzygote}())
  absμ = factored * absA₁
  #absμ = sum(Aexpanded, dims=1) |> vec
  return absμ
end
#cpu version of above
function (Θ::AbstractVolumePartsIter{TM,TV,T})(G::Matrix{T}, Xv::XYIterLevel,
  RegressionType::Val{:interceptzygote}) where {TM<:Matrix,TV<:Vector,T}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2), Nexpanded=size(Xv.ws,1))

  someones = ones(T, dims.K, 1)
  prodmatG = cumprodbyrow(G)
  Ã = hcat(someones, prodmatG)

  Ãexpanded = view(Ã, Θ.expand.index.Ã...)
  LÃexpanded = view(Ã, Θ.expand.index.LÃ...)
  factored = hcat(Xv.tv1s, genabsμₖ.(Ãexpanded, LÃexpanded, Xv.ws', Xv.Rtws', Xv.RLws')') #0.4ms

  absA₁ =reg(factored, Xv.v, Val{:choleskyzygote}())
  absμ = factored * absA₁
  #absμ = sum(Aexpanded, dims=1) |> vec
  return absμ
end

function lossforzygote(G::Matrix{T}, Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel;
  RegressionType::Val=Val{PARAM[:iterregressiontypezygote]}()) where {TM,TV,T}

  #ε = factored * Aₜ .- Xv.ṽ
  Σε2 = sum((Θ(G, Xv, RegressionType) .- Xv.v).^2) #0.3ms

  return Σε2

end

function lossforzygote(G::Matrix{T}, Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterDemean;
  RegressionType::Val=Val{PARAM[:iterregressiontypezygote]}()) where {TM,TV,T}

  #ε = factored * Aₜ .- Xv.ṽ
  Σε2 = sum((Θ(G, Xv, RegressionType) .- Xv.ṽ).^2) #0.3ms

  return Σε2
end

#computs the predicted values for each strategy k and time t and demeans the result
#via the projection matrix
function projectvolumefromA(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterDemean) where {TM,TV,T}

  @unpack A, LA = Θ.expand
  @unpack ws, Rtws, RLws, M = Xv
  #first map each coefficient to its appropriate

  absμₖ = genabsμₖ.(A', LA', ws, Rtws, RLws)

  multiply!(M, absμₖ)
  return absμₖ
end

#same as the above without demeaning
function projectvolumefromA(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel) where {TM,TV,T}

  @unpack A, LA = Θ.expand
  @unpack ws, Rtws, RLws= Xv
  #first map each coefficient to its appropriate

  absμₖ =genabsμₖ.(A', LA', ws, Rtws, RLws)

  return absμₖ
end

function projectvolumefromA(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel
    ) where {TM<:CuMatrix, TV<:CuVector, T}

  @unpack A, LA = Θ.expand
  @unpack ws, Rtws, RLws = Xv
  #first map each coefficient to its appropriate

  dims = (K=size(Θ.A,1), T = size(Θ.A,2))

  genabsμₖ.(A', LA', ws, Rtws, RLws)
  #absμₖ .= cudagenabsμₖ.(A', LA', ws, RLws)

  return absμₖ
end

##############Iter loss methods
#again, mostly legacy

function iterloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel, RegressionType::Val;
  ) where {TM<:CuMatrix, TV<:CuVector, T}

  Δ = Θ(Xv, RegressionType) .- Xv.v
  λ = CUDA.sum(Δ .* Δ)# / Ncleanpanel

  return λ
end

iterloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel, RegressionType::Val;
  ) where {TM, TV, T} = sum(
    (Θ(Xv, RegressionType) .- Xv.v).^2)

iterloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterDemean, RegressionType::Val;
  ) where {TM, TV, T} = sum(
    (Θ(Xv, RegressionType) .- Xv.ṽ).^2)

function currentloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterDemean, absμₖ::TM) where {TM, TV, T}

  someones::Vector{T} = ones(T, size(Θ.A,1))
  return sum((absμₖ*someones .- Xv.ṽ).^2)
end

function currentloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterDemean,
   absμₖ::TM) where {TM<:CuMatrix, TV<:CuVector, T}

  someones::CuVector{T} = CUDA.ones(T, size(Θ.A,1))
  return sum((absμₖ*someones .- Xv.ṽ).^2)
end

function currentloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel,
   absμₖ::TM) where {TM<:CuMatrix, TV<:CuVector, T}

  someones::CuVector{T} = CUDA.ones(T, size(Θ.A,1))
  return sum((absμₖ*someones .- Xv.v).^2)
end

function testiterloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel;
    RegressionType::Val=Val(PARAM[:iterregressiontype])) where {TM, TV, T}

  vhattest = iterloss(deepcopy(Θ), Xv, RegressionType)

  vhatver = sum((Θ(Xv, RegressionType) .- Xv.v) .^ 2)

  @assert vhatver ≈ vhattest
end

function testiterloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterDemean;
  RegressionType::Val=Val(:choleskyzygote)) where {TM, TV, T}

  vhattest = iterloss(deepcopy(Θ), Xv, RegressionType)

  vhatver = sum((Θ(Xv, RegressionType) .- Xv.ṽ) .^ 2)

  @assert vhatver ≈ vhattest
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA