text:
"""
    distance(x1,x2,gridspacing :: Number)

Compute the cartesian distance between two points given their indices and the gridspacing. It asummes uniform grid.

"""
function distance(x1,x2,gridspacing :: Number,weight=1)
    return gridspacing*hypot( x2[1]-x1[1], x2[2]-x1[2] )
end

"""
    compute_N2(xBar_Tv,z)
Take a (1,1,size(z),size(t)) profile of temperature or virtual temperature and return the Brunt - Väisälä frequency at each z level and at each t.
"""
function compute_N2(xBar_Tv,z)
    T = typeof(ustrip(g)/z[1])
    N2 = zeros(T,length(z),size(xBar_Tv,4))
    factor = ustrip(g/Dryair.cp)
    @views N2[1:end-1,:]  .= ustrip(g) .* ( (xBar_Tv[1,1,2:end,:] .- xBar_Tv[1,1,1:end-1,:]) ./(z[2:end] .- z[1:end-1]) .+ factor) ./ xBar_Tv[1,1,1:end-1,:]
    @views N2[end,:]      .= N2[end-1,:]
    ind_smallN2 = findall(abs.(N2) .< 1e-6)
    one_z = CartesianIndex(1,0)
    @inbounds for ind in ind_smallN2
        if 1 < ind[1] < length(z)
           N2[ind] =  (N2[ind + one_z] + N2[ind - one_z])/2 # If N2 is small, substite by mean of neighbours
        elseif ind[1] == 1
           N2[ind] = N2[ind + one_z]
        elseif ind[1] == length(z)
            N2[ind] = N2[ind - one_z]
        end
    end
    return N2
end

function compute_N2(xBar_Tv :: Array{ <:Quantity }, z :: Array{ <:Quantity })
    T = typeof(g/z[1])
    N2 = zeros(T,length(z),size(xBar_Tv,4))
    factor = g/Dryair.cp
    @views N2[1:end-1,:]  .= g .* ( (xBar_Tv[1,1,2:end,:] .- xBar_Tv[1,1,1:end-1,:]) ./ (z[2:end] .- z[1:end-1]) .+ factor) ./ xBar_Tv[1,1,1:end-1,:]
    @views N2[end,:]      .= N2[end-1,:]
    ind_smallN2 = findall(abs.(N2) .< 1e-6u"s^-2")
    one_z = CartesianIndex(1,0)
    @inbounds for ind in ind_smallN2
        if 1 < ind[1] < length(z)
           N2[ind] = (N2[ind + one_z] + N2[ind - one_z])/2 # If N2 is small, substite by mean of neighbours
        elseif ind[1] == 1
           N2[ind] = N2[ind + one_z]
        elseif ind[1] == length(z)
            N2[ind] = N2[ind - one_z]
        end
    end
    return N2
end

#WIP
function compute_mse(T,z,qv)
    sz = size(T)
    return  Dryair.cp*T .+ g*reshape(z,(1,1,sz[3],1)) .+ liquidwater.Lv*qv

end

function get_tendency(field :: AbstractArray{T,4}; dt = error("dt is required for the budget computation")) where T
    dfield_dt = similar(field)
    @. dfield_dt[:,:,:,1:end-1] = @views (field[:,:,:,2:end] - field[:,:,:,1:end-1]) / dt
    @. dfield_dt[:,:,:,end] = dfield_dt[:,:,:,end-1]
 return dfield_dt
end

function get_tendency(field :: AbstractArray{T,1}; dt = error("dt is required for the budget computation")) where T
    dfield_dt = similar(field)
    @. dfield_dt[1:end-1] = (field[2:end] - field[1:end-1]) / dt
    dfield_dt[end] = dfield_dt[end-1]
 return dfield_dt
end

function get_advection_asresidual(tendency,sources...)
    return reduce(+,sources) .- tendency
end

# integrate_vertically(field :: AbstractArray{T,4};dz = 1, weight = 1) where {T} = reduce(+,dz*weight.*field,dims=3)

# integrate_vertically(field :: AbstractArray{T,3};dz = 1, weight = 1) where {T} = reduce(+,dz*weight.*field,dims=2)

# integrate_vertically(field :: AbstractArray{T,2};dz = 1, weight = 1) where {T} = reduce(+,dz*weight.*field,dims=2)

# integrate_vertically(field :: AbstractArray{T,1};dz = 1 , weight = 1) where {T} = reduce(+,dz*weight*field)

# function integrate_vertically(field :: AbstractArray{T,4}; coord :: AbstractArray{T,1},weight = 1) where T
#     sz = size(field)
#     sc = size(coord)
#     integral = zeros(eltype(field),(sz[1],sz[2],1,sz[4]))
#     @inbounds for ind in CartesianIndices(field)
#         if ind[3] != sz[3]
#             integral[ind[1],ind[2],1,ind[4]] = weight*(coord[ind[3]+1] - coord[ind[3]]) * field[ind]
#         end
#     end
#     return integral
# end

# function integrate_vertically(field :: AbstractArray{T,4}; coord :: AbstractArray{T,1},weight :: AbstractArray{T,1}) where T
#     sz = size(field)
#     sc = size(coord)
#     integral = zeros(eltype(field),(sz[1],sz[2],1,sz[4]))
#     @inbounds for ind in CartesianIndices(field)
#         if ind[3] != sz[3]
#             integral[ind[1],ind[2],1,ind[4]] += weight[ind[3]]*(coord[ind[3]+1] - coord[ind[3]]) * field[ind]
#         end
#     end
#     return integral
# end

function integrate_horizontally(field :: AbstractArray{T,4}; darea) where T
    return darea*reduce(+,field)

end

function integrate_vertically(field :: AbstractArray{T,4}; coord :: AbstractArray{T,1},weight :: AbstractArray{T,4}) where T
    sz = size(field)
    sc = size(coord)
    integral = zeros(eltype(field),(sz[1],sz[2],1,sz[4]))
    for ind in CartesianIndices(field)
        if ind[3] != sz[3]
            integral[ind[1],ind[2],1,ind[4]] += weight[ind]*(coord[ind[3]+1] - coord[ind[3]]) * field[ind]
        end
    end
    return integral
end

"""
    compute_virtual_temp(temperature, specific_humidity)
Compute the virtual temperature considering only water vapor mixing ratio.
"""

function compute_virtual_temp(temperature,specific_humidity)
    return temperature*(1 + epsilon*specific_humidity)
end

function spatial_derivative!(output,field,dx,dim)
    if dim == 1
        onex = CartesianIndex((1, ntuple(_->0, ndims(u) - 1)...))

        @inbounds for ind in CartesianIndices(field)
            if ind[dim] == 1
                output[ind] = (field[ind+onex] - field[ind])/dx
            elseif ind[dim] == size(field,dim)
                output[ind] = (field[ind] - field[ind-onex])/dx
            else
                output[ind] = (field[ind+onex] - field[ind-onex])/2dx
            end
        end
    elseif dim == 2
        oney = CartesianIndex((0,1, ntuple(_->0, ndims(u) - 2)...))
        @inbounds for ind in CartesianIndices(field)
            if ind[dim] == 1
                output[ind] = (field[ind+oney] - field[ind])/dx
            elseif ind[dim] == size(field,dim)
                output[ind] = (field[ind] - field[ind-oney])/dx
            else
                output[ind] = (field[ind+oney] - field[ind-oney])/2dx
            end
        end

    end
    return output
end

function spatial_derivative(field, dx, dim)
    output = similar(field)
    return spatial_derivative!(output,field,dx,dim)
end

function get_vorticity!(output,u,v,dx,dy)
    ## This implementation assumes 4d input
    onex = CartesianIndex((1, ntuple(_->0, ndims(u) - 1)...))
    oney = CartesianIndex((0,1, ntuple(_->0, ndims(u) - 2)...))
    for ind in CartesianIndices(output)
        if ind[1] == 1
            output[ind] = (v[ind+onex] - v[ind])/dx
            elseif ind[1] == size(v,1)
            output[ind] = (v[ind] - v[ind-onex])/dx
        else
            output[ind] = (v[ind+onex] - v[ind-onex])/2dx
            end
    end
    for ind in CartesianIndices(output)
        if ind[2] == 1
            output[ind] -= (u[ind+oney] - u[ind])/dy
        elseif ind[2] == size(u,2)
            output[ind] -= (u[ind] - u[ind-oney])/dy
        else
            output[ind] -= (u[ind+oney] - u[ind-oney])/2dy
        end
    end
    return output
end

function get_divergence!(output,u,v,dx,dy)
    ## This implementation assumes 4d input
    onex = CartesianIndex((1, ntuple(_->0, ndims(u) - 1)...))
    oney = CartesianIndex((0,1, ntuple(_->0, ndims(u) - 2)...))
    for ind in CartesianIndices(output)
        if ind[1] == 1
            output[ind] = (v[ind+onex] - v[ind])/dy
            elseif ind[1] == size(v,1)
            output[ind] = (v[ind] - v[ind-onex])/dy
        else
            output[ind] = (v[ind+onex] - v[ind-onex])/2dy
            end
    end
    for ind in CartesianIndices(output)
        if ind[2] == 1
            output[ind] += (u[ind+oney] - u[ind])/dx
        elseif ind[2] == size(u,2)
            output[ind] += (u[ind] - u[ind-oney])/dx
        else
            output[ind] += (u[ind+oney] - u[ind-oney])/2dx
        end
    end
    return output
end

function get_divergence(u,v,dx,dy)
    output = similar(u)
    return get_divergence!(output,u,v,dx,dy)
end

function get_vorticity(u,v,dx,dy)
    output = similar(u)
    return get_vorticity!(output,u,v,dx,dy)
end

function get_okubo_weiss!(output,u,v,dx,dy)
    onex = CartesianIndex((1, ntuple(_->0, ndims(u) - 1)...))
    oney = CartesianIndex((0,1, ntuple(_->0, ndims(u) - 2)...))
    sz = size(u)
#    dvdx
    for ind in CartesianIndices(u)
        ##dudx and dvdx
        if ind[1] == 1
            dudx = (u[ind+onex] - u[ind])/dx
            dvdx = (v[ind+onex] - v[ind])/dx
        elseif ind[1] == sz[1]
            dudx = (u[ind] - u[ind-onex])/dx
            dvdx = (v[ind] - v[ind-onex])/dx
        else
            dudx = (u[ind+onex] - u[ind-onex])/2dx
            dvdx = (v[ind+onex] - v[ind-onex])/2dx
        end

        #dudy and dvdy
        if ind[2] == 1
            dvdy = (v[ind+oney] - v[ind])/dy
            dudy = (u[ind+oney] - u[ind])/dy
        elseif ind[2] == sz[2]
            dvdy = (v[ind] - v[ind-oney])/dy
            dudy = (u[ind] - u[ind-oney])/dy
        else
            dvdy = (v[ind+oney] - v[ind-oney])/2dy
            dudy = (u[ind+oney] - u[ind-oney])/2dy
        end
        output[ind] = dudx*dudx + dvdy*dvdy - 2dudx*dvdy + 2dvdx*dudy +
            2dvdx*dudy
    end
    return output
end

function get_okubo_weiss(u,v,dx,dy)
    ow = similar(u)
    get_okubo_weiss!(ow,u,v,dx,dy)
    return ow
end

"""
    get_saturation_vapor_pressure(T)
Receive temperature T in Kelvin and compute the saturation vapor pressure in hPa from the August-Roche-Magnus formula that approximates the solution to the Clausius-Clapeyron relationship (Wikipedia contributors. (2020, December 19). Clausius–Clapeyron relation. In Wikipedia, The Free Encyclopedia. Retrieved 06:57, December 20, 2020, from https://en.wikipedia.org/w/index.php?title=Clausius%E2%80%93Clapeyron_relation&oldid=995159175)
"""
function get_saturation_vapor_pressure(T)
    return 6.112*exp(17.67 * (T-273.15) / (243.5 + (T - 273.15)))
end

function get_saturation_vapor_pressure(T :: Quantity)
    return 6.112u"hPa"*exp(17.67 * (T-273.15u"K") / (243.5u"K" + (T - 273.15u"K")))
end

"""
    get_partial_vapor_pressure(mixing_ratio,pressure)
Receive a water vapor mixing ratio (unitless g/g) and environmental pressure and compute the partial pressure of water vapor in the same units as the input pressure.
"""
function get_partial_vapor_pressure(mixing_ratio,pressure)
    return mixing_ratio*pressure/(epsilon + mixing_ratio)
end

"""
    get_mixing_ratio(water_vapor_partial_pressure,env_pressure)
Receive a water vapor mixing ratio (unitless g/g) and environmental pressure and compute the partial pressure of water vapor in the same units as the incoming pressure.
"""
function get_mixing_ratio(water_vapor_partial_pressure,env_pressure)
    return epsilon*water_vapor_partial_pressure/(env_pressure - water_vapor_partial_pressure)
end

"""
    get_specific_entropy(temperature,mixing_ratio,pressure)
Receive temperature in Kelvin, water vapor mixing ratio (unitless g/g) and pressure (hPa) and compute the specific entropy of a parcel using equation in Emmanuel's (E94, EQN. 4.5.9)
"""
function get_specific_entropy(temperature,mixing_ratio,pressure)
    vapor_pressure = get_partial_vapor_pressure(mixing_ratio,pressure)
    saturation_vapor_pressure = get_saturation_vapor_pressure(temperature)
    RH = min(vapor_pressure/saturation_vapor_pressure,1.0)
    specific_entropy =  (Dryair.cp + mixing_ratio * Liquidwater.cp) *
        log(temperature/unit(temperature)) - Dryair.R * log((pressure - vapor_pressure)/unit(pressure)) +
        Liquidwater.Lv * mixing_ratio / temperature - mixing_ratio * Watervapor.R * log(RH)
end

"""
    get_lifted_condensation_level(temperature,relative_humidity,pressure)
Receive temperature in Kelvin, relative humidity (unitless) and pressure (hPa) and compute the lifted condensation level based on Emanuel's E94 "calcsound.f" code at http://texmex.mit.edu/pub/emanuel/BOOK/
"""
function get_lifted_condensation_level(temperature,relative_humidity,pressure)
    return pressure * (relative_humidity^(temperature/(1669.0-122.0*relative_humidity-temperature)))
end

function get_lifted_condensation_level(temperature :: Quantity ,relative_humidity :: Quantity ,pressure :: Quantity)
    return pressure * (relative_humidity^(temperature/(1669.0u"K"-122.0u"K"*relative_humidity-temperature)))
end
#we need temperature to celsius
#saturation vapor pressure

"""
    specific_humidity_to_mixing_ratio(specific_humidity)
Take a specific humidity (unitless g/g) and return a mixing ratio
"""
function specific_humidity_to_mixing_ratio(specific_humidity)
return mixing_ratio = specific_humidity / (1 - specific_humidity)

end

"""
    mixing_ratio_to_specific_humidity(mixing_ratio)
Take a mixing ratio (unitless g/g) and return a specific humidity
"""
function mixing_ratio_to_specific_humidity(mixing_ratio)
    return q = mixing_ratio / (1 + mixing_ratio)
end

"""
    get_virtual_temperature(temperature,mixing_ratio_total_water,mixing_ratio_water_vapor)
Receive temperature (K) and mixing ratios of total water and water vapor (unitless g/g) and compute the virtual temperature
"""
function get_virtual_temperature(temperature,mixing_ratio_total_water,mixing_ratio_water_vapor)
    return temperature*(1 + mixing_ratio_water_vapor/epsilon)/(1 + mixing_ratio_total_water)
end

"""
    get_buoyancy_of_lifted_parcel(tparcel,rparcel,pparcel,t,r,p,ptop=50)
"""
function get_buoyancy_of_lifted_parcel(tparcel,rparcel,pparcel,t,r,p,ptop=50u"hPa")
    n_valid_levels = findfirst(<(ptop),p)
    p = p[begin:n_valid_levels]
    t = t[begin:n_valid_levels]
    r = r[begin:n_valid_levels]
    tvirtual_diff_parcel_env = similar(t)
    parcel_sat_vapor_pressure = get_saturation_vapor_pressure(tparcel)
    parcel_vapor_pressure = get_partial_vapor_pressure(rparcel,pparcel)
    parcel_rh = min(parcel_vapor_pressure/parcel_sat_vapor_pressure  , 1.0)
    parcel_specific_entropy = get_specific_entropy(tparcel,rparcel,pparcel)
    parcel_lcl = get_lifted_condensation_level(tparcel,parcel_rh,pparcel)
    @show parcel_lcl
    below_lcl = findall(>=(parcel_lcl),p)
    above_lcl = findall(<(parcel_lcl),p)

    #These two must populate buoyancy of lifted parcel_get_vapor_pressure
    #this would be adiabatic lifting, easy enough
    for level in below_lcl
        tlifted = tparcel*(p[level]/pparcel)^(Dryair.R/Dryair.cp)
        rlifted = rparcel
        tvirtual_lifted = get_virtual_temperature(tlifted,rlifted,rlifted)
        tvirtual_env = get_virtual_temperature(t[level],r[level],r[level])
        tvirtual_diff_parcel_env[level] = tvirtual_lifted - tvirtual_env
    end

    #We start with environmental values of temperature, mixing ratio, entropy etc
    for level in above_lcl
        niter = 0
        t_previousiter = t[level]
        saturation_vapor_pressure_previousiter = get_saturation_vapor_pressure(t_previousiter)
        mixing_ratio_previousiter = get_mixing_ratio(saturation_vapor_pressure_previousiter,p[level])
        t_currentiter = 0.0u"K"
        mixing_ratio_currentiter = 0.0u"g/g"
        while (abs(t_previousiter - t_currentiter) > 0.001u"K" )
            niter += 1
            t_currentiter = t_previousiter
            saturation_vapor_pressure_currentiter = get_saturation_vapor_pressure(t_currentiter)
            mixing_ratio_currentiter = get_mixing_ratio(saturation_vapor_pressure_currentiter,p[level])
            dsdt = (Dryair.cp + rparcel*Liquidwater.cp + Liquidwater.Lv*Liquidwater.Lv*mixing_ratio_currentiter/
            (Watervapor.R*t_currentiter*t_currentiter))/t_currentiter
            vapor_pressure_currentiter = get_partial_vapor_pressure(mixing_ratio_currentiter,p[level])
            entropy_currentinter = (Dryair.cp+rparcel*Liquidwater.cp)*log(t_currentiter/unit(t_currentiter)) -
            Dryair.R*log((p[level]-vapor_pressure_currentiter)/unit(p[level])) + Liquidwater.Lv*mixing_ratio_currentiter / t_currentiter

            if niter < 3
                temperature_step = 0.3
            else
                temperature_step = 1
            end
            t_previousiter = t_currentiter + temperature_step*(parcel_specific_entropy - entropy_currentinter)/dsdt

            if (niter > 500 ) | (vapor_pressure_currentiter > ( p[level] - 1.0u"hPa") )
                error("Temperature didn't converge during lift")
            end
        end
        tvirtual_lifted = get_virtual_temperature(t_currentiter,rparcel,mixing_ratio_currentiter)
        tvirtual_env = get_virtual_temperature(t[level],r[level],r[level])
        tvirtual_diff_parcel_env[level] = tvirtual_lifted - tvirtual_env
    end
    return tvirtual_diff_parcel_env
end

"""
    get_potential_temperature(temperature, pressure, reference_pressure)
Compute potential temperature from temperature and pressure.
"""
function get_potential_temperature(temperature, pressure, reference_pressure)
    exponent = ustrip(Dryair.R / Dryair.cp)
    return temperature * (reference_pressure/pressure)^exponent
end

function get_potential_temperature(temperature :: Quantity, pressure :: Quantity, reference_pressure :: Quantity)
    exponent = Dryair.R / Dryair.cp
    return temperature * (reference_pressure/pressure)^exponent
end

"""
    get_virtual_temperature(temperature, specific_humidity)
Compute virtual temperature from temperature and specific humidity.
"""
function get_virtual_temperature(temperature, specific_humidity)
    return (one(temperature) + one(temperature)/1000*epsilon*specific_humidity)*temperature
end

function get_virtual_temperature(temperature :: Quantity, specific_humidity :: Quantity)
    return (one(temperature) + one(temperature)/1000*epsilon*specific_humidity)*temperature
end

"""
   function surface_sensible_heat_flux_to_buoyancy(SST , sensible_heat_flux ; rho = 1)
Convert surface energy fluxes in units of W/m^2 to units of buoyancy m^2/s^3).
"""
function surface_sensible_heat_flux_to_buoyancy(SST, sensible_heat_flux; rho = 1.0)
    return ustrip(g) /(ustrip(Dryair.cp)*SST) * sensible_heat_flux
end

function surface_sensible_heat_flux_to_buoyancy(SST :: Quantity, sensible_heat_flux :: Quantity; rho = 1u"kg/m^3")
    return g /(1u"kg/m^3"*Dryair.cp*SST) * sensible_heat_flux
end

"""
   function surface_latent_heat_flux_to_buoyancy(SST , sensible_heat_flux ; rho = 1)
Convert surface energy fluxes in units of W/m^2 to units of buoyancy m^2/s^3).
"""
function surface_latent_heat_flux_to_buoyancy(SST, latent_heat_flux; rho = 1.0)
    return ustrip(g)/(one(SST)*ustrip(Dryair.cp)*SST)*(epsilon*ustrip(Dryair.cp)*SST/ustrip(Liquidwater.Lv)*latent_heat_flux)
end

function surface_latent_heat_flux_to_buoyancy(SST :: Quantity, latent_heat_flux :: Quantity; rho = 1u"kg/m^3")
    return g/(1u"kg/m^3"*Dryair.cp*SST)*(epsilon*Dryair.cp*SST/Liquidwater.Lv*latent_heat_flux)
end

"""
   get_buoyancy(temperature_anomaly,mean_temperature)
Compute buoyancy from a temperature anomaly and a temperature profile as in RamirezReyes and Yang 2021
"""
function get_buoyancy(temperature_anomaly,mean_temperature)
    return ustrip(g) * ustrip(temperature_anomaly)/ustrip(mean_temperature)
end

function get_buoyancy(temperature_anomaly :: Quantity ,mean_temperature :: Quantity)
    return g * temperature_anomaly/mean_temperature
end

"""
   radiative_heating_rate_to_buoyancy(mean_temperature,radiative_heating_rate)
Convert radiative heating rate in K/s to buoyancy production rate in m/s^2 using the a mean sounding as reference profile
"""
function radiative_heating_rate_to_buoyancy(mean_temperature,radiative_heating_rate)
    return ustrip(g) * radiative_heating_rate / mean_temperature
end

function radiative_heating_rate_to_buoyancy(mean_temperature :: Quantity,radiative_heating_rate :: Quantity)
    return g * radiative_heating_rate / mean_temperature
end

"""
   get_density_ideal_gas(pressure, temperature)
Compute density of dry air from pressure and temperature using the Ideal gas law.
"""
function get_density_ideal_gas(pressure, temperature)
    return pressure / temperature / eltype(temperature)(ustrip(Dryair.R))
end

function get_density_ideal_gas(pressure :: Quantity, temperature :: Quantity)
    return pressure / temperature / Dryair.R
end


text:
function eigsolve(A, x₀, howmany::Int, which::Selector, alg::Lanczos)
    krylovdim = alg.krylovdim
    maxiter = alg.maxiter
    howmany > krylovdim && error("krylov dimension $(krylovdim) too small to compute $howmany eigenvalues")

    ## FIRST ITERATION: setting up
    numiter = 1
    # Compute Lanczos factorization
    iter = LanczosIterator(A, x₀, alg.orth)
    fact = initialize(iter; verbosity = alg.verbosity - 2)
    numops = 1
    sizehint!(fact, krylovdim)
    β = normres(fact)
    tol::eltype(β) = alg.tol
    if normres(fact) > tol || howmany > 1
        while length(fact) < krylovdim
            fact = expand!(iter, fact; verbosity = alg.verbosity-2)
            numops += 1
            normres(fact) <= tol && length(fact) >= howmany && break
        end
    end

    # Process
    # allocate storage
    HH = fill(zero(eltype(fact)), krylovdim+1, krylovdim)
    UU = fill(zero(eltype(fact)), krylovdim, krylovdim)

    # initialize
    β = normres(fact)
    m = length(fact)
    U = copyto!(view(UU, 1:m, 1:m), I)
    f = view(HH, m+1, 1:m)
    T = rayleighquotient(fact) # symtridiagonal

    # compute eigenvalues
    D, U = tridiageigh!(T, U)
    by, rev = eigsort(which)
    p = sortperm(D, by = by, rev = rev)
    D, U = permuteeig!(D, U, p)
    mul!(f, view(U, m, :), β)
    converged = 0
    while converged < length(fact) && abs(f[converged+1]) < tol
        converged += 1
    end

    if alg.verbosity > 1
        msg = "Lanczos eigsolve in iter $numiter: "
        msg *= "$converged values converged, normres = ("
        msg *= @sprintf("%.2e", abs(f[1]))
        for i = 2:howmany
            msg *= ", "
            msg *= @sprintf("%.2e", abs(f[i]))
        end
        msg *= ")"
        @info msg
    end

    ## OTHER ITERATIONS: recycle
    while numiter < maxiter && converged < howmany
        numiter += 1

        # Determine how many to keep
        keep = div(3*krylovdim + 2*converged, 5) # strictly smaller than krylovdim since converged < howmany <= krylovdim, at least equal to converged

        # Restore Lanczos form in the first keep columns
        H = fill!(view(HH, 1:keep+1, 1:keep), 0)
        @inbounds for j = 1:keep
            H[j,j] = D[j]
            H[keep+1,j] = f[j]
        end
        @inbounds for j = keep:-1:1
            h, ν = householder(H, j+1, 1:j, j)
            H[j+1,j] = ν
            H[j+1,1:j-1] .= 0
            lmul!(h, H)
            rmul!(view(H, 1:j, :), h')
            rmul!(U, h')
        end
        @inbounds for j = 1:keep
            fact.αs[j] = H[j,j]
            fact.βs[j] = H[j+1,j]
        end

        # Update B by applying U using Householder reflections
        B = basis(fact)
        basistransform!(B, view(U, :, 1:keep))
        # for j = 1:m
        #     h, ν = householder(U, j:m, j)
        #     lmul!(h, view(U, :, j+1:krylovdim))
        #     rmul!(B, h')
        # end
        r = residual(fact)
        B[keep+1] = rmul!(r, 1/β)

        # Shrink Lanczos factorization
        fact = shrink!(fact, keep)

        # Lanczos factorization: recylce fact
        while length(fact) < krylovdim
            fact = expand!(iter, fact; verbosity = alg.verbosity-2)
            numops += 1
            normres(fact) < tol && length(fact) >= howmany && break
        end

        # post process
        β = normres(fact)
        m = length(fact)
        U = copyto!(view(UU, 1:m, 1:m), I)
        f = view(HH, m+1, 1:m)
        T = rayleighquotient(fact) # symtridiagonal

        # compute eigenvalues
        D, U = tridiageigh!(T, U)
        by, rev = eigsort(which)
        p = sortperm(D, by = by, rev = rev)
        D, U = permuteeig!(D, U, p)
        mul!(f, view(U,m,:), β)
        converged = 0
        while converged < length(fact) && abs(f[converged+1]) < tol
            converged += 1
        end

        if alg.verbosity > 1
            msg = "Lanczos eigsolve in iter $numiter: "
            msg *= "$converged values converged, normres = ("
            msg *= @sprintf("%.2e", abs(f[1]))
            for i = 2:howmany
                msg *= ", "
                msg *= @sprintf("%.2e", abs(f[i]))
            end
            msg *= ")"
            @info msg
        end
    end

    if converged > howmany
        howmany = converged
    end
    values = D[1:howmany]

    # Compute eigenvectors
    V = view(U, :, 1:howmany)

    # Compute convergence information
    vectors = let B = basis(fact)
        [B*v for v in cols(V)]
    end
    residuals = let r = residual(fact)
        [r*last(v) for v in cols(V)]
    end
    normresiduals = let f = f
        map(i->abs(f[i]), 1:howmany)
    end

    if alg.verbosity > 0
        if converged < howmany
            @warn """Lanczos eigsolve finished without convergence after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        else
            @info """Lanczos eigsolve finished after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        end
    end

    return values, vectors, ConvergenceInfo(converged, residuals, normresiduals, numiter, numops)
end

ext: jl


text:
module Glu_ED_Model
using JuMP
using KNITRO

function glu_ed_model(glc_ext_input)
    #glc_ext_input = 0.05
    glu_ed = Model(
        optimizer_with_attributes(
            KNITRO.Optimizer,
            "ms_enable" => 1,
            "ms_maxsolves" => 200,
        ),
    )
    JuMP.set_silent(glu_ed)

    # Thermodynamic function
    thermo_factor(x) = tanh(-0.7 * x) # x = ΔrG/RT
    register(glu_ed, :thermo_factor, 1, thermo_factor; autodiff = true)

    # Bounds
    LB_enz = 0.0      # [g enz / g DW cell]
    UB_enz = 1.0      # [g enz / g DW cell]
    LB_met = log(1e-6)   # [log(1 uM)]
    UB_met = log(100e-3) # [log(100 mM)]
    LB_v = -50        # [mmol/gDW/h]
    UB_v = 50        # [mmol/gDW/h]
    LB_dg = -100        # [kJ/mol]
    UB_dg = 100        # [kJ/mol]

    @variables glu_ed begin
        # biomass [1/h] (g6p + ATP -> ADP)
        0 <= mu <= 10

        # Enzymes concentrations [g enz / g DW]
        LB_enz <= pts <= UB_enz # enzyme for: glc external + ATP -> ADP + g6p
        LB_enz <= ed <= UB_enz # enzyme for: g6p + NAD + ADP + NADP -> NADH + ATP + NADPH + pep + pyr
        LB_enz <= pyk <= UB_enz # enzyme for: pep + ADP -> ATP + pyr
        LB_enz <= ldh <= UB_enz # enzyme for: pyr + NADH -> NAD + lac
        LB_enz <= ppc <= UB_enz # enzyme for: pep + CO2 -> oac + P
        LB_enz <= akgsyn <= UB_enz # enzyme for: pyr + oac -> akg
        LB_enz <= gdhm <= UB_enz # enzyme for: akg + NH3 + NADPH -> glu + NADP
        LB_enz <= burn <= UB_enz # enzyme for: ATP -> ADP
        LB_enz <= nadtrdh <= UB_enz # enzyme for: NAD + NADPH -> NADH + NADP
        LB_enz <= lp <= UB_enz # enzyme for: lac internal -> lac external
        LB_enz <= glnsyn <= UB_enz # enzyme for: glu + ATP + NH3 -> gln + ADP
        LB_enz <= sink <= UB_enz # enzyme for: gln ->

        # Metabolite concentrations [log(M)]
        LB_met <= g6p <= UB_met # glucose-6-phosphate
        LB_met <= pyr <= UB_met # pyruvate
        LB_met <= pep <= UB_met # phosphoenolpyrovate
        LB_met <= oac <= UB_met # oxalacetate
        LB_met <= akg <= UB_met # alpha-ketogluterate
        LB_met <= glu <= UB_met # glutamate
        LB_met <= gln <= UB_met # glutamine
        LB_met <= lac <= UB_met # lactate
        LB_met <= atp <= UB_met # ATP
        LB_met <= adp <= UB_met # ADP
        LB_met <= nad <= UB_met # NAD
        LB_met <= nadh <= UB_met # NADH
        LB_met <= nadp <= UB_met # NADP
        LB_met <= nadph <= UB_met # NADPH
        LB_met <= co2 <= UB_met # CO2
        LB_met <= phos <= UB_met # phosphate
        LB_met <= nh3 <= UB_met # NH3

        # Fluxes [mmol/gDW/h]
        LB_v <= v_pts <= UB_v
        LB_v <= v_ed <= UB_v
        LB_v <= v_pyk <= UB_v
        LB_v <= v_ldh <= UB_v
        LB_v <= v_ppc <= UB_v
        LB_v <= v_akgsyn <= UB_v
        LB_v <= v_gdhm <= UB_v
        LB_v <= v_burn <= UB_v
        LB_v <= v_nadtrdh <= UB_v
        LB_v <= v_lp <= UB_v
        LB_v <= v_glnsyn <= UB_v
        LB_v <= v_sink <= UB_v

        # Thermodynamic variables
        LB_dg <= dg_pts <= UB_dg
        LB_dg <= dg_ed <= UB_dg
        LB_dg <= dg_pyk <= UB_dg
        LB_dg <= dg_ldh <= UB_dg
        LB_dg <= dg_ppc <= UB_dg
        LB_dg <= dg_akgsyn <= UB_dg
        LB_dg <= dg_gdhm <= UB_dg
        LB_dg <= dg_burn <= UB_dg
        LB_dg <= dg_nadtrdh <= UB_dg
        LB_dg <= dg_lp <= UB_dg
        LB_dg <= dg_glnsyn <= UB_dg
        LB_dg <= dg_sink <= UB_dg
    end

    @NLparameters glu_ed begin
        RT == 8.3145e-3 * 298.15 # [kJ/mol]

        # Enzyme rates
        kcat_pts == 213.75
        kcat_ed == 16.0100104283114
        kcat_pyk == 117.269621287688
        kcat_ldh == 8.65643947924285
        kcat_ppc == 46
        kcat_akgsyn == 4
        kcat_gdhm == 342
        kcat_burn == 22
        kcat_nadtrdh == 51.9313402919591
        kcat_lp == 100
        kcat_glnsyn == 33
        kcat_sink == 50 # arbitrary for now

        # Gibbs energy of reaction
        dG0_pts == -16.7
        dG0_ed == -4.71
        dG0_pyk == -31.7
        dG0_ldh == -23.7
        dG0_ppc == -40.3
        dG0_akgsyn == -60.7
        dG0_gdhm == -33.4
        dG0_burn == -29.6
        dG0_nadtrdh == 0
        dG0_lp == -10.3
        dG0_glnsyn == -15.3
        dG0_sink == -10   # arbitrary for now

        # media conditions
        glc_ext == log(glc_ext_input) # [log(50 mM)]
        lac_ext == log(10e-6) # [log(10 uM)]

        # capacity constraint
        total_proteome_mass_fraction == 0.26 # [g enz/gDW]

        # maintenance requirement
        min_burn_flux == 1 # [mmol/gDW/h]
        min_sink_flux == 1
    end

    @NLconstraints glu_ed begin
        # ΔG at concentrations
        dg_pts == dG0_pts + RT * (adp + g6p - atp - glc_ext)
        dg_ed == dG0_ed + RT * (nadph + nadh + atp + pep + pyr - adp - nadp - nad - g6p)
        dg_pyk == dG0_pyk + RT * (pyr + atp - adp - pep)
        dg_ldh == dG0_ldh + RT * (lac + nad - nadh - pyr)
        dg_ppc == dG0_ppc + RT * (pep + co2 - oac - phos)
        dg_akgsyn == dG0_akgsyn + RT * (pyr + oac - akg - 2 * co2)
        dg_gdhm == dG0_gdhm + RT * (akg + nadh + nh3 - glu - nad)
        dg_burn == dG0_burn + RT * (adp - atp)
        dg_nadtrdh == dG0_nadtrdh + RT * (nad + nadph - nadh - nadp)
        dg_lp == dG0_lp + RT * (lac_ext - lac)
        dg_glnsyn == dG0_glnsyn + RT * (gln + adp - glu - atp - nh3)
        dg_sink == dG0_sink + RT * (-gln)

        # flux bounds due to kinetics and thermodynamics
        v_pts == pts * kcat_pts * thermo_factor(dg_pts / RT)
        v_ed == ed * kcat_ed * thermo_factor(dg_ed / RT)
        v_pyk == pyk * kcat_pyk * thermo_factor(dg_pyk / RT)
        v_ldh == ldh * kcat_ldh * thermo_factor(dg_ldh / RT)
        v_ppc == ppc * kcat_ppc * thermo_factor(dg_ppc / RT)
        v_akgsyn == akgsyn * kcat_akgsyn * thermo_factor(dg_akgsyn / RT)
        v_gdhm == gdhm * kcat_gdhm * thermo_factor(dg_gdhm / RT)
        v_burn == burn * kcat_burn * thermo_factor(dg_burn / RT)
        v_nadtrdh == nadtrdh * kcat_nadtrdh * thermo_factor(dg_nadtrdh / RT)
        v_lp == lp * kcat_lp * thermo_factor(dg_lp / RT)
        v_glnsyn == glnsyn * kcat_glnsyn * thermo_factor(dg_glnsyn / RT)
        v_sink == sink * kcat_sink * thermo_factor(dg_sink / RT)

        # mass balance constraints
        v_pts - v_ed == mu # g6p
        v_ed - v_pyk - v_ppc == 0  # pep
        v_ed + v_pyk - v_ldh - v_akgsyn == 0  # pyr
        v_ldh - v_lp == 0  # lac
        v_ppc - v_akgsyn == 0  # oac
        v_akgsyn - v_gdhm == 0  # akg
        v_gdhm - v_glnsyn == 0  # glu
        v_glnsyn - v_sink == 0  # gln
        v_ed + v_pyk - v_pts - v_burn - v_glnsyn == mu # atp
        v_pts + v_burn + v_glnsyn - v_ed - v_pyk == -mu # adp
        v_ldh - v_ed - v_nadtrdh - v_akgsyn == 0  # nad
        v_ed + v_nadtrdh + v_akgsyn - v_ldh == 0  # nadh
        v_gdhm + v_nadtrdh - v_akgsyn - v_ed == 0  # nadp
        v_ed + v_akgsyn - v_nadtrdh - v_gdhm == 0  # nadph

        # density constraint(s) (can add more, e.g. membrane)
        pts + ed + pyk + ldh + ppc + akgsyn + gdhm + burn + nadtrdh + lp + glnsyn + sink <=
        total_proteome_mass_fraction

        # minimum maintenance
        min_burn_flux <= v_burn
        min_sink_flux <= v_sink

        # measured ratio constraints (not strictly necessary)
        # atp   == log(10)  + adp  # atp/adp = 10, but remember that the concentration variables are logged
        # nadh  == log(0.2) + nad
        # nadph == log(0.2) + nadp
    end

    @objective(glu_ed, Max, gln)
    optimize!(glu_ed)
    objective_value(glu_ed)

    results_ed = Dict(
        "mu" => value(mu),
        "pts" => value(pts),
        "ed" => value(ed),
        "pyk" => value(pyk),
        "ldh" => value(ldh),
        "ppc" => value(ppc),
        "akgsyn" => value(akgsyn),
        "gdhm" => value(gdhm),
        "burn" => value(burn),
        "nadtrdh" => value(nadtrdh),
        "lp" => value(lp),
        "glc_ext" => value(glc_ext),
        "g6p" => value(g6p),
        "pyr" => value(pyr),
        "pep" => value(pep),
        "oac" => value(oac),
        "akg" => value(akg),
        "lac" => value(lac),
        "glu" => value(glu),
        "atp" => value(atp),
        "adp" => value(adp),
        "nad" => value(nad),
        "nadh" => value(nadh),
        "nadp" => value(nadp),
        "nadph" => value(nadph),
        "v_pts" => value(v_pts),
        "v_ed" => value(v_ed),
        "v_pyk" => value(v_pyk),
        "v_ldh" => value(v_ldh),
        "v_ppc" => value(v_ppc),
        "v_akgsyn" => value(v_akgsyn),
        "v_gdhm" => value(v_gdhm),
        "v_burn" => value(v_burn),
        "v_nadtrdh" => value(v_nadtrdh),
        "v_lp" => value(v_lp),
        "v_glnsyn" => value(v_glnsyn),
        "v_sink" => value(v_sink),
    )
    return results_ed
end

end #module


text:
function PTCKrylovinit(
    x0,
    dx,
    F!,
    Jvec,
    delta0,
    Pvec,
    PvecKnowsdelta,
    pside,
    lsolver,
    eta,
    fixedeta,
    lmaxit,
    maxit,
    printerr,
    pdata,
)
    #
    #   Initialize the PTC-Krylov iteration.
    #
    n = length(x0)
    x = copy(x0)
#    tmp1 = zeros(n)
#    tmp2 = zeros(n)
#    tmp3 = zeros(n)
#    tmp4 = zeros(n)
#    kl_store = (tmp1, tmp2, tmp3, tmp4)
    kl_store = kstore(n, lsolver)
    ItRules = (
        dx = dx,
        f = F!,
        Jvec = Jvec,
        delta0 = delta0,
        Pvec = Pvec,
        PvecKnowsdelta = PvecKnowsdelta,
        pside = pside,
        lsolver = lsolver,
        kl_store = kl_store,
        eta = eta,
        fixedeta = fixedeta,
        lmaxit = lmaxit,
        maxit = maxit,
        printerr = printerr,
        pdata = pdata,
    )
    return (ItRules, x, n)
end

"""
PTCUpdatei(FPS::AbstractArray, FS, x, ItRules, step, residm, delta)

Updates the PTC-Krylov iteration. This is a much simpler algorithm
than Newton-Krylov. In particular, there is no line search to manage.

Do not mess with this function!

"""
#function PTCUpdatei(FPS::AbstractArray, FS, x, ItRules, step, residm, delta, etag)
function PTCUpdatei(FPS, FS, x, ItRules, step, residm, delta, etag)
    T = eltype(FPS)
    F! = ItRules.f
    pdata = ItRules.pdata
    #
    #
    #    step .= -(FPF \ FS)
    step .*= 0.0
    #
    #   If the preconditioner can use delta, tell it what delta is.
    #
    PvecKnowsdelta = ItRules.PvecKnowsdelta
    delta2pdata(PvecKnowsdelta, delta, pdata)
    #
    kout = Krylov_Step!(step, x, FS, FPS, ItRules, etag, delta)
    Lstats = kout.Lstats
    step = kout.step
    #
    # update solution/function value
    #
    x .= x + step
    EvalF!(F!, FS, x, pdata)
    resnorm = norm(FS)
    #
    # Update delta
    #
    delta *= (residm / resnorm)
    return (x, delta, FS, resnorm, Lstats)
end

"""
delta2pdata(PvecKnowsdelta, delta, pdata)

If your preconditioner is aware of the pseuto time step (delta)
put the value where it's supposed to be inside the precomputed data.

I also check that this has been done right and complain if not.
"""
function delta2pdata(PvecKnowsdelta, delta, pdata)
    PvecKnowsdelta || return
    # Once you're here you've told me that the preconditioner is delta-aware.
    # I will look for the array deltaval before I write to it.
    T = typeof(pdata)
    Pnames = fieldnames(T)
    valok = false
    for ip in Pnames
        valok = valok || :deltaval == ip
    end
    valok ? (pdata.deltaval[1] = delta) :
    error(
        "PvecKnowsdelta is set to true, but there the array
         deltaval is not a field of pdata. Check the docstrings.",
    )
end


text:
function assemPotentialKandVwithDuffy(BEMesh,quadRuleA,quadRuleB)

	#constant
	mu0=4*pi*10^-7; mu0inv=1/mu0

	# dimension
	n=size(BEMesh.edges,1)
	n2=n*n

	#matrix to be computed in sparse format
	VE=zeros(n2); VI=zeros(n2); VJ=zeros(n2) #P0-P0, one contribution by pair of edges

	KE1=zeros(n,2); KI1=zeros(n,2); KJ1=zeros(n,2) #P1-P0, two contributions by edge
	KE2=zeros(n2,2); KI2=zeros(n2,2); KJ2=zeros(n2,2) #P1-P0, two contributions by pair of edges
	#reshape at the end vector size(2*n+2*n2)

	#quadrature rules
	xGaussA=quadRuleA.x
	wGaussA=quadRuleA.w
	nQuadA=size(quadRuleA.w,1)
	xGaussB=quadRuleB.x
	wGaussB=quadRuleB.w
	nQuadB=size(quadRuleB.w,1)

	## single integral
	e=collect(1:n)
	edgeToNode=BEMesh.edges
	h=BEMesh.length
	for i in 1:nQuadA
		KE1=KE1+[0.5*h*0.5*(wGaussA[i]*(1-(0.5+xGaussA[i]*0.5)))  0.5*h*0.5*(wGaussA[i]*( (0.5+xGaussA[i]*0.5))) ]
		#The famous 1/2 for the angular coefficient is needed
		#However the BEM solution evaluated at a point not on the boundary but close to it is wrong, I guess because this anglar coefficient
		# which should be 1 must be wrong if we are too close ...
		#without the 1/2:
		#KE1=KE1+[h*0.5*(wGaussA[i]*(1-(0.5+xGaussA[i]*0.5)))  h*0.5*(wGaussA[i]*( (0.5+xGaussA[i]*0.5))) ]
	end

	KI1=[e e]
	KJ1=edgeToNode

	## double integrals
	# pairing of each edge with each other edge
	#ind=collect(1:n) #collect not needed
	ind=1:n
	e1=repeat(ind,1,n)
	#e1=reshape(transpose(e1),n*n,1) # wrong, this gives a (n^2)x1 matrix
									# contrary to matlab a Nx1 matrix is not the same as a vector of size N !!!
									# and then a dim is added in when using it as an index in edges
									# we need a vector of length n^2:
	e1=reshape(transpose(e1),n*n)

	#e2=repeat(ind,n,1)	# this gives a n^2x1 matrix and then a dim is added in when using it as index
						# we need a vector of length n^2:
	e2=repeat(ind,n)
	#@show e1,e2

	#the two nodes of each edge
	edgeToNode1=BEMesh.edges[e1,:]
	edgeToNode2=BEMesh.edges[e2,:]
	#@show edgeToNode1
	#@show size(edgeToNode1)

	#in order to remap the quadrature we need the midpoints
	#and length of first and second edge
    PA=0.5*(BEMesh.coordinates[edgeToNode1[:,1],:]+BEMesh.coordinates[edgeToNode1[:,2],:])
    PB=0.5*(BEMesh.coordinates[edgeToNode2[:,1],:]+BEMesh.coordinates[edgeToNode2[:,2],:])
	#@show PA

	hA=BEMesh.length[e1]
	hB=BEMesh.length[e2]
	#@show hA,hB

	#quadrature points on first and second edge
	PAG=zeros(size(PA,1),2,nQuadA)
	PBG=zeros(size(PB,1),2,nQuadB)
	for i in 1:nQuadA
		PAG[:,:,i]=PA + xGaussA[i]*0.5*(BEMesh.coordinates[edgeToNode1[:,2],:]-BEMesh.coordinates[edgeToNode1[:,1],:])
	end
	for i in 1:nQuadB
		PBG[:,:,i]=PB + xGaussB[i]*0.5*(BEMesh.coordinates[edgeToNode2[:,2],:]-BEMesh.coordinates[edgeToNode2[:,1],:])
	end

	#@show size(PAG)
	#@show PAG
	#@show size(PBG)
	#@show PBG

	#normal on the second edge
	normal2=BEMesh.normalsIn[e2,:]

	# pairing of each quadrature point of the first edge with each of the second edge
	idA=reshape(transpose(repeat(1:nQuadA,1,nQuadB)),nQuadA*nQuadB)
	idB=repeat(1:nQuadB,nQuadA)

	#@show idA
	#@show idB

	#@show size(e1,1)
	#@show size(e2)

	#for e in 1:size(e1,1) # ==size(e2,1)
	for e in 1:n2 #loop over all pairings of edges

		test=e1[e]-e2[e]

		if(test != 0) #not the same 2 edges, we use classic quadrature
			for i in 1:(nQuadA*nQuadB)
				#Evaluate Green function and first derivative (in y=B)
				G = Green(PAG[e,1,idA[i]],PAG[e,2,idA[i]],PBG[e,1,idB[i]],PBG[e,2,idB[i]],derivative=1)
				#@show G
				#@show G.G_xy
				#Green function value
				G_P1GP2G = G.G_xy
				#@show G_P1GP2G

				#Neumann trace
				dnG_P1GP2G = G.dyrG_xy * normal2[e,1] +  G.dyzG_xy * normal2[e,2]
				#@show dnG_P1GP2G

				VE[e]=VE[e]+wGaussA[idA[i]]*wGaussB[idB[i]]*0.5*0.5*hA[e]*hB[e]*G_P1GP2G

				#@show VE[e]
				#@show size(KE2)
				#@show size(KE2[e,:])

				KE2[e:e,:]=KE2[e:e,:]+mu0inv*0.25*wGaussA[idA[i]]*wGaussB[idB[i]]*
					[ hA[e]*hB[e]*(1.0./PBG[e,1,idB[i]]).*dnG_P1GP2G*(1-(0.5+xGaussB[idB[i]]*0.5))  hA[e]*hB[e]*(1.0./PBG[e,1,idB[i]]).*dnG_P1GP2G*((0.5+xGaussB[idB[i]]*0.5)) ]

				#@show KE2[e,:]
				# strange in julia: KE2[e,:] is a column vector !!!! use KE2[e:e,:] to get the natural row vector
			end
		else #test == 0, same 2 edges, we need to do something for the singularity. Let's try Duffy trick
			ee=e1[e]
			edgeToNode=BEMesh.edges[ee,:]
			#@show edgeToNode
			#@show edgeToNode[1]
			#@show edgeToNode[2]
			#readline()

			Pmid=0.5*(BEMesh.coordinates[edgeToNode[1],:]+BEMesh.coordinates[edgeToNode[2],:])
			Pvec=0.5*(BEMesh.coordinates[edgeToNode[1],:]-BEMesh.coordinates[edgeToNode[2],:])
			h=BEMesh.length[ee]

			normal22=BEMesh.normalsIn[ee,:]

			##### first version not vectorized
			#for i in 1:nQuadA
			#	for j in 1:nQuadB
			#		PaG=Pmid + xGaussA[i]*Pvec
			#		PbG=Pmid + xGaussB[j]*Pvec

			#		u=0.5*(xGaussA[i]+1)*(xGaussB[j]+1)-1
			#		PuG=Pmid + u*Pvec

			#		Gau = Green(PaG[1],PaG[2],PuG[1],PuG[2],derivative=1)
			#		Gub = Green(PuG[1],PuG[2],PbG[1],PbG[2],derivative=1)

			#		G_au=Gau.G_xy
			#		G_ub=Gub.G_xy

			#		VE[e]=VE[e]+wGaussA[i]*wGaussB[j]*0.5*0.5*h*h*(0.5*(xGaussA[i]+1)*G_au + 0.5*(xGaussB[j]+1)*G_ub)

			#		dnG_au = Gau.dyrG_xy * normal22[1] +  Gau.dyzG_xy * normal22[2]
			#		dnG_ub = Gub.dyrG_xy * normal22[1] +  Gub.dyzG_xy * normal22[2]

			#		tmp1 = ( 0.5*(xGaussA[i]+1)*(1.0/PuG[1])*dnG_au*(1-(0.5+u*0.5)) +  0.5*(xGaussB[j]+1)*(1.0/PbG[1])*dnG_ub*(1-(0.5+xGaussB[j]*0.5)) )
			#		tmp2 = ( 0.5*(xGaussA[i]+1)*(1.0/PuG[1])*dnG_au*((0.5+u*0.5)) + 0.5*(xGaussB[j]+1)*(1.0/PbG[1])*dnG_ub*((0.5+xGaussB[j]*0.5)) )
			#		KE2[e:e,:]=KE2[e:e,:]+mu0inv*0.25*wGaussA[i]*wGaussB[j]*h*h*[tmp1  tmp2]

			#		#@show Gau
			#		#readline()

			#	end
			#end

			###### vectorized version
			PaG = (Pmid .+ (xGaussA[idA]'.*Pvec))'
			PbG = (Pmid .+ (xGaussB[idB]'.*Pvec))'
			u=0.5*(xGaussA[idA].+1).*(xGaussB[idB].+1).-1
			PuG = (Pmid .+ (u'.*Pvec))'

			#@show size(PaG)

			Gau = @.Green(PaG[:,1],PaG[:,2],PuG[:,1],PuG[:,2],derivative=1)
			Gub = @.Green(PuG[:,1],PuG[:,2],PbG[:,1],PbG[:,2],derivative=1)

			G_au=getindex.(Gau,:G_xy)
			G_ub=getindex.(Gub,:G_xy)

			#@show size(G_au)

			VE[e]=VE[e]+0.5*0.5*h*h*sum(wGaussA[idA].*wGaussB[idB].*(0.5*(xGaussA[idA].+1).*G_au .+ 0.5*(xGaussB[idB].+1).*G_ub))
			#@show VE[e]
			#readline()

			dnG_au = getindex.(Gau,:dyrG_xy) * normal22[1] +  getindex.(Gau,:dyzG_xy) * normal22[2]
			dnG_ub = getindex.(Gub,:dyrG_xy) * normal22[1] +  getindex.(Gub,:dyzG_xy) * normal22[2]

			tmp1 = mu0inv*0.25*h*h*wGaussA[idA].*wGaussB[idB].*( 0.5*(xGaussA[idA].+1).*(1.0./PuG[:,1]).*dnG_au.*(1.0 .- (0.5.+u*0.5)) +  0.5*(xGaussB[idB].+1).*(1.0./PbG[:,1]).*dnG_ub.*(1.0 .- (0.5.+xGaussB[idB]*0.5)) )
			#@show size(tmp1)
			tmp2 = mu0inv*0.25*h*h*wGaussA[idA].*wGaussB[idB].*( 0.5*(xGaussA[idA].+1).*(1.0./PuG[:,1]).*dnG_au.*(0.5 .+ u*0.5)  + 0.5*(xGaussB[idB].+1).*(1.0./PbG[:,1]).*dnG_ub.*(0.5 .+ xGaussB[idB]*0.5) )
			#@show size(tmp2)

			KE2[e:e,:]=KE2[e:e,:]+[sum(tmp1)  sum(tmp2)]
			#readline()

		end

	end

	VI=e1; VJ=e2;
	KI2=[e1 e1]; KJ2=[edgeToNode2[:,1] edgeToNode2[:,2]]

	KI=[KI1;KI2]; KJ=[KJ1;KJ2]; KE=[KE1;KE2]
	KI=reshape(KI,size(KI,1)*size(KI,2))
	KJ=reshape(KJ,size(KJ,1)*size(KJ,2))
	KE=reshape(KE,size(KE,1)*size(KE,2))

	K=(E=KE,I=KI,J=KJ)
	V=(E=VE,I=VI,J=VJ)

	return K,V

end


text:
module Geometric

import Plots
import Images

#=
This function creates a Poincaré plot
:param n: is the array that contains the NN-inetrvals
:return: a plot object
=#
function poincare(n)
    x=[]
    y=[]
    for i in 1:length(n)-1
        push!(x,n[i])
        push!(y,n[i+1])
    end
    p=Plots.scatter(x,y,xlabel="RRn",ylabel="RRn+1",legend=false);
    return p;
end # poincare

#=
This function creates a recurrence plot
:param n: is the array that contains the NN-inetrvals
:param e: the maximum distance between two intervals, default="mean" == the mean value of succsessive differences
:return: a plot object
=#
function recurrence(n,e)
    if e=="mean"
        diff=[]
        for i in 1:length(n)-1
            push!(diff,abs(n[i+1]-n[i]))
        end
        e=sum(diff)/length(diff)
    end
    x=zeros(length(n),length(n))
    for i in 1:length(n)
        for j in i:length(n)
            if sqrt((n[i]-n[j])^2)<=e
                x[i,j]=1
                x[j,i]=1
            end
        end
    end
    img=Images.Gray.(x)
    r=Plots.plot(img);
    return r;
end # recurrence

end # module


text:
"""
    mean_square_disaplcement(df; method=[:ensemble_average, time_average], non_averaged::Bool) -> DataFrame, [DataFrame]

Compute mean_square_disaplcement over files.

**Example**
```julia
using SPT

```
"""
function _ensemble_time_average_msd(
    df::DataFrame,
    id::Symbol,
    x::Symbol,
    y::Symbol;
    return_tamsd::Bool = false,
)
    ta_msd = time_average_msd(df, id, x, y)
    return_tamsd ? (ensemble_tamsd(ta_msd), ta_msd) : ensemble_tamsd(ta_msd)
end

function _time_average_msd(
    df::DataFrame,
    id::Symbol = :TrackID,
    x::Symbol = :POSITION_X,
    y::Symbol = :POSITION_Y;
)
    tamsd = DataFrame(TrackID = Int64[], cell_name = String[], msd = Float64[], delta_t = Int64[], n = Int64[])
    @inbounds Threads.@threads for cell = sort(unique(df.cell_name))
        _data = df[df.cell_name .== cell, :]
        for n = sort(unique(_data[!, id]))
            data = extract(_data, Int(n), id, [x, y])
            T = size(data, 1)
            for δ = 1:T-1
                r² = 0.0
                @simd for t = 1:(T-δ)
                    r² += squared_displacement(data, t, δ)
                end
                r² /= (T - δ)
                push!(tamsd, [n, cell, r², δ, T - δ])
            end
        end
    end
    tamsd
end

function _ensemble_msd(df::DataFrame, id::Symbol, x::Symbol, y::Symbol)
    eamsd = DataFrame(TrackID = Int64[], cell_name = String[], msd = Float64[], n = Int[], delta_t = Int64[])
    @inbounds Threads.@threads for cell = sort(unique(df.cell_name))
        _data = df[df.cell_name .== cell, :]
        for n = sort(unique(_data[!, id]))
            data = extract(_data, Int(n), id, [x, y])
            T = size(data, 1)
            @simd for δ = 1:T-1
                push!(eamsd, [n, cell, squared_displacement(data, 1, δ), 1, δ])
            end
        end
    end
    eamsd
end

function _ensemble_tamsd(df::DataFrame)
    eatamsd = DataFrame(
        cell_name = String[],
        delta_t = Float64[],
        msd = Float64[],
        n = Int64[],
        std = Float64[],
        sem = Float64[],
        ci = Float64[],
    )
    @inbounds Threads.@threads for cell = sort(unique(df.cell_name))
        _data = df[df.cell_name .== cell, :]
        for i = 1:maximum(_data.delta_t)
            data = _data[_data.delta_t .== i, :]
            push!(
                eatamsd,
                [
                    cell,
                    i,
                    mean(data.msd),
                    sum(data.n),
                    std(data.msd),
                    sem(data.msd),
                    1.96 * sem(data.msd),
                ],
            )
        end
    end
    eatamsd
end

max_stars_repo_name: SosUts/DiffusionMoments

max_issues_repo_name: SosUts/DiffusionMoments

max_forks_repo_name: SosUts/DiffusionMoments


text:
"""
    GroupedDataFrame

The result of a [`groupby`](@ref) operation on an `AbstractDataFrame`; a
view into the `AbstractDataFrame` grouped by rows.

Not meant to be constructed directly, see `groupby`.
"""
struct GroupedDataFrame{T<:AbstractDataFrame}
    parent::T
    cols::Vector{Int}    # columns used for grouping
    groups::Vector{Int}  # group indices for each row
    idx::Vector{Int}     # indexing vector when grouped by the given columns
    starts::Vector{Int}  # starts of groups
    ends::Vector{Int}    # ends of groups
end

Base.broadcastable(::GroupedDataFrame) =
    throw(ArgumentError("broadcasting over `GroupedDataFrame`s is reserved"))

"""
    parent(gd::GroupedDataFrame)

Return the parent data frame of `gd`.
"""
Base.parent(gd::GroupedDataFrame) = getfield(gd, :parent)

#
# Split
#

"""
A view of an `AbstractDataFrame` split into row groups

```julia
groupby(d::AbstractDataFrame, cols; sort=false, skipmissing=false)
```

### Arguments

* `df` : an `AbstractDataFrame` to split
* `cols` : data frame columns to group by
* `sort` : whether to sort rows according to the values of the grouping columns `cols`
* `skipmissing` : whether to skip rows with `missing` values in one of the grouping columns `cols`

### Returns

A `GroupedDataFrame` : a grouped view into `df`

### Details

An iterator over a `GroupedDataFrame` returns a `SubDataFrame` view
for each grouping into `df`.
Within each group, the order of rows in `df` is preserved.

`cols` can be any valid data frame indexing expression.
In particular if it is an empty vector then a single-group `GroupedDataFrame`
is created.

A `GroupedDataFrame` also supports
indexing by groups, `map` (which applies a function to each group)
and `combine` (which applies a function to each group
and combines the result into a data frame).

See the following for additional split-apply-combine operations:

* [`by`](@ref) : split-apply-combine using functions
* [`aggregate`](@ref) : split-apply-combine; applies functions in the form of a cross product
* [`map`](@ref) : apply a function to each group of a `GroupedDataFrame` (without combining)
* [`combine`](@ref) : combine a `GroupedDataFrame`, optionally applying a function to each group

### Examples

```julia
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> gd = groupby(df, :a)
GroupedDataFrame with 4 groups based on key: a
First Group (2 rows): a = 1
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 1     │
│ 2   │ 1     │ 2     │ 5     │
⋮
Last Group (2 rows): a = 4
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 4     │ 1     │ 4     │
│ 2   │ 4     │ 1     │ 8     │

julia> gd[1]
2×3 SubDataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 1     │
│ 2   │ 1     │ 2     │ 5     │

julia> last(gd)
2×3 SubDataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 4     │ 1     │ 4     │
│ 2   │ 4     │ 1     │ 8     │

julia> for g in gd
           println(g)
       end
2×3 SubDataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 1     │
│ 2   │ 1     │ 2     │ 5     │
2×3 SubDataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 2     │ 1     │ 2     │
│ 2   │ 2     │ 1     │ 6     │
2×3 SubDataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 3     │ 2     │ 3     │
│ 2   │ 3     │ 2     │ 7     │
2×3 SubDataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 4     │ 1     │ 4     │
│ 2   │ 4     │ 1     │ 8     │
```

"""
function groupby(df::AbstractDataFrame, cols;
                 sort::Bool=false, skipmissing::Bool=false)
    _check_consistency(df)
    idxcols = index(df)[cols]
    intcols = idxcols isa Int ? [idxcols] : convert(Vector{Int}, idxcols)
    if isempty(intcols)
        return GroupedDataFrame(df, intcols, ones(Int, nrow(df)),
                                collect(axes(df, 1)), [1], [nrow(df)])
    end
    sdf = df[!, intcols]
    df_groups = group_rows(sdf, false, sort, skipmissing)
    GroupedDataFrame(df, intcols, df_groups.groups, df_groups.rperm,
                     df_groups.starts, df_groups.stops)
end

function Base.iterate(gd::GroupedDataFrame, i=1)
    if i > length(gd.starts)
        nothing
    else
        (view(gd.parent, gd.idx[gd.starts[i]:gd.ends[i]], :), i+1)
    end
end

Base.length(gd::GroupedDataFrame) = length(gd.starts)
Compat.lastindex(gd::GroupedDataFrame) = length(gd.starts)
Base.first(gd::GroupedDataFrame) = gd[1]
Base.last(gd::GroupedDataFrame) = gd[end]

Base.getindex(gd::GroupedDataFrame, idx::Integer) =
    view(gd.parent, gd.idx[gd.starts[idx]:gd.ends[idx]], :)

function Base.getindex(gd::GroupedDataFrame, idxs::AbstractArray)
    new_starts = gd.starts[idxs]
    new_ends = gd.ends[idxs]
    if !allunique(new_starts)
        throw(ArgumentError("duplicates in idxs argument are not allowed"))
    end
    new_groups = zeros(Int, length(gd.groups))
    for idx in eachindex(new_starts)
        @inbounds for j in new_starts[idx]:new_ends[idx]
            new_groups[gd.idx[j]] = idx
        end
    end
    GroupedDataFrame(gd.parent, gd.cols, new_groups, gd.idx, new_starts, new_ends)
end

Base.getindex(gd::GroupedDataFrame, idxs::Colon) =
    GroupedDataFrame(gd.parent, gd.cols, gd.groups, gd.idx, gd.starts, gd.ends)

function Base.:(==)(gd1::GroupedDataFrame, gd2::GroupedDataFrame)
    gd1.cols == gd2.cols &&
        length(gd1) == length(gd2) &&
        all(x -> ==(x...), zip(gd1, gd2))
end

function Base.isequal(gd1::GroupedDataFrame, gd2::GroupedDataFrame)
    isequal(gd1.cols, gd2.cols) &&
        isequal(length(gd1), length(gd2)) &&
        all(x -> isequal(x...), zip(gd1, gd2))
end

Base.names(gd::GroupedDataFrame) = names(gd.parent)
_names(gd::GroupedDataFrame) = _names(gd.parent)

"""
    map(cols => f, gd::GroupedDataFrame)
    map(f, gd::GroupedDataFrame)

Apply a function to each group of rows and return a [`GroupedDataFrame`](@ref).

If the first argument is a `cols => f` pair, `cols` must be a column name or index, or
a vector or tuple thereof, and `f` must be a callable. If `cols` is a single column index,
`f` is called with a `SubArray` view into that column for each group; else, `f` is called
with a named tuple holding `SubArray` views into these columns.

If the first argument is a vector, tuple or named tuple of such pairs, each pair is
handled as described above. If a named tuple, field names are used to name
each generated column.

If the first argument is a callable `f`, it is passed a [`SubDataFrame`](@ref) view for each group,
and the returned `DataFrame` then consists of the returned rows plus the grouping columns.
If the returned data frame contains columns with the same names as the grouping columns,
they are required to be equal.
Note that this second form is much slower than the first one due to type instability.

`f` can return a single value, a row or multiple rows. The type of the returned value
determines the shape of the resulting data frame:
- A single value gives a data frame with a single column and one row per group.
- A named tuple of single values or a [`DataFrameRow`](@ref) gives a data frame with one column
  for each field and one row per group.
- A vector gives a data frame with a single column and as many rows
  for each group as the length of the returned vector for that group.
- A data frame, a named tuple of vectors or a matrix gives a data frame
  with the same columns and as many rows for each group as the rows returned for that group.

`f` must always return the same kind of object (as defined in the above list) for
all groups, and if a named tuple or data frame, with the same fields or columns.
Named tuples cannot mix single values and vectors.
Due to type instability, returning a single value or a named tuple is dramatically
faster than returning a data frame.

As a special case, if a tuple or vector of pairs is passed as the first argument, each function
is required to return a single value or vector, which will produce each a separate column.

In all cases, the resulting `GroupedDataFrame` contains all the grouping columns in addition
to those generated by the application of `f`.
Column names are automatically generated when necessary: for functions
operating on a single column and returning a single value or vector, the function name is
appended to the input column name; for other functions, columns are called `x1`, `x2`
and so on.

Optimized methods are used when standard summary functions (`sum`, `prod`,
`minimum`, `maximum`, `mean`, `var`, `std`, `first`, `last` and `length`)
are specified using the pair syntax (e.g. `col => sum`).
When computing the `sum` or `mean` over floating point columns, results will be less
accurate than the standard [`sum`](@ref) function (which uses pairwise summation). Use
`col => x -> sum(x)` to avoid the optimized method and use the slower, more accurate one.

### Examples

```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> gd = groupby(df, :a);

julia> map(:c => sum, gd)
GroupedDataFrame{DataFrame} with 4 groups based on key: :a
First Group: 1 row
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
⋮
Last Group: 1 row
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 4     │ 12    │

julia> map(df -> sum(df.c), gd) # Slower variant
GroupedDataFrame{DataFrame} with 4 groups based on key: :a
First Group: 1 row
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
⋮
Last Group: 1 row
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 4     │ 12    │
```

See [`by`](@ref) for more examples.

### See also

`combine(f, gd)` returns a `DataFrame` rather than a `GroupedDataFrame`

"""
function Base.map(f::Any, gd::GroupedDataFrame)
    if length(gd) > 0
        idx, valscat = _combine(f, gd)
        keys = _names(gd.parent)[gd.cols]
        for key in keys
            if hasproperty(valscat, key) &&
               !isequal(valscat[!, key], view(gd.parent[!, key], idx))
               throw(ArgumentError("column :$key in returned data frame " *
                                   "is not equal to grouping key :$key"))
            end
        end
        parent = hcat!(gd.parent[idx, gd.cols],
                       without(valscat, intersect(keys, _names(valscat))))
        if length(idx) == 0
            return GroupedDataFrame(parent, collect(1:length(gd.cols)), idx,
                                    Int[], Int[], Int[])
        end
        starts = Vector{Int}(undef, length(gd))
        ends = Vector{Int}(undef, length(gd))
        starts[1] = 1
        j = 2
        @inbounds for i in 2:length(idx)
            if idx[i] != idx[i-1]
                starts[j] = i
                ends[j-1] = i - 1
                j += 1
            end
        end
        # In case some groups have to be dropped
        resize!(starts, j-1)
        resize!(ends, j-1)
        ends[end] = length(idx)
        return GroupedDataFrame(parent, collect(1:length(gd.cols)), idx,
                                collect(1:length(idx)), starts, ends)
    else
        return GroupedDataFrame(gd.parent[1:0, gd.cols], collect(1:length(gd.cols)),
                                Int[], Int[], Int[], Int[])
    end
end

"""
    combine(gd::GroupedDataFrame, cols => f...)
    combine(gd::GroupedDataFrame; (colname = cols => f)...)
    combine(gd::GroupedDataFrame, f)
    combine(f, gd::GroupedDataFrame)

Transform a [`GroupedDataFrame`](@ref) into a `DataFrame`.

If the last argument(s) consist(s) in one or more `cols => f` pair(s), or if
`colname = cols => f` keyword arguments are provided, `cols` must be
a column name or index, or a vector or tuple thereof, and `f` must be a callable.
A pair or a (named) tuple of pairs can also be provided as the first or last argument.
If `cols` is a single column index, `f` is called with a `SubArray` view into that
column for each group; else, `f` is called with a named tuple holding `SubArray`
views into these columns.

If the last argument is a callable `f`, it is passed a [`SubDataFrame`](@ref) view for each group,
and the returned `DataFrame` then consists of the returned rows plus the grouping columns.
If the returned data frame contains columns with the same names as the grouping columns,
they are required to be equal.
Note that this second form is much slower than the first one due to type instability.
A method is defined with `f` as the first argument, so do-block
notation can be used.

`f` can return a single value, a row or multiple rows. The type of the returned value
determines the shape of the resulting data frame:
- A single value gives a data frame with a single column and one row per group.
- A named tuple of single values or a [`DataFrameRow`](@ref) gives a data frame with one column
  for each field and one row per group.
- A vector gives a data frame with a single column and as many rows
  for each group as the length of the returned vector for that group.
- A data frame, a named tuple of vectors or a matrix gives a data frame
  with the same columns and as many rows for each group as the rows returned for that group.

`f` must always return the same kind of object (as defined in the above list) for
all groups, and if a named tuple or data frame, with the same fields or columns.
Named tuples cannot mix single values and vectors.
Due to type instability, returning a single value or a named tuple is dramatically
faster than returning a data frame.

As a special case, if a tuple or vector of pairs is passed as the first argument, each function
is required to return a single value or vector, which will produce each a separate column.

In all cases, the resulting data frame contains all the grouping columns in addition
to those generated by the application of `f`.
Column names are automatically generated when necessary: for functions
operating on a single column and returning a single value or vector, the function name is
appended to the input column name; for other functions, columns are called `x1`, `x2`
and so on. The resulting data frame will be sorted if `sort=true` was passed to the
[`groupby`](@ref) call from which `gd` was constructed. Otherwise, ordering of rows
is undefined.

Optimized methods are used when standard summary functions (`sum`, `prod`,
`minimum`, `maximum`, `mean`, `var`, `std`, `first`, `last` and `length`)
are specified using the pair syntax (e.g. `col => sum`).
When computing the `sum` or `mean` over floating point columns, results will be less
accurate than the standard [`sum`](@ref) function (which uses pairwise summation). Use
`col => x -> sum(x)` to avoid the optimized method and use the slower, more accurate one.

### Examples

```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> gd = groupby(df, :a);

julia> combine(gd, :c => sum)
4×2 DataFrame
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> combine(:c => sum, gd)
4×2 DataFrame
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> combine(df -> sum(df.c), gd) # Slower variant
4×2 DataFrame
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │
```

See [`by`](@ref) for more examples.

### See also

[`by(f, df, cols)`](@ref) is a shorthand for `combine(f, groupby(df, cols))`.

[`map`](@ref): `combine(f, groupby(df, cols))` is a more efficient equivalent
of `combine(map(f, groupby(df, cols)))`.

"""
function combine(f::Any, gd::GroupedDataFrame)
    if length(gd) > 0
        idx, valscat = _combine(f, gd)
        keys = _names(gd.parent)[gd.cols]
        for key in keys
            if hasproperty(valscat, key) &&
               !isequal(valscat[!, key], view(gd.parent[!, key], idx))
               throw(ArgumentError("column :$key in returned data frame " *
                                   "is not equal to grouping key :$key"))
            end
        end
        return hcat!(gd.parent[idx, gd.cols],
                     without(valscat, intersect(keys, _names(valscat))))
    else
        return gd.parent[1:0, gd.cols]
    end
end
combine(gd::GroupedDataFrame, f::Any) = combine(f, gd)
combine(gd::GroupedDataFrame, f::Pair...) = combine(f, gd)
combine(gd::GroupedDataFrame, f::Pair) = combine(f, gd)

function combine(gd::GroupedDataFrame; f...)
    if length(f) == 0
        Base.depwarn("combine(gd) is deprecated, use DataFrame(gd) instead", :combine)
        combine(identity, gd)
    else
        combine(values(f), gd)
    end
end

# Wrapping automatically adds column names when the value returned
# by the user-provided function lacks them
wrap(x::Union{AbstractDataFrame, NamedTuple, DataFrameRow}) = x
wrap(x::AbstractMatrix) =
    NamedTuple{Tuple(gennames(size(x, 2)))}(Tuple(view(x, :, i) for i in 1:size(x, 2)))
wrap(x::Any) = (x1=x,)

function do_call(f::Any, gd::GroupedDataFrame, incols::AbstractVector, i::Integer)
    idx = gd.idx[gd.starts[i]:gd.ends[i]]
    f(view(incols, idx))
end
function do_call(f::Any, gd::GroupedDataFrame, incols::NamedTuple, i::Integer)
    idx = gd.idx[gd.starts[i]:gd.ends[i]]
    f(map(c -> view(c, idx), incols))
end
do_call(f::Any, gd::GroupedDataFrame, incols::Nothing, i::Integer) =
    f(gd[i])

_nrow(df::AbstractDataFrame) = nrow(df)
_nrow(x::NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}) =
    isempty(x) ? 0 : length(x[1])
_ncol(df::AbstractDataFrame) = ncol(df)
_ncol(x::Union{NamedTuple, DataFrameRow}) = length(x)

abstract type AbstractAggregate end

struct Reduce{O, C, A} <: AbstractAggregate
    op::O
    condf::C
    adjust::A
end
Reduce(f, condf=nothing) = Reduce(f, condf, nothing)

check_aggregate(f::Any) = f
check_aggregate(::typeof(sum)) = Reduce(Base.add_sum)
check_aggregate(::typeof(prod)) = Reduce(Base.mul_prod)
check_aggregate(::typeof(maximum)) = Reduce(max)
check_aggregate(::typeof(minimum)) = Reduce(min)
check_aggregate(::typeof(mean)) = Reduce(Base.add_sum, nothing, /)
check_aggregate(::typeof(sum∘skipmissing)) = Reduce(Base.add_sum, !ismissing)
check_aggregate(::typeof(prod∘skipmissing)) = Reduce(Base.mul_prod, !ismissing)
check_aggregate(::typeof(maximum∘skipmissing)) = Reduce(max, !ismissing)
check_aggregate(::typeof(minimum∘skipmissing)) = Reduce(min, !ismissing)
check_aggregate(::typeof(mean∘skipmissing)) = Reduce(Base.add_sum, !ismissing, /)

# Other aggregate functions which are not strictly reductions
struct Aggregate{F, C} <: AbstractAggregate
    f::F
    condf::C
end
Aggregate(f) = Aggregate(f, nothing)

check_aggregate(::typeof(var)) = Aggregate(var)
check_aggregate(::typeof(var∘skipmissing)) = Aggregate(var, !ismissing)
check_aggregate(::typeof(std)) = Aggregate(std)
check_aggregate(::typeof(std∘skipmissing)) = Aggregate(std, !ismissing)
check_aggregate(::typeof(first)) = Aggregate(first)
check_aggregate(::typeof(first∘skipmissing)) = Aggregate(first, !ismissing)
check_aggregate(::typeof(last)) = Aggregate(last)
check_aggregate(::typeof(last∘skipmissing)) = Aggregate(last, !ismissing)
check_aggregate(::typeof(length)) = Aggregate(length)
# SkipMissing does not support length

# Find first value matching condition for each group
# Optimized for situations where a matching value is typically encountered
# among the first rows for each group
function fillfirst!(condf, outcol::AbstractVector, incol::AbstractVector,
                    gd::GroupedDataFrame; rev::Bool=false)
    nfilled = 0
    @inbounds for i in eachindex(outcol)
        s = gd.starts[i]
        offsets = rev ? (nrow(gd[i])-1:-1:0) : (0:nrow(gd[i])-1)
        for j in offsets
            x = incol[gd.idx[s+j]]
            if !condf === nothing || condf(x)
                outcol[i] = x
                nfilled += 1
                break
            end
        end
    end
    if nfilled < length(outcol)
        throw(ArgumentError("some groups contain only missing values"))
    end
    outcol
end

# Use reducedim_init to get a vector of the right type,
# but trick it into using the expected length
groupreduce_init(op, condf, incol, gd) =
    Base.reducedim_init(identity, op, view(incol, 1:length(gd)), 2)
for (op, initf) in ((:max, :typemin), (:min, :typemax))
    @eval begin
        function groupreduce_init(::typeof($op), condf, incol::AbstractVector{T}, gd) where T
            # !ismissing check is purely an optimization to avoid a copy later
            outcol = similar(incol, condf === !ismissing ? nonmissingtype(T) : T, length(gd))
            # Comparison is possible only between CatValues from the same pool
            if incol isa CategoricalVector
                U = Union{CategoricalArrays.leveltype(outcol),
                          eltype(outcol) >: Missing ? Missing : Union{}}
                outcol = CategoricalArray{U, 1}(outcol.refs, incol.pool)
            end
            # It is safe to use a non-missing init value
            # since missing will poison the result if present
            S = nonmissingtype(T)
            if isconcretetype(S) && hasmethod($initf, Tuple{S})
                fill!(outcol, $initf(S))
            elseif condf !== nothing
                fillfirst!(condf, outcol, incol, gd)
            else
                @inbounds for i in eachindex(outcol)
                    outcol[i] = incol[gd.idx[gd.starts[i]]]
                end
            end
            return outcol
        end
    end
end

function copyto_widen!(res::AbstractVector{T}, x::AbstractVector) where T
    @inbounds for i in eachindex(res, x)
        val = x[i]
        S = typeof(val)
        if S <: T || promote_type(S, T) <: T
            res[i] = val
        else
            newres = Tables.allocatecolumn(promote_type(S, T), length(x))
            return copyto_widen!(newres, x)
        end
    end
    return res
end

function groupreduce!(res, f, op, condf, adjust,
                      incol::AbstractVector{T}, gd::GroupedDataFrame) where T
    n = length(gd)
    if adjust !== nothing
        counts = zeros(Int, n)
    end
    @inbounds for i in eachindex(incol, gd.groups)
        gix = gd.groups[i]
        x = incol[i]
        if gix > 0 && (condf === nothing || condf(x))
            res[gix] = op(res[gix], f(x, gix))
            adjust !== nothing && (counts[gix] += 1)
        end
    end
    outcol = adjust === nothing ? res : map(adjust, res, counts)
    # Undo pool sharing done by groupreduce_init
    if outcol isa CategoricalVector
        U = Union{CategoricalArrays.leveltype(outcol),
                  eltype(outcol) >: Missing ? Missing : Union{}}
        outcol = CategoricalArray{U, 1}(outcol.refs, incol.pool)
    end
    if isconcretetype(eltype(outcol))
        return outcol
    else
        copyto_widen!(Tables.allocatecolumn(typeof(first(outcol)), n), outcol)
    end
end

# Function barrier works around type instability of _groupreduce_init due to applicable
groupreduce(f, op, condf, adjust, incol::AbstractVector, gd::GroupedDataFrame) =
    groupreduce!(groupreduce_init(op, condf, incol, gd),
                 f, op, condf, adjust, incol, gd)
# Avoids the overhead due to Missing when computing reduction
groupreduce(f, op, condf::typeof(!ismissing), adjust,
            incol::AbstractVector, gd::GroupedDataFrame) =
    groupreduce!(disallowmissing(groupreduce_init(op, condf, incol, gd)),
                 f, op, condf, adjust, incol, gd)

(r::Reduce)(incol::AbstractVector, gd::GroupedDataFrame) =
    groupreduce((x, i) -> x, r.op, r.condf, r.adjust, incol, gd)

function (agg::Aggregate{typeof(var)})(incol::AbstractVector, gd::GroupedDataFrame)
    means = groupreduce((x, i) -> x, Base.add_sum, agg.condf, /, incol, gd)
    # !ismissing check is purely an optimization to avoid a copy later
    if eltype(means) >: Missing && agg.condf !== !ismissing
        T = Union{Missing, real(eltype(means))}
    else
        T = real(eltype(means))
    end
    res = zeros(T, length(gd))
    groupreduce!(res, (x, i) -> @inbounds(abs2(x - means[i])), +,
                 agg.condf, (x, l) -> x / (l-1), incol, gd)
end

function (agg::Aggregate{typeof(std)})(incol::AbstractVector, gd::GroupedDataFrame)
    outcol = Aggregate(var, agg.condf)(incol, gd)
    map!(sqrt, outcol, outcol)
end

for f in (first, last)
    function (agg::Aggregate{typeof(f)})(incol::AbstractVector, gd::GroupedDataFrame)
        n = length(gd)
        outcol = similar(incol, n)
        if agg.condf === !ismissing
            fillfirst!(agg.condf, outcol, incol, gd, rev=agg.f === last)
        else
            v = agg.f === first ? gd.starts : gd.ends
            map!(i -> incol[gd.idx[v[i]]], outcol, 1:n)
        end
        if isconcretetype(eltype(outcol))
            return outcol
        else
            return copyto_widen!(Tables.allocatecolumn(typeof(first(outcol)), n), outcol)
        end
    end
end

(agg::Aggregate{typeof(length)})(incol::AbstractVector, gd::GroupedDataFrame) =
    gd.ends .- gd.starts .+ 1

function do_f(f, x...)
    @inline function fun(x...)
        res = f(x...)
        if res isa Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix}
            throw(ArgumentError("a single value or vector result is required when passing " *
                                "a vector or tuple of functions (got $(typeof(res)))"))
        end
        res
    end
end

# Avoid recompilation of larger function for each set of names
function _combine(@nospecialize(f::NamedTuple{<:Any, <:Tuple{Vararg{Pair}}}),
                  gd::GroupedDataFrame)
    idx, valscat = _combine(collect(f), gd)
    rename!(valscat, collect(Symbol, propertynames(f)))
    idx, valscat
end

# Avoid recompilation of larger function for each length
_combine(@nospecialize(f::Tuple{Vararg{Pair}}), gd::GroupedDataFrame) =
    _combine(collect(f), gd)

function _combine(f::AbstractVector{<:Pair}, gd::GroupedDataFrame)
    res = map(f) do p
        agg = check_aggregate(last(p))
        if agg isa AbstractAggregate && p isa Pair{<:ColumnIndex}
            incol = gd.parent[!, first(p)]
            idx = gd.idx[gd.starts]
            outcol = agg(incol, gd)
            return idx, outcol
        else
            fun = do_f(last(p))
            if p isa Pair{<:ColumnIndex}
                incols = gd.parent[!, first(p)]
            else
                df = gd.parent[!, collect(first(p))]
                incols = NamedTuple{Tuple(names(df))}(eachcol(df))
            end
            firstres = do_call(fun, gd, incols, 1)
            idx, outcols, _ = _combine_with_first(wrap(firstres), fun, gd, incols)
            return idx, outcols[1]
        end
    end
    # TODO: avoid recomputing idx for each pair
    idx = res[1][1]
    outcols = map(x -> x[2], res)
    if !all(x -> length(x) == length(outcols[1]), outcols)
        throw(ArgumentError("all functions must return values of the same length"))
    end
    nams = [f[i] isa Pair{<:ColumnIndex} ?
                Symbol(names(gd.parent)[index(gd.parent)[first(f[i])]],
                       '_', funname(last(f[i]))) :
                Symbol('x', i)
            for i in 1:length(f)]
    valscat = DataFrame(collect(AbstractVector, outcols), nams, makeunique=true)
    return idx, valscat
end

function _combine(f::Any, gd::GroupedDataFrame)
    if f isa Pair{<:ColumnIndex}
        incols = gd.parent[!, first(f)]
        fun = last(f)
    elseif f isa Pair
        df = gd.parent[! , collect(first(f))]
        incols = NamedTuple{Tuple(names(df))}(eachcol(df))
        fun = last(f)
    else
        incols = nothing
        fun = f
    end
    agg = check_aggregate(fun)
    if agg isa AbstractAggregate && f isa Pair{<:ColumnIndex}
        idx = gd.idx[gd.starts]
        outcols = (agg(incols, gd),)
        # nms is set below
    else
        firstres = do_call(fun, gd, incols, 1)
        idx, outcols, nms = _combine_with_first(wrap(firstres), fun, gd, incols)
    end
    if f isa Pair{<:ColumnIndex} &&
        (agg isa AbstractAggregate ||
         !isa(firstres, Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix}))
         nms = [Symbol(names(gd.parent)[index(gd.parent)[first(f)]], '_', funname(fun))]
    end
    valscat = DataFrame(collect(AbstractVector, outcols), collect(Symbol, nms))
    return idx, valscat
end

function _combine_with_first(first::Union{NamedTuple, DataFrameRow, AbstractDataFrame},
                             f::Any, gd::GroupedDataFrame,
                             incols::Union{Nothing, AbstractVector, NamedTuple})
    if first isa AbstractDataFrame
        n = 0
        eltys = eltype.(eachcol(first))
    elseif first isa NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}
        n = 0
        eltys = map(eltype, first)
    elseif first isa DataFrameRow
        n = length(gd)
        eltys = eltype.(eachcol(parent(first)))
    else # NamedTuple giving a single row
        n = length(gd)
        eltys = map(typeof, first)
        if any(x -> x <: AbstractVector, eltys)
            throw(ArgumentError("mixing single values and vectors in a named tuple is not allowed"))
        end
    end
    idx = Vector{Int}(undef, n)
    local initialcols
    let eltys=eltys, n=n # Workaround for julia#15276
        initialcols = ntuple(i -> Tables.allocatecolumn(eltys[i], n), _ncol(first))
    end
    outcols = _combine_with_first!(first, initialcols, idx, 1, 1, f, gd, incols,
                                   tuple(propertynames(first)...))
    idx, outcols, propertynames(first)
end

function fill_row!(row, outcols::NTuple{N, AbstractVector},
                   i::Integer, colstart::Integer,
                   colnames::NTuple{N, Symbol}) where N
    if !isa(row, Union{NamedTuple, DataFrameRow})
        throw(ArgumentError("return value must not change its kind " *
                            "(single row or variable number of rows) across groups"))
    elseif row isa NamedTuple && any(x -> x isa AbstractVector, row)
        throw(ArgumentError("mixing single values and vectors in a named tuple is not allowed"))
    elseif _ncol(row) != N
        throw(ArgumentError("return value must have the same number of columns " *
                            "for all groups (got $N and $(length(row)))"))
    end
    @inbounds for j in colstart:length(outcols)
        col = outcols[j]
        cn = colnames[j]
        local val
        try
            val = row[cn]
        catch
            throw(ArgumentError("return value must have the same column names " *
                                "for all groups (got $colnames and $(propertynames(row)))"))
        end
        S = typeof(val)
        T = eltype(col)
        if S <: T || promote_type(S, T) <: T
            col[i] = val
        else
            return j
        end
    end
    return nothing
end

function _combine_with_first!(first::Union{NamedTuple, DataFrameRow},
                              outcols::NTuple{N, AbstractVector},
                              idx::Vector{Int}, rowstart::Integer, colstart::Integer,
                              f::Any, gd::GroupedDataFrame,
                              incols::Union{Nothing, AbstractVector, NamedTuple},
                              colnames::NTuple{N, Symbol}) where N
    len = length(gd)
    # Handle first group
    j = fill_row!(first, outcols, rowstart, colstart, colnames)
    @assert j === nothing # eltype is guaranteed to match
    idx[rowstart] = gd.idx[gd.starts[rowstart]]
    # Handle remaining groups
    @inbounds for i in rowstart+1:len
        row = wrap(do_call(f, gd, incols, i))
        j = fill_row!(row, outcols, i, 1, colnames)
        if j !== nothing # Need to widen column type
            local newcols
            let i = i, j = j, outcols=outcols, row=row # Workaround for julia#15276
                newcols = ntuple(length(outcols)) do k
                    S = typeof(row[k])
                    T = eltype(outcols[k])
                    U = promote_type(S, T)
                    if S <: T || U <: T
                        outcols[k]
                    else
                        copyto!(Tables.allocatecolumn(U, length(outcols[k])),
                                1, outcols[k], 1, k >= j ? i-1 : i)
                    end
                end
            end
            return _combine_with_first!(row, newcols, idx, i, j, f, gd, incols, colnames)
        end
        idx[i] = gd.idx[gd.starts[i]]
    end
    outcols
end

# This needs to be in a separate function
# to work around a crash due to JuliaLang/julia#29430
if VERSION >= v"1.1.0-DEV.723"
    @inline function do_append!(do_it, col, vals)
        do_it && append!(col, vals)
        return do_it
    end
else
    @noinline function do_append!(do_it, col, vals)
        do_it && append!(col, vals)
        return do_it
    end
end

function append_rows!(rows, outcols::NTuple{N, AbstractVector},
                      colstart::Integer, colnames::NTuple{N, Symbol}) where N
    if !isa(rows, Union{AbstractDataFrame, NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}})
        throw(ArgumentError("return value must not change its kind " *
                            "(single row or variable number of rows) across groups"))
    elseif _ncol(rows) != N
        throw(ArgumentError("return value must have the same number of columns " *
                            "for all groups (got $N and $(_ncol(rows)))"))
    end
    @inbounds for j in colstart:length(outcols)
        col = outcols[j]
        cn = colnames[j]
        local vals
        try
            vals = getproperty(rows, cn)
        catch
            throw(ArgumentError("return value must have the same column names " *
                                "for all groups (got $(Tuple(colnames)) and $(Tuple(names(rows))))"))
        end
        S = eltype(vals)
        T = eltype(col)
        if !do_append!(S <: T || promote_type(S, T) <: T, col, vals)
            return j
        end
    end
    return nothing
end

function _combine_with_first!(first::Union{AbstractDataFrame,
                                           NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}},
                              outcols::NTuple{N, AbstractVector},
                              idx::Vector{Int}, rowstart::Integer, colstart::Integer,
                              f::Any, gd::GroupedDataFrame,
                              incols::Union{Nothing, AbstractVector, NamedTuple},
                              colnames::NTuple{N, Symbol}) where N
    len = length(gd)
    # Handle first group
    j = append_rows!(first, outcols, colstart, colnames)
    @assert j === nothing # eltype is guaranteed to match
    append!(idx, Iterators.repeated(gd.idx[gd.starts[rowstart]], _nrow(first)))
    # Handle remaining groups
    @inbounds for i in rowstart+1:len
        rows = wrap(do_call(f, gd, incols, i))
        j = append_rows!(rows, outcols, 1, colnames)
        if j !== nothing # Need to widen column type
            local newcols
            let i = i, j = j, outcols=outcols, rows=rows # Workaround for julia#15276
                newcols = ntuple(length(outcols)) do k
                    S = eltype(rows isa AbstractDataFrame ? rows[!, k] : rows[k])
                    T = eltype(outcols[k])
                    U = promote_type(S, T)
                    if S <: T || U <: T
                        outcols[k]
                    else
                        copyto!(Tables.allocatecolumn(U, length(outcols[k])), outcols[k])
                    end
                end
            end
            return _combine_with_first!(rows, newcols, idx, i, j, f, gd, incols, colnames)
        end
        append!(idx, Iterators.repeated(gd.idx[gd.starts[i]], _nrow(rows)))
    end
    outcols
end

"""
    by(df::AbstractDataFrame, keys, cols=>f...;
       sort::Bool=false, skipmissing::Bool=false)
    by(df::AbstractDataFrame, keys; (colname = cols => f)...,
       sort::Bool=false, skipmissing::Bool=false)
    by(df::AbstractDataFrame, keys, f;
       sort::Bool=false, skipmissing::Bool=false)
    by(f, df::AbstractDataFrame, keys;
       sort::Bool=false, skipmissing::Bool=false)

Split-apply-combine in one step: apply `f` to each grouping in `df`
based on grouping columns `keys`, and return a `DataFrame`.

`keys` can be either a single column index, or a vector thereof.

If the last argument(s) consist(s) in one or more `cols => f` pair(s), or if
`colname = cols => f` keyword arguments are provided, `cols` must be
a column name or index, or a vector or tuple thereof, and `f` must be a callable.
A pair or a (named) tuple of pairs can also be provided as the first or last argument.
If `cols` is a single column index, `f` is called with a `SubArray` view into that
column for each group; else, `f` is called with a named tuple holding `SubArray`
views into these columns.

If the last argument is a callable `f`, it is passed a [`SubDataFrame`](@ref) view for each group,
and the returned `DataFrame` then consists of the returned rows plus the grouping columns.
If the returned data frame contains columns with the same names as the grouping columns,
they are required to be equal.
Note that this second form is much slower than the first one due to type instability.
A method is defined with `f` as the first argument, so do-block
notation can be used.

`f` can return a single value, a row or multiple rows. The type of the returned value
determines the shape of the resulting data frame:
- A single value gives a data frame with a single column and one row per group.
- A named tuple of single values or a [`DataFrameRow`](@ref) gives a data frame with one column
  for each field and one row per group.
- A vector gives a data frame with a single column and as many rows
  for each group as the length of the returned vector for that group.
- A data frame, a named tuple of vectors or a matrix gives a data frame
  with the same columns and as many rows for each group as the rows returned for that group.

`f` must always return the same kind of object (as defined in the above list) for
all groups, and if a named tuple or data frame, with the same fields or columns.
Named tuples cannot mix single values and vectors.
Due to type instability, returning a single value or a named tuple is dramatically
faster than returning a data frame.

As a special case, if multiple pairs are passed as last arguments, each function
is required to return a single value or vector, which will produce each a separate column.

In all cases, the resulting data frame contains all the grouping columns in addition
to those generated by the application of `f`.
Column names are automatically generated when necessary: for functions
operating on a single column and returning a single value or vector, the function name is
appended to the input colummn name; for other functions, columns are called `x1`, `x2`
and so on. The resulting data frame will be sorted on `keys` if `sort=true`.
Otherwise, ordering of rows is undefined.
If `skipmissing=true` then the resulting data frame will not contain groups
with `missing` values in one of the `keys` columns.

Optimized methods are used when standard summary functions (`sum`, `prod`,
`minimum`, `maximum`, `mean`, `var`, `std`, `first`, `last` and `length)
are specified using the pair syntax (e.g. `col => sum`).
When computing the `sum` or `mean` over floating point columns, results will be less
accurate than the standard [`sum`](@ref) function (which uses pairwise summation). Use
`col => x -> sum(x)` to avoid the optimized method and use the slower, more accurate one.

`by(d, cols, f)` is equivalent to `combine(f, groupby(d, cols))` and to the
less efficient `combine(map(f, groupby(d, cols)))`.

### Examples

```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> by(df, :a, :c => sum)
4×2 DataFrame
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> by(df, :a, d -> sum(d.c)) # Slower variant
4×2 DataFrame
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> by(df, :a) do d # do syntax for the slower variant
           sum(d.c)
       end
4×2 DataFrame
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> by(df, :a, :c => x -> 2 .* x)
8×2 DataFrame
│ Row │ a     │ c_function │
│     │ Int64 │ Int64      │
├─────┼───────┼────────────┤
│ 1   │ 1     │ 2          │
│ 2   │ 1     │ 10         │
│ 3   │ 2     │ 4          │
│ 4   │ 2     │ 12         │
│ 5   │ 3     │ 6          │
│ 6   │ 3     │ 14         │
│ 7   │ 4     │ 8          │
│ 8   │ 4     │ 16         │

julia> by(df, :a, c_sum = :c => sum, c_sum2 = :c => x -> sum(x.^2))
4×3 DataFrame
│ Row │ a     │ c_sum │ c_sum2 │
│     │ Int64 │ Int64 │ Int64  │
├─────┼───────┼───────┼────────┤
│ 1   │ 1     │ 6     │ 26     │
│ 2   │ 2     │ 8     │ 40     │
│ 3   │ 3     │ 10    │ 58     │
│ 4   │ 4     │ 12    │ 80     │

julia> by(df, :a, (:b, :c) => x -> (minb = minimum(x.b), sumc = sum(x.c)))
4×3 DataFrame
│ Row │ a     │ minb  │ sumc  │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 6     │
│ 2   │ 2     │ 1     │ 8     │
│ 3   │ 3     │ 2     │ 10    │
│ 4   │ 4     │ 1     │ 12    │
```

"""
by(d::AbstractDataFrame, cols::Any, f::Any;
   sort::Bool=false, skipmissing::Bool=false) =
    combine(f, groupby(d, cols, sort=sort, skipmissing=skipmissing))
by(f::Any, d::AbstractDataFrame, cols::Any;
   sort::Bool=false, skipmissing::Bool=false) =
    by(d, cols, f, sort=sort, skipmissing=skipmissing)
by(d::AbstractDataFrame, cols::Any, f::Pair;
   sort::Bool=false, skipmissing::Bool=false) =
    combine(f, groupby(d, cols, sort=sort, skipmissing=skipmissing))
by(d::AbstractDataFrame, cols::Any, f::Pair...;
   sort::Bool=false, skipmissing::Bool=false) =
    combine(f, groupby(d, cols, sort=sort, skipmissing=skipmissing))
by(d::AbstractDataFrame, cols::Any;
   sort::Bool=false, skipmissing::Bool=false, f...) =
    combine(values(f), groupby(d, cols, sort=sort, skipmissing=skipmissing))

#
# Aggregate convenience functions
#

# Applies a set of functions over a DataFrame, in the from of a cross-product
"""
Split-apply-combine that applies a set of functions over columns of an
`AbstractDataFrame` or [`GroupedDataFrame`](@ref)

```julia
aggregate(df::AbstractDataFrame, fs)
aggregate(df::AbstractDataFrame, cols, fs; sort=false, skipmissing=false)
aggregate(gd::GroupedDataFrame, fs; sort=false)
```

### Arguments

* `df` : an `AbstractDataFrame`
* `gd` : a `GroupedDataFrame`
* `cols` : a column indicator (`Symbol`, `Int`, `Vector{Symbol}`, etc.)
* `fs` : a function or vector of functions to be applied to vectors
  within groups; expects each argument to be a column vector
* `sort` : whether to sort rows according to the values of the grouping columns
* `skipmissing` : whether to skip rows with `missing` values in one of the grouping columns `cols`

Each `fs` should return a value or vector. All returns must be the
same length.

### Returns

* `::DataFrame`

### Examples

```jldoctest
julia> using Statistics

julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> aggregate(df, :a, sum)
4×3 DataFrame
│ Row │ a     │ b_sum │ c_sum │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 4     │ 6     │
│ 2   │ 2     │ 2     │ 8     │
│ 3   │ 3     │ 4     │ 10    │
│ 4   │ 4     │ 2     │ 12    │

julia> aggregate(df, :a, [sum, x->mean(skipmissing(x))])
4×5 DataFrame
│ Row │ a     │ b_sum │ c_sum │ b_function │ c_function │
│     │ Int64 │ Int64 │ Int64 │ Float64    │ Float64    │
├─────┼───────┼───────┼───────┼────────────┼────────────┤
│ 1   │ 1     │ 4     │ 6     │ 2.0        │ 3.0        │
│ 2   │ 2     │ 2     │ 8     │ 1.0        │ 4.0        │
│ 3   │ 3     │ 4     │ 10    │ 2.0        │ 5.0        │
│ 4   │ 4     │ 2     │ 12    │ 1.0        │ 6.0        │

julia> aggregate(groupby(df, :a), [sum, x->mean(skipmissing(x))])
4×5 DataFrame
│ Row │ a     │ b_sum │ c_sum │ b_function │ c_function │
│     │ Int64 │ Int64 │ Int64 │ Float64    │ Float64    │
├─────┼───────┼───────┼───────┼────────────┼────────────┤
│ 1   │ 1     │ 4     │ 6     │ 2.0        │ 3.0        │
│ 2   │ 2     │ 2     │ 8     │ 1.0        │ 4.0        │
│ 3   │ 3     │ 4     │ 10    │ 2.0        │ 5.0        │
│ 4   │ 4     │ 2     │ 12    │ 1.0        │ 6.0        │
```
"""
aggregate(d::AbstractDataFrame, fs::Any; sort::Bool=false) =
    aggregate(d, [fs], sort=sort)
function aggregate(d::AbstractDataFrame, fs::AbstractVector; sort::Bool=false)
    headers = _makeheaders(fs, _names(d))
    _aggregate(d, fs, headers, sort)
end

# Applies aggregate to non-key cols of each SubDataFrame of a GroupedDataFrame
aggregate(gd::GroupedDataFrame, f::Any; sort::Bool=false) = aggregate(gd, [f], sort=sort)
function aggregate(gd::GroupedDataFrame, fs::AbstractVector; sort::Bool=false)
    headers = _makeheaders(fs, setdiff(_names(gd), _names(gd.parent)[gd.cols]))
    res = combine(x -> _aggregate(without(x, gd.cols), fs, headers), gd)
    sort && sort!(res, headers)
    res
end

# Groups DataFrame by cols before applying aggregate
function aggregate(d::AbstractDataFrame, cols, fs::Any;
                   sort::Bool=false, skipmissing::Bool=false)
    aggregate(groupby(d, cols, sort=sort, skipmissing=skipmissing), fs)
end

function funname(f)
    n = nameof(f)
    String(n)[1] == '#' ? :function : n
end

_makeheaders(fs::AbstractVector, cn::AbstractVector{Symbol}) =
    [Symbol(colname, '_', funname(f)) for f in fs for colname in cn]

function _aggregate(d::AbstractDataFrame, fs::AbstractVector,
                    headers::AbstractVector{Symbol}, sort::Bool=false)
    res = DataFrame(AbstractVector[vcat(f(d[!, i])) for f in fs for i in 1:size(d, 2)], headers, makeunique=true)
    sort && sort!(res, headers)
    res
end

function DataFrame(gd::GroupedDataFrame; copycols::Bool=true)
    if !copycols
        throw(ArgumentError("It is not possible to construct a `DataFrame`" *
                            "from GroupedDataFrame with `copycols=false`"))
    end
    length(gd) == 0 && return similar(parent(gd), 0)
    idx = similar(gd.idx)
    doff = 1
    for (s,e) in zip(gd.starts, gd.ends)
        n = e - s + 1
        copyto!(idx, doff, gd.idx, s, n)
        doff += n
    end
    resize!(idx, doff - 1)
    parent(gd)[idx, :]
end

"""
    groupindices(gd::GroupedDataFrame)

Return a vector of group indices for each row of `parent(gd)`.

Rows appearing in group `gd[i]` are attributed index `i`. Rows not present in
any group are attributed `missing` (this can happen if `skipmissing=true` was
passed when creating `gd`, or if `gd` is a subset from a larger [`GroupedDataFrame`](@ref)).
"""
groupindices(gd::GroupedDataFrame) = replace(gd.groups, 0=>missing)

"""
    groupvars(gd::GroupedDataFrame)

Return a vector of column names in `parent(gd)` used for grouping.
"""
groupvars(gd::GroupedDataFrame) = _names(gd)[gd.cols]


text:
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux

using Plots
plotly()

include("OdeModel.jl")
## ODE MODEL
hill = OdeModel.hill
@parameters t k₁ k₂ k₃ k₄ k₅ k₆ KM₁ KM₂ KM₃ KM₄ n₁ n₂ n₃ n₄ α β S_tot

## Fokker planck
@parameters H, S, θ
@parameters D_H, D_S
@derivatives Dt'~t
@derivatives Dss''~S
@derivatives Dhh''~H
@derivatives Ds'~S
@derivatives Dh'~H
@variables p(..)

F_H = k₁ * (1-α)*hill(S, KM₁, n₁) * hill(H, KM₂, n₂) + k₂ - k₃*H
F_S = k₄ * (1-β)*hill(H, KM₃,n₃) * hill(S, KM₄, n₄)*(S_tot - S) + k₅ - k₆*S

# Replace with number
p_ =[k₁ => 40,
    k₂ => 1,
    k₃ => 0.004,
    k₄ => 0.02,
    k₅ => 0.04,
    k₆ => 0.004,
    KM₁ => 280,
    KM₂ => 4250,
    KM₃ => 2200,
    KM₄ => 90,
    n₁ => 2,
    n₂ => 4,
    n₃ => 2,
    n₄ => 4,
    α => 0.35,
    β => 0.95,
    S_tot => 225,
    D_H=>500,
    D_S=>0.45]

F_H = substitute.(F_H, (p_,))[1]
F_S = substitute.(F_S, (p_,))[1]

σ_H, σ_S = 5000.0, 150.0
H0, S0 = 3275.0, 112.0

eq = Dt(p(S,H,t,θ)) ~ - Dh(F_H*p(S,H,t,θ)) - Ds(F_S*p(S,H,t,θ)) + 500.0 * Dhh(p(S,H,t,θ)) + 0.45 *Dss(p(S,H,t,θ))

# Domains
domains = [ S ∈ IntervalDomain(0.0,112.0),
            H ∈ IntervalDomain(0.0,3275.0),
            t ∈ IntervalDomain(0.0,300.0)]
S_max = domains[1].domain.upper
H_max = domains[2].domain.upper

# boundary conditions
norm2(X, a, b) = (X-a)^2 / b^2
bcs = [
    p(S,H, 0, θ) ~  (2*pi*σ_H*σ_S)^-1 * exp( -norm2(H,H0, σ_H) - norm2(S,S0, σ_S)  ),
    p(0,H,t,θ) ~ 0.f0,
    p(S,0,t,θ) ~ 0.f0,
    p(S_max,H,t,θ) ~ 0.f0,
    p(S,H_max,t,θ) ~ 0.f0,
]

# Discretization
dS = 10.0; dH=100.0; dt = 10.0

# Nerual Network
dim = 1

chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

initθ = initial_params(chain)

discretization = NeuralPDE.PhysicsInformedNN([dS,dH,dt],
                                             chain,
                                             initθ,
                                             strategy = NeuralPDE.StochasticTraining(include_frac=0.9))

pde_system = PDESystem(eq,bcs,domains,[S,H,t],[p])

prob = NeuralPDE.discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

@time res = GalacticOptim.solve(prob, ADAM(0.1), progress = false; cb = cb, maxiters=3000)

phi = discretization.phi

## Create Terminal PDE Problem

TerminalPDEProblem

max_stars_repo_name: stevengogogo/ProgrammableAging

max_issues_repo_name: stevengogogo/ProgrammableAging

max_forks_repo_name: stevengogogo/ProgrammableAging