text:
module Conductivity

# This implements different methods to specify conductivity of rocks
#
# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
using ..MaterialParameters: MaterialParamsInfo
import Base.show, GeoParams.param_info

abstract type AbstractConductivity{T} <: AbstractMaterialParam end

export  compute_conductivity,                       # calculation routines
        compute_conductivity!,
        param_info,
        ConstantConductivity,                       # constant
        T_Conductivity_Whittington,                 # T-dependent heat capacity
        T_Conductivity_Whittington_parameterised,   # parameterised T-dependent heat capacity
        TP_Conductivity,                            # TP dependent conductivity
        Set_TP_Conductivity                         # Routine to set pre-defined parameters

include("../Computations.jl")
include("../Utils.jl")

# Constant Conductivity -------------------------------------------------------
"""
    ConstantConductivity(k=3.0W/m/K)

Set a constant conductivity
```math
    k  = cst
```
where ``k`` is the thermal conductivity [``W/m/K``].
"""
@with_kw_noshow struct ConstantConductivity{T,U} <: AbstractConductivity{T}
    k::GeoUnit{T,U}           =   3.0Watt/m/K
end
ConstantConductivity(args...) = ConstantConductivity(convert.(GeoUnit,args)...)

function param_info(s::ConstantConductivity) # info about the struct
    return MaterialParamsInfo(Equation = L"k = cst")
end

# Calculation routine
function (s::ConstantConductivity{_T})(;kwargs...) where _T
    @unpack_val k = s

    return k
end

# (s::ConstantConductivity{_T})(args) where _T = s(;args...)
compute_conductivity(s::ConstantConductivity{_T}; kwargs...) where _T = s()

function (s::ConstantConductivity{_T})(I::Integer...) where _T
    @unpack_val k = s

    return fill(k, I...)
end

compute_conductivity(s::ConstantConductivity{_T}, I::Integer...) where _T = s(I...)

"""
    compute_conductivity(k_array::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::ConstantConductivity) where N

In-place routine to compute constant conductivity
"""
# function compute_conductivity!(k_array::AbstractArray{_T,N}, s::ConstantConductivity{_T}, P::AbstractArray{_T,N}, T::AbstractArray{_T,N}) where {_T,N}
function compute_conductivity!(k_array::AbstractArray{_T,N}, s::ConstantConductivity{_T}; kwargs...) where {_T,N}
    @unpack_val k = s
    for i in eachindex(k_array)
        @inbounds k_array[i] = k
    end
    return nothing
end

# compute_conductivity!(k_array::AbstractArray{_T,N}, s::ConstantConductivity{_T}, args) where {_T,N} = compute_conductivity!(k_array, s; args...)

# Print info
function show(io::IO, g::ConstantConductivity)
    print(io, "Constant conductivity: k=$(g.k.val)")
end
#-------------------------------------------------------------------------

# Temperature dependent conductivity -------------------------------
"""
    T_Conductivity_Whittington()

Sets a temperature-dependent conductivity following the parameterization of *Whittington, A.G., Hofmeister, A.M., Nabelek, P.I., 2009. Temperature-dependent thermal diffusivity of the Earth’s crust and implications for magmatism. Nature 458, 319–321. https://doi.org/10.1038/nature07818.*
Their parameterization is originally given for the thermal diffusivity, together with a parameterization for thermal conductivity, which allows us to compute
```math
    Cp = a + b T - c/T^2
```
```math
    \\kappa = d/T - e \\textrm{ if } T<=846K
```
```math
    \\kappa = f - g*T \\textrm{ if } T>846K
```
```math
    \\rho = 2700 kg/m^3
```
```math
    k = \\kappa \\rho Cp
```

where ``Cp`` is the heat capacity [``J/mol/K``], and ``a,b,c`` are parameters that dependent on the temperature `T`:
- a = 199.50 J/mol/K    if T<= 846 K
- a = 199.50 J/mol/K    if T> 846 K
- b = 0.0857J/mol/K^2   if T<= 846 K
- b = 0.0323J/mol/K^2   if T> 846 K
- c = 5e6J/mol*K        if T<= 846 K
- c = 47.9e-6J/mol*K    if T> 846 K
- d = 576.3m^2/s*K
- e = 0.062m^2/s
- f = 0.732m^2/s
- g = 0.000135m^2/s/K
"""
@with_kw_noshow struct T_Conductivity_Whittington{T,U1,U2,U3,U4,U5,U6,U7,U8,U9} <: AbstractConductivity{T}
    # Note: the resulting curve of k was visually compared with Fig. 2 of the paper
    a0::GeoUnit{T,U1}             =   199.5J/mol/K                # prefactor for low T       (T<= 846 K)
    a1::GeoUnit{T,U1}             =   229.32J/mol/K               # prefactor for high T      (T>  846 K)
    b0::GeoUnit{T,U2}             =   0.0857J/mol/K^2             # linear term for low T     (T<= 846 K)
    b1::GeoUnit{T,U2}             =   0.0323J/mol/K^2             # linear term for high T    (T>  846 K)
    c0::GeoUnit{T,U3}             =   5e6J/mol*K                  # quadratic term for low T  (T<= 846 K)
    c1::GeoUnit{T,U3}             =   47.9e-6J/mol*K              # quadratic term for high T (T>  846 K)
    molmass::GeoUnit{T,U4}        =   0.22178kg/mol               # average molar mass
    Tcutoff::GeoUnit{T,U5}        =   846.0K                      # cutoff temperature
    rho::GeoUnit{T,U6}            =   2700kg/m^3                  # Density they use for an average crust
    d::GeoUnit{T,U7}              =   576.3*1e-6m^2/s*K           # diffusivity parameterization
    e::GeoUnit{T,U8}              =   0.062*1e-6m^2/s             # diffusivity parameterization
    f::GeoUnit{T,U8}              =   0.732*1e-6m^2/s             # diffusivity parameterization
    g::GeoUnit{T,U9}              =   0.000135*1e-6m^2/s/K        # diffusivity parameterization
end
T_Conductivity_Whittington(args...) = T_Conductivity_Whittington(convert.(GeoUnit,args)...)

function param_info(s::T_Conductivity_Whittington) # info about the structwhere {_T,N}
    return MaterialParamsInfo(Equation = L"k = f(T) ")
end

# Calculation routine
function (s::T_Conductivity_Whittington{_T})(; T::_T=zero(_T), kwargs...) where _T
    if T isa Quantity
        @unpack_units a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g = s
    else
        @unpack_val  a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g = s
    end

    if T ≤ Tcutoff
        return (a0 + b0*T - c0/T^2)/molmass * (d/T - e) * rho
    else
        return (a1 + b1*T - c1/T^2)/molmass * (f - g*T) * rho
    end
end

compute_conductivity(s::T_Conductivity_Whittington{_T}; T::_T=zero(_T)) where _T = s(; T=T)

function (s::T_Conductivity_Whittington)(T::AbstractArray; kwargs...)
    if T isa Quantity
        @unpack_units a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g = s
    else
        @unpack_val  a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g = s
    end

    k = similar(T)   #creating an array makes 1 allocation
    inv_molmass = 1/molmass # multiplication is considerably faster than division
    @inbounds for i in eachindex(T)
        if T[i] <= Tcutoff
            a,b,c = a0,b0,c0
            κ     = d/T[i] - e
        else
            a,b,c = a1,b1,c1
            κ     = f - g*T[i]
        end
        cp = (a + b*T[i] - c/T[i]^2)*inv_molmass # conductivity
        k[i] = κ*rho*cp       # compute conductivity from diffusivity
    end

    return k
end

(s::T_Conductivity_Whittington)(T::AbstractArray, args...) = s(T; args...)
compute_conductivity(s::T_Conductivity_Whittington, T::AbstractArray, args...) = s(T; args...)

"""
    compute_conductivity!(k_array::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::T_Conductivity_Whittington) where N

In-place routine to compute temperature-dependent conductivity
"""
function compute_conductivity!(k::AbstractArray{_T,N}, s::T_Conductivity_Whittington{_T}; T::AbstractArray{_T,N}, kwargs...) where {_T,N}
    if T isa Quantity
        @unpack_units a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g = s
    else
        @unpack_val  a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g = s
    end

    inv_molmass = 1/molmass # multiplication is considerably faster than division
    @inbounds for i in eachindex(T)
        if T[i] <= Tcutoff
            k[i] = (a0 + b0*T[i] - c0/T[i]^2)/inv_molmass * (d/T[i] - e) * rho
        else
            k[i] = (a1 + b1*T[i] - c1/T[i]^2)/inv_molmass * (f - g*T[i]) * rho
        end
    end
    return nothing
end

# Print info
function show(io::IO, g::T_Conductivity_Whittington) #info about the struct
    print(io, "T-dependent conductivity following the original Whittington et al. (2009) formulation for average crust). \n");
end
#-------------------------------------------------------------------------

# Temperature dependent conductivity -------------------------------
"""
    T_Conductivity_Whittington_parameterised()

Sets a temperature-dependent conductivity that is  parameterization after *Whittington, et al.  2009*

The original parameterization involves quite a few parameters; this is a polynomial fit that is roughly valid from 0-1000Celcius
```math
    k [W/m/K] = -2 10^{-9} (T-Ts)^3 + 6 10^{-6} (T-Ts)^2 - 0.0062 (T-Ts) + 4
```
```math
    Ts = 273.15 K
```
where `T[K]` is the temperature in Kelvin (or the nondimensional equivalent of it).
"""
@with_kw_noshow struct T_Conductivity_Whittington_parameterised{T,U1,U2,U3,U4,U5} <: AbstractConductivity{T}
    # Note: the resulting curve of k was visually compared with Fig. 2 of the paper
    a::GeoUnit{T,U1}             =  -2e-09Watt/m/K^4                #
    b::GeoUnit{T,U2}             =     6e-06Watt/m/K^3              #
    c::GeoUnit{T,U3}             =   -0.0062Watt/m/K^2              #
    d::GeoUnit{T,U4}             =         4Watt/m/K                #
    Ts::GeoUnit{T,U5}            =   273.15*K                    # Temperature shift [C->K]
end
T_Conductivity_Whittington_parameterised(args...) = T_Conductivity_Whittington_parameterised(convert.(GeoUnit,args)...)

function param_info(s::T_Conductivity_Whittington_parameterised) # info about the struct
    return MaterialParamsInfo(Equation = L"k = a*(T[K]-Ts)^3 + b*(T[K]-Ts)^2  + c*(T[K]-Ts) + d, Ts=273.15")
end

# Calculation routine
function (s::T_Conductivity_Whittington_parameterised{_T})(; T::_T=zero(_T), kwargs...) where _T
    if T isa Quantity
        @unpack_units a,b,c,d,Ts = s
    else
        @unpack_val  a,b,c,d,Ts = s
    end

    T_C  = T - Ts
    return a*T_C^3 + b*T_C^2 + c*T_C + d
end

function (s::T_Conductivity_Whittington_parameterised)(T::AbstractArray; kwargs...)
    if T isa Quantity
        @unpack_units a,b,c,d,Ts = s
    else
        @unpack_val  a,b,c,d,Ts = s
    end

    k = similar(T)   # creating an array makes 1 allocation

    @inbounds for i in eachindex(T)
        # Note: in general, we operate with SI units or the non-dimensional equivalent of that
        # This parameterisation was developed
        T_C  = T[i] - Ts
        k[i] = a*T_C^3 + b*T_C^2 + c*T_C + d
    end

    return k
end

(s::T_Conductivity_Whittington_parameterised)(T::AbstractArray, args) = s(T; args...)
compute_conductivity(s::T_Conductivity_Whittington_parameterised, T::AbstractArray, args)   = s(T, args...)

"""
    compute_conductivity!(k::AbstractArray{_T,N}, s::T_Conductivity_Whittington_parameterised{_T}, P::AbstractArray{_T,N}, T::AbstractArray{_T,N})

In-place routine to compute temperature-dependent conductivity
"""
function compute_conductivity!(k::AbstractArray{_T,N}, s::T_Conductivity_Whittington_parameterised{_T}; T::AbstractArray{_T,N}, kwargs...) where {_T,N}
    if T isa Quantity
        @unpack_units a,b,c,d,Ts = s
    else
        @unpack_val  a,b,c,d,Ts = s
    end

    @inbounds for i in eachindex(T)
        T_C  = T[i] - Ts
        k[i] = a*T_C^3 + b*T_C^2 + c*T_C + d
    end
    return nothing
end

# Print info
function show(io::IO, g::T_Conductivity_Whittington_parameterised) #info about the struct
    print(io, "T-dependent conductivity parameterized from Whittington et al. (2009): k=$(NumValue(g.a))*(T-Ts)^3 + $(NumValue(g.b))*(T-Ts)^2 + $(NumValue(g.c))*(T-Ts) + $(NumValue(g.d)) with T in K and Ts=$(Value(g.Ts))\n");
end
#-------------------------------------------------------------------------

# Temperature (& Pressure) dependent conductivity -------------------------------
"""
    TP_Conductivity()

Sets a temperature (and pressure)-dependent conductivity parameterization as described in Gerya, Numerical Geodynamics (2nd edition, Table 21.2).
The general for

```math
    k = \\left( a_k +  {b_k \\over {T + c_k}} \\right) (1 + d_k P)
```

where ``k`` is the conductivity [``W/K/m``], and ``a_k,b_k,c_k,d_k`` are parameters that dependent on the temperature `T` and pressure `P`:
- ``a_k`` = 1.18Watt/K/m
- ``b_k`` = 474Watt/m
- ``c_k`` = 77K
- ``d_k`` = 0/MPa
"""
@with_kw_noshow struct TP_Conductivity{T,N,U1,U2,U3,U4} <: AbstractConductivity{T}
    Name::NTuple{N,Char}          =   ""                  # The name is encoded as a NTuple{Char} to make it isbits
    a::GeoUnit{T,U1}              =   1.18Watt/K/m        # empirical fitting term
    b::GeoUnit{T,U2}              =   474.0Watt/m         # empirical fitting term
    c::GeoUnit{T,U3}              =   77.0K               # empirical fitting term
    d::GeoUnit{T,U4}              =   0.0/MPa             # empirical fitting term
end
TP_Conductivity(args...) = TP_Conductivity(NTuple{length(args[1]), Char}(collect.(args[1])), convert.(GeoUnit,args[2:end])...)

function param_info(s::TP_Conductivity)
    name = String(collect(s.Name))
    eq = L"k = \left(a_k + {b_k/{T + c_k}} \right)*(1 + d_k*P) "
    if name == ""
        return MaterialParamsInfo(Equation = eq)
    end
    return MaterialParamsInfo(Equation = eq, Comment = TP_Conductivity_info[name][2].Comment)
end

"""
    Set_TP_Conductivity["Name of temperature(-pressure) dependent conductivity"]

This is a dictionary with pre-defined laws:
- "UpperCrust"
- "LowerCrust"
- "OceanicCrust"
- "Mantle"

# Example
```julia
julia> k=Set_TP_Conductivity["Mantle"]
T/P dependent conductivity: k = (0.73 W K⁻¹ m⁻¹ + 1293 W m⁻¹/(T + 77 K))*(1 + 4.0e-5 MPa⁻¹*P)
```

"""
Set_TP_Conductivity(name::String) = TP_Conductivity_info[name][1]

TP_Conductivity_info = Dict([
    ("UpperCrust",
        (TP_Conductivity(Name="UpperCrust", a=0.64Watt/K/m, b=807Watt/m, c=77K, d=0/MPa),
        MaterialParamsInfo(Comment="Sediment/upper crust T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!"))
    )

    ("LowerCrust",
        (TP_Conductivity(Name="LowerCrust", a=1.18Watt/K/m, b=474Watt/m, c=77K, d=0/MPa),
        MaterialParamsInfo(Comment="Lower crust T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!"))
    )

    ("OceanicCrust",
        (TP_Conductivity(Name="OceanicCrust", a=1.18Watt/K/m, b=474Watt/m, c=77K, d=0/MPa),
        MaterialParamsInfo(Comment="Oceanic crust T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!"))
    )

    ("Mantle",
        (TP_Conductivity(Name="Mantle", a=0.73Watt/K/m, b=1293Watt/m, c=77K, d=0.00004/MPa),
        MaterialParamsInfo(Comment="Mantle T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!"))
    )

])

# Calculation routine
function (s::TP_Conductivity{_T})(; P::_T=zero(_T), T::_T=zero(_T), kwargs...) where _T
    if T isa Quantity
        @unpack_units a,b,c,d = s
    else
        @unpack_val  a,b,c,d = s
    end

    if ustrip(d)==0
        return a + b/(T + c)
    else
        return (a + b/(T + c))*(1 + d*P)
    end
end

function (s::TP_Conductivity{_T})(P::AbstractArray{_T,N}, T::AbstractArray{_T,N}; kwargs...) where {_T,N}
    if T isa Quantity
        @unpack_units a,b,c,d = s
    else
        @unpack_val  a,b,c,d = s
    end

    k = similar(T)

    @inbounds if ustrip(d)==0
        for i in eachindex(T)
            k[i] = a + b/(T[i] + c)
        end
    else
        if size(T) != size(P)
            error("Size of P and T arrays should be the same")
        end

        for i in eachindex(T)
            k[i] = (a + b/(T[i] + c))*(1 + d*P[i])
        end
    end

    return k
end

(s::TP_Conductivity)(P::AbstractArray, T::AbstractArray) = s(P, T)
compute_conductivity(s::TP_Conductivity, P::AbstractArray, T::AbstractArray) = s(P, T)

# Calculation routine
function compute_conductivity!(K::AbstractArray{_T, N}, s::TP_Conductivity{_T}; P::AbstractArray{_T, N}, T::AbstractArray{_T, N}, kwargs...) where{_T, N}
    if T isa Quantity
        @unpack_units a,b,c,d = s
    else
        @unpack_val  a,b,c,d = s
    end

    for i in eachindex(T)
        @inbounds if d==0
            K[i] = a + b/(T[i] + c)
        else
            K[i] = (a + b/(T[i] + c))*(one(_T) + d*P[i])
        end
    end

    return nothing
end

# add methods programatically
for myType in (:T_Conductivity_Whittington_parameterised, :T_Conductivity_Whittington, :ConstantConductivity, :TP_Conductivity)
    @eval begin
        (s::$(myType))(args)= s(; args...)
        compute_conductivity(s::$(myType), args) = s(args)
        compute_conductivity!(k::AbstractArray, s::$(myType), args) = compute_conductivity!(k, s; args...)
    end
end

# Print info
function show(io::IO, g::TP_Conductivity)
    if ustrip(Value(g.d))==0
        print(io, "T/P dependent conductivity: Name = $(String(collect(g.Name))), k = $(g.a.val) + $(g.b.val)/(T + $(g.c.val))  \n");
    else
        print(io, "T/P dependent conductivity: Name = $(String(collect(g.Name))), k = ($(g.a.val) + $(g.b.val)/(T + $(g.c.val)))*(1 + $(g.d.val)*P)  \n");
    end
end
#-------------------------------------------------------------------------

# Help info for the calculation routines
"""
    k = compute_conductivity(P, T, s:<AbstractConductivity)

Returns the thermal conductivity `k` at any temperature `T` and pressure `P` using any of the parameterizations implemented.

Currently available:
- ConstantConductivity
- T\\_Conductivity_Whittington
- TP\\_Conductivity

# Example
Using dimensional units
```julia
julia> T  = (250:100:1250)*K;
julia> cp = T_HeatCapacity_Whittington()
julia> Cp = ComputeHeatCapacity(0,T,cp)
```

"""
compute_conductivity()

"""
    k = compute_conductivity(T::Any, s::AbstractConductivity)

Computes conductivity if only temperature (and not pressure) is specified
"""
# compute_conductivity(s::AbstractConductivity, T::AbstractArray{_T}) where _T =  compute_conductivity(s,similar(T),T)

"""
    k = compute_conductivity(s::ConstantConductivity)

Returns conductivity if we are sure that we will only employ constant values throughout the simulation
"""
#compute_conductivity(s::ConstantConductivity) =  compute_conductivity(s,0,0)

# Computational routines needed for computations with the MaterialParams structure
function compute_conductivity(s::AbstractMaterialParamsStruct, args)
    return s.Conductivity[1](args)
end

"""
    compute_conductivity!(K::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},Temp::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of conductivity `K` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `Temp` arrays.
This assumes that the `Phase` of every point is specified as an Integer in the `Phases` array.

_________________________________________________________________________________________________________

compute_conductivity!(k::AbstractArray{T,N}, PhaseRatios::AbstractArray{T, M}, P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)
"""
compute_conductivity(args...) = compute_param(compute_conductivity, args...)
compute_conductivity!(args...) = compute_param!(compute_conductivity, args...)

#= these routines are now computed above
function compute_conductivity!(K::AbstractArray{T, N}, Phases::AbstractArray{<:Integer, N}, P::AbstractArray{T, N},Temp::AbstractArray{T, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where {T<:AbstractFloat,N}

    for i = 1:length(MatParam)

        if !isnothing(MatParam[i].Conductivity)
            # Create views into arrays (so we don't have to allocate)
            ind = Phases .== MatParam[i].Phase;
            K_local     =   view(K   , ind )
            P_local     =   view(P   , ind )
            T_local     =   view(Temp, ind )

            compute_conductivity!(K_local, MatParam[i].Conductivity[1], P_local, T_local )
        end

    end

end

function compute_conductivity!(k::AbstractArray{T, N}, PhaseRatios::AbstractArray{T, M}, P::AbstractArray{T, N},Temp::AbstractArray{T, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where {T<:AbstractFloat, N,M}

    if M!=(N+1)
        error("The PhaseRatios array should have one dimension more than the other arrays")
    end

    k .= 0.0;
    k_local     = zeros(size(k))
    for i = 1:length(MatParam)
        k_local .= 0.0
        Fraction    = selectdim(PhaseRatios,M,i);
        if (maximum(Fraction)>0.0) & (!isnothing(MatParam[i].Conductivity))

            compute_conductivity!(k_local, MatParam[i].Conductivity[1], P, Temp )

            k .= k .+ k_local.*Fraction
        end

    end

end
=#

end


text:
global n_repeats = 10

function upload_tripdata_large_to_s3(bucket_name)
    # # Construct directory of 12 csv files for NYC taxi records from 2012 (total 30 GB) duplicated ncopies times
    # ncopies = 3
    # for month in ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]
    #     for ncopy in 1:ncopies
    #         verify_file_in_s3(
    #             bucket_name,
    #             "tripdata_large/tripdata_$(month)_copy$(ncopy).csv",
    #             "https://s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2012-$(month).csv",
    #         )
    #     end
    # end

    # Construct directory with 4 csv files (~10 GB) replicated n_repeats times
    global n_repeats
    for month in ["01", "02", "03", "04"]
        for ncopy = 1:n_repeats
            verify_file_in_s3(
                bucket_name,
                "tripdata_large/tripdata_$(month)_copy$(ncopy).csv",
                "https://s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2012-$(month).csv",
            )
        end
    end
end

@testset "Stress test of BanyanDataFrames on a large dataset" begin
    run_with_job("Stress test on multi-file dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        upload_tripdata_large_to_s3(bucket)
        tripdata_large = read_csv("s3://$(bucket)/tripdata_large")
        global n_repeats

        # Test Coverage: dropmissing, filter, setindex!, getindex, rename, sort,
        #                length, groupcols, valuecols, groupby, transform, combine,
        #                unique
        # Missing: nonunique, innerjoin, select, subset, allowmissing,
        #          disallowmissing

        # Properties
        @test nrow(tripdata_large) == 61577490 * n_repeats
        @test ncol(tripdata_large) == 18
        @test size(tripdata_large) == (61577490 * n_repeats, 18)
        @test sort(names(tripdata_large)) == [
            "dropoff_datetime",
            "dropoff_latitude",
            "dropoff_longitude",
            "fare_amount",
            "mta_tax",
            "passenger_count",
            "payment_type",
            "pickup_datetime",
            "pickup_latitude",
            "pickup_longitude",
            "rate_code",
            "store_and_fwd_flag",
            "surcharge",
            "tip_amount",
            "tolls_amount",
            "total_amount",
            "trip_distance",
            "vendor_id",
        ]
        @test sort(propertynames(tripdata_large)) == [
            :dropoff_datetime,
            :dropoff_latitude,
            :dropoff_longitude,
            :fare_amount,
            :mta_tax,
            :passenger_count,
            :payment_type,
            :pickup_datetime,
            :pickup_latitude,
            :pickup_longitude,
            :rate_code,
            :store_and_fwd_flag,
            :surcharge,
            :tip_amount,
            :tolls_amount,
            :total_amount,
            :trip_distance,
            :vendor_id,
        ]

        # Dropmissing
        tripdata_large = dropmissing(tripdata_large)
        @test nrow(tripdata_large) == 31210304 * n_repeats

        # Select trips with more than 1 passenger and group by the passenger count.
        # 1 - Compute the fraction of trips for each group that have a trip distance longer
        # than the average trip distance for the group.
        # 2 - Compute the min-max normalized

        tripdata_large = dropmissing(tripdata_large)
        tripdata_large_filtered = filter(row -> (row.passenger_count > 1), tripdata_large)
        res = transform(
            groupby(tripdata_large_filtered, :passenger_count),
            :trip_distance => d -> d .> mean(d),
        )
        res = rename(res, :trip_distance_function => :is_distance_longer_than_mean)
        gdf = groupby(res, :passenger_count)
        res = rename(
            combine(
                groupby(res, :passenger_count),
                nrow,
                :is_distance_longer_than_mean => sum,
            ),
            :nrow => :n_trips,
            :is_distance_longer_than_mean_sum => :total_trip_distance,
        )
        res[:, :fraction] = map(
            (n, s) -> round(s / n, digits = 2),
            res[:, :n_trips],
            res[:, :total_trip_distance],
        )
        res = sort(res, [:fraction, :n_trips])

        res_collected = collect(res[:, [:n_trips, :total_trip_distance, :fraction]])

        @test nrow(tripdata_large_filtered) == 6242979 * n_repeats
        @test nrow(unique(tripdata_large_filtered)) == 6242979
        @test length(gdf) == 6
        @test collect(groupcols(gdf)) == [:passenger_count]
        @test Set(collect(valuecols(gdf))) == Set([
            :vendor_id,
            :pickup_datetime,
            :dropoff_datetime,
            :trip_distance,
            :pickup_longitude,
            :pickup_latitude,
            :rate_code,
            :store_and_fwd_flag,
            :dropoff_longitude,
            :dropoff_latitude,
            :payment_type,
            :fare_amount,
            :surcharge,
            :mta_tax,
            :tip_amount,
            :tolls_amount,
            :total_amount,
            :is_distance_longer_than_mean,
        ])
        @test res_collected[:, :n_trips] ==
              [15368, 1116069, 4510158, 601379, 2, 3] * n_repeats
        @test res_collected[:, :total_trip_distance] ==
              [4239, 325916, 1320006, 178700, 1, 2] * n_repeats
        @test res_collected[:, :fraction] == [0.28, 0.29, 0.3, 0.5, 0.67]

    end
end

ext: jl

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
"""
    interactions(from::MangalNode, ::Colon, query::Pair...)

Returns interactions established *by* the species given as its first argument.
"""
function interactions(from::MangalNode, ::Colon, query::Pair...)
    return interactions(Pair("node_from", string(from.id)), query...)
end

"""
    interactions(::Colon, to::MangalNode, query::Pair...)

Returns interactions established *to* the species given as its second argument.
"""
function interactions(::Colon, to::MangalNode, query::Pair...)
    return interactions(Pair("node_to", string(to.id)), query...)
end

"""
    interactions(with::MangalNode, query::Pair...)

Returns interactions established *around* the species given as its first argument.
"""
function interactions(with::MangalNode, query::Pair...)
    to = interactions(with, :, query...)
    fr = interactions(:, with, query...)
    append!(to, fr)
    return unique(to)
end

"""
    interactions(from::MangalNode, to::MangalNode, query::Pair...)

Returns interactions between two nodes.
"""
function interactions(from::MangalNode, to::MangalNode, query::Pair...)
    return interactions(Pair("node_from", string(from.id)), Pair("node_to", string(to.id)), query...)
end

"""
    interactions(n::MangalNetwork, query::Pair...)

Returns interactions within a network.
"""
function interactions(n::MangalNetwork, query::Pair...)
    return interactions("network_id" => n.id, query...)
end


text:
@testset "FilterModifier" begin

    df = DataFrame(x = [1, 2, 3, 4, 5, 6], y = [6, 5, 4, 3, 2, 1], z = [6, 5, 4, 3, 2, 1])
    conn_df = EasyStream.TablesConnector(df);
    stream = EasyStream.BatchStream(conn_df; batch = 2);

    filter = EasyStream.FilterModifier([:x, :y])
    push!(stream, filter)

    stream_filtered = EasyStream.listen(stream)

    @test (:x in propertynames(stream_filtered)) == true
    @test (:y in propertynames(stream_filtered)) == true
    @test (:z in propertynames(stream_filtered)) == false

    filter = EasyStream.FilterModifier(:x)
    push!(stream, filter)

    stream_filtered = EasyStream.listen(stream)

    @test (:x in propertynames(stream_filtered)) == true
    @test (:y in propertynames(stream_filtered)) == false
    @test (:z in propertynames(stream_filtered)) == false

    @test_logs (:warn, "There are duplicate columns.") EasyStream.FilterModifier([:x, :x, :y])
    @test_logs (:warn, "There are duplicate columns.") EasyStream.FilterModifier(:x, :x, :y)
end


text:
@i @inline function :(*=)(+)(z::Tropical, x::Tropical, y::Tropical)
    @invcheckoff if (content(x) > content(y), ~)
        (z |> content) += (x |> content)
    else
        (z |> content) += (y |> content)
    end
end

@i @inline function (:*=(identity))(x::Tropical, y::Tropical)
    (x |> content) += (y |> content)
end

@i @inline function (:*=(*))(out!::Tropical, x::Tropical, y::Tropical)
    (out! |> content) += (x |> content) + (y |> content)
end

"""
add a number to out!.
"""
@i @inline function i_unsafe_addto(out!::Tropical, x!::Tropical, branch!::Bool)
	@invcheckoff if (content(out!) < content(x!), branch!)
		FLIP(branch!)
		NiLang.SWAP(out!, x!)
	end
end

# branch should be initialized to false.
@i @inline function i_muladd(out!::Tropical, x::Tropical, y::Tropical, branch::Bool)
	x *= y
	@invcheckoff if (out! < x, branch)
		FLIP(branch)
		NiLang.SWAP(out!, x)
	end
end

maxloc(v::AbstractVector) = findmax(v)[2]


text:
"""
    ibcleft(O, AL, C; kwargs...)

Compute left infinite boundary condition, defined by
```
   ┌────────┐   ┌────────┐     ┌─            ┌────────┐    ─┐ P
┌──1   AL   2───1   AL   2─── ┌┘          ───1   AL   2───  └┐  ───(2)
│  └────3───┘   └────3───┘    │   ┌┐         └────3───┘      │
│       3────────────4        │  ─┘│              │          │
│      (       O      )       │    │  ───         │          │
│       1────────────2        │    │              │          │
│  ┌────3───┐   ┌────3───┐    │  ──┴──       ┌────3───┐      │
└──1conj(AL)2───1conj(AL)2─── └┐          ───1conj(AL)2───  ┌┘  ───(1)
   └────────┘   └────────┘     └─            └────────┘    ─┘
```

### Return values:
`ibc, info = ibcleft(O, AL, C)`

### Arguments:
*   `O`: Hamiltonian (only two-site operator is supported)
*   `AL`: left-canonical tensor of mixed-canonical uniform MPS representation
*   `C`: center matrix of mixed-canonical uniform MPS representation

If `O`, `AL` and `C` is `AbstractArray` oblect, nothing have to be done. The following
methods should be defined:

*   `TangentMps.transfer_from_left(X, O, AL)`
*   `TangentMps.transfer_from_left(X, O, (AL, AL))`
*   `TangentMps.similar_normalized_eye(X)`
*   `Base.*(X, Y)`
*   `Base.*(x::Number, X)`
*   `Base.+(X, Y)`
*   `Base.adjoint(X)`
*   `Base.eltype(X)`
*   `Base.similar(X, [T::Type<:Number])`
*   `Base.fill!(X, α::Number)`
*   `Base.copyto!(X, Y)`
*   `LinearAlgebra.mul!(X, Y, α::Number)`
*   `LinearAlgebra.rmul!(X, α::Number)`
*   `LinearAlgebra.axpy!(α::Number, X, Y)`
*   `LinearAlgebra.axpby!(α::Number, X, β::Number, Y)`
*   `LinearAlgebra.tr(X)`
*   `LinearAlgebra.dot(X,Y)`
*   `LinearAlgebra.norm(X)`

where `typeof(X) == typeof(Y) == typeof(C)` is satisfied.

### Keyword arguments:
Keyword arguments are passed to `KrylovKit.linsolve` used in `ibc_left`.
*   `atol::Real`: the requested accuracy, i.e. absolute tolerance, on the norm of the
    residual.
*   `rtol::Real`: the requested accuracy on the norm of the residual, relative to the norm
    of the right hand side `b`.
*   `tol::Real`: the requested accuracy on the norm of the residual which is actually used,
    but which defaults to `max(atol, rtol*norm(b))`. So either `atol` and `rtol` or directly
    use `tol` (in which case the value of `atol` and `rtol` will be ignored).
*   `krylovdim::Integer`: the maximum dimension of the Krylov subspace that will be
    constructed.
*   `maxiter::Integer`: the number of times the Krylov subspace can be rebuilt; see below for
    further details on the algorithms.
"""
function ibcleft(
    O,
    AL,
    C;
    atol = KrylovDefaults.tol,
    rtol = KrylovDefaults.tol,
    krylovdim = KrylovDefaults.krylovdim,
    maxiter = KrylovDefaults.maxiter,
)
    L = transfer_from_left(I, O, (AL, AL))
    eye = similar_normalized_eye(C)
    rhs = L - tr(C' * L * C) * eye
    ibc, inf = linsolve(
        rhs;
        atol = atol,
        rtol = rtol,
        krylovdim = krylovdim,
        maxiter = maxiter,
    ) do Y
        Y - transfer_from_left(Y, AL) + tr(C' * Y * C) * eye
    end
    inf.converged == 1 || @warn "L un-converged"
    ibc, inf
end

"""
    ibcright(O, AR, C; kwargs...)

Compute left infinite boundary condition, defined by
```
        ┌─            ┌────────┐    ─┐ P    ┌────────┐   ┌────────┐
(1)─── ┌┘          ───1   AR   2───  └┐  ───1   AR   2───1   AR   2──┐
       │   ┌┐         └────3───┘      │     └────3───┘   └────3───┘  │
       │  ─┘│              │          │          3────────────4      │
       │    │  ───         │          │         (      O       )     │
       │    │              │          │          1────────────2      │
       │  ──┴──       ┌────3───┐      │     ┌────3───┐   ┌────3───┐  │
(2)─── └┐          ───1conj(AR)2───  ┌┘  ───1conj(AR)2───1conj(AR)2──┘
        └─            └────────┘    ─┘      └────────┘   └────────┘
```

### Return values:
`ibc, info = ibcright(O, AR, C)`

### Arguments:
*   `O`: Hamiltonian (only two-site operator is supported)
*   `AR`: right-canonical tensor of mixed-canonical uniform MPS representation
*   `C`: center matrix of mixed-canonical uniform MPS representation

If `O`, `AR` and `C` is `AbstractArray` oblect, nothing have to be done. The following
methods should be defined:
*   `TangentMps.transfer_from_right(X, O, AR)`
*   `TangentMps.transfer_from_right(X, O, (AR, AR))`
*   `TangentMps.similar_normalized_eye(X)`
*   `Base.*(X, Y)`
*   `Base.*(x::Number, X)`
*   `Base.+(X, Y)`
*   `Base.adjoint(X)`
*   `Base.eltype(X)`
*   `Base.similar(X, [T::Type<:Number])`
*   `Base.fill!(X, α::Number)`
*   `Base.copyto!(X, Y)`
*   `LinearAlgebra.mul!(X, Y, α::Number)`
*   `LinearAlgebra.rmul!(X, α::Number)`
*   `LinearAlgebra.axpy!(α::Number, X, Y)`
*   `LinearAlgebra.axpby!(α::Number, X, β::Number, Y)`
*   `LinearAlgebra.tr(X)`
*   `LinearAlgebra.dot(X,Y)`
*   `LinearAlgebra.norm(X)`

where `typeof(X) == typeof(Y) == typeof(C)` is satisfied.

### Keyword arguments:
Keyword arguments are passed to `KrylovKit.linsolve` used in `ibc_left`.
*   `atol::Real`: the requested accuracy, i.e. absolute tolerance, on the norm of the
    residual.
*   `rtol::Real`: the requested accuracy on the norm of the residual, relative to the norm
    of the right hand side `b`.
*   `tol::Real`: the requested accuracy on the norm of the residual which is actually used,
    but which defaults to `max(atol, rtol*norm(b))`. So either `atol` and `rtol` or directly
    use `tol` (in which case the value of `atol` and `rtol` will be ignored).
*   `krylovdim::Integer`: the maximum dimension of the Krylov subspace that will be
    constructed.
*   `maxiter::Integer`: the number of times the Krylov subspace can be rebuilt; see below for
    further details on the algorithms.
"""
function ibcright(
    O,
    AR,
    C;
    atol = KrylovDefaults.tol,
    rtol = KrylovDefaults.tol,
    krylovdim = KrylovDefaults.krylovdim,
    maxiter = KrylovDefaults.maxiter,
)
    R = transfer_from_right(I, O, (AR, AR))
    eye = similar_normalized_eye(C)
    rhs = R - tr(C * R * C') * eye
    ibc, inf = linsolve(
        rhs;
        atol = atol,
        rtol = rtol,
        krylovdim = krylovdim,
        maxiter = maxiter,
    ) do Y
        Y - transfer_from_right(Y, AR) + tr(C * Y * C') * eye
    end
    inf.converged == 1 || @warn "R un-converged"
    ibc, inf
end

"""
    vumpsstep(O, AL, AR, AC, C; kwargs...)
    vumpsstep(O, AL, AC, C; kwargs...) # with inversion symmetry

### Return values:
    `(AL, AR, AC, C), (ϵL, ϵR), info = vumpsstep(O, AL, AR, AC, C; kwargs...)`
    `(AL, AC, C), (ϵL,), info = vumpsstep(O, AL, AC, C; kwargs...)`
"""
function vumpsstep(
    O,
    AL,
    AR,
    AC,
    C;
    ibcatol = KrylovDefaults.tol,
    ibcrtol = KrylovDefaults.tol,
    ibckrylovdim = KrylovDefaults.krylovdim,
    ibcmaxiter = KrylovDefaults.maxiter,
    eigtol = KrylovDefaults.tol,
    eigkrylovdim = KrylovDefaults.krylovdim,
    eigmaxiter = KrylovDefaults.maxiter,
)

    (L, infL) = ibcleft(
        O,
        AL,
        C;
        atol = ibcatol,
        rtol = ibcrtol,
        krylovdim = ibckrylovdim,
        maxiter = ibcmaxiter,
    )
    (R, infR) = ibcright(
        O,
        AR,
        C;
        atol = ibcatol,
        rtol = ibcrtol,
        krylovdim = ibckrylovdim,
        maxiter = ibcmaxiter,
    )

    (valAC,), (vecAC,), infAC = eigsolve(
        AC,
        1,
        :SR;
        ishermitian = true,
        tol = eigtol,
        krylovdim = eigkrylovdim,
        maxiter = eigmaxiter,
    ) do X
        mul_matrix_from_left(X, L) +
        mul_matrix_from_right(X, R) +
        mul_operator_with_left(X, O, AL) +
        mul_operator_with_right(X, O, AR)
    end
    infAC.converged ≥ 1 || @warn "AC un-converged"
    norm(vecAC) ≈ 1 && (vecAC /= norm(vecAC))

    (valC,), (vecC,), infC = eigsolve(
        C,
        1,
        :SR;
        ishermitian = true,
        tol = eigtol,
        krylovdim = eigkrylovdim,
        maxiter = eigmaxiter,
    ) do X
        Y = mul_matrix_from_left(AR, X)
        transfer_from_right(
            I,
            mul_matrix_from_left(Y, L) +
            mul_matrix_from_right(Y, R) +
            mul_operator_with_left(Y, O, AL) +
            mul_operator_with_right(Y, O, AR),
            AR,
        )
    end
    infC.converged ≥ 1 || @warn "C un-converged"
    norm(vecC) ≈ 1 && (vecC /= norm(vecC))

    AL_, ϵL = al_from_ac_and_c(vecAC, vecC)
    AR_, ϵR = ar_from_ac_and_c(vecAC, vecC)

    (AL_, AR_, vecAC, vecC), (ϵL, ϵR), (infL, infR, infAC, infC)
end
function vumpsstep(
    O,
    AL,
    AC,
    C;
    ibc_atol = KrylovDefaults.tol,
    ibc_rtol = KrylovDefaults.tol,
    ibc_krylovdim = KrylovDefaults.krylovdim,
    ibc_maxiter = KrylovDefaults.maxiter,
    eig_tol = KrylovDefaults.tol,
    eig_krylovdim = KrylovDefaults.krylovdim,
    eig_maxiter = KrylovDefaults.maxiter,
)

    AR = permutedims(AL, (2, 1, 3))

    (L, infL) = ibcleft(
        O,
        AL,
        C;
        atol = ibc_atol,
        rtol = ibc_rtol,
        krylovdim = ibc_krylovdim,
        maxiter = ibc_maxiter,
    )
    R = transpose(L)

    (valAC,), (vecAC,), infAC = eigsolve(
        AC,
        1,
        :SR;
        ishermitian = true,
        tol = eig_tol,
        krylovdim = eig_krylovdim,
        maxiter = eig_maxiter,
    ) do X
        mul_matrix_from_left(X, L) +
        mul_matrix_from_right(X, R) +
        mul_operator_with_left(X, O, AL) +
        mul_operator_with_right(X, O, AR)
    end
    infAC.converged ≥ 1 || @warn "AC un-converged"
    norm(vecAC) ≈ 1 && (vecAC /= norm(vecAC))

    (valC,), (vecC,), infC = eigsolve(
        C,
        1,
        :SR;
        ishermitian = true,
        tol = eig_tol,
        krylovdim = eig_krylovdim,
        maxiter = eig_maxiter,
    ) do X
        Y = mul_matrix_from_left(AR, X)
        transfer_from_right(
            I,
            mul_matrix_from_left(Y, L) +
            mul_matrix_from_right(Y, R) +
            mul_operator_with_left(Y, O, AL) +
            mul_operator_with_right(Y, O, AR),
            AR,
        )
    end
    infC.converged ≥ 1 || @warn "C un-converged"
    norm(vecC) ≈ 1 && (vecC /= norm(vecC))

    AL_, ϵL = al_from_ac_and_c(vecAC, vecC)

    (AL_, vecAC, vecC), (ϵL,), (infL, infAC, infC)
end

"""
    tdvpstep(O, dt, AL, AR, AC, C; kwargs...)
    tdvpstep(O, dt, AL, AC, C; kwargs...) # with inversion symmetry

### Return values:
    `(AL, AR, AC, C), norm, (ϵL, ϵR), info = tdvpstep(O, AL, AR, AC, C; kwargs...)`
    `(AL, AC, C), norm, (ϵL,), info = tdvpstep(O, AL, AC, C; kwargs...)`
where `norm` denotes `|| exp(dt*O)|Ψ(AL, AR, AC, C)⟩ || / || |Ψ(AL, AR, AC, C)⟩ ||`
"""
function tdvpstep(
    O,
    dt,
    AL,
    AR,
    AC,
    C;
    ishermitian,
    ibc_atol = KrylovDefaults.tol,
    ibc_rtol = KrylovDefaults.tol,
    ibc_krylovdim = KrylovDefaults.krylovdim,
    ibc_maxiter = KrylovDefaults.maxiter,
    eig_tol = KrylovDefaults.tol,
    eig_krylovdim = KrylovDefaults.krylovdim,
    eig_maxiter = KrylovDefaults.maxiter,
)
    EO = tr(transfer_from_left(I, O, (AL, AC)))
    O_ = reshape(reshape(O, size(O, 1) * size(O, 2), :) - EO * I, size(O))

    (L, infL) = ibcleft(
        O,
        AL,
        C;
        atol = ibc_atol,
        rtol = ibc_rtol,
        krylovdim = ibc_krylovdim,
        maxiter = ibc_maxiter,
    )
    (R, infR) = ibcright(
        O,
        AR,
        C;
        atol = ibc_atol,
        rtol = ibc_rtol,
        krylovdim = ibc_krylovdim,
        maxiter = ibc_maxiter,
    )
    L = L - tr(transfer_from_left(L, AC)) * I
    R = R - tr(transfer_from_right(R, AC)) * I

    vecAC, infAC = exponentiate(
        dt,
        AC;
        ishermitian = ishermitian,
        tol = eig_tol,
        krylovdim = eig_krylovdim,
        maxiter = eig_maxiter,
    ) do X
        mul_matrix_from_left(X, L) +
        mul_matrix_from_right(X, R) +
        mul_operator_with_left(X, O, AL) +
        mul_operator_with_right(X, O, AR)
    end
    infAC.converged == 1 || @warn "AC un-converged"
    normAC = norm(vecAC)
    vecAC /= normAC
    normAC *= exp(dt * EO)

    vecC, infC = exponentiate(
        dt,
        C;
        ishermitian = ishermitian,
        tol = eig_tol,
        krylovdim = eig_krylovdim,
        maxiter = eig_maxiter,
    ) do X
        Y = mul_matrix_from_left(AR, X)
        transfer_from_right(
            I,
            mul_matrix_from_left(Y, L) +
            mul_matrix_from_right(Y, R) +
            mul_operator_with_left(Y, O, AL) +
            mul_operator_with_right(Y, O, AR),
            AR,
        )
    end
    infC.converged == 1 || @warn "C un-converged"
    normC = norm(vecC)
    vecC /= normC
    normC *= exp(dt * EO)

    AL_, ϵL = al_from_ac_and_c(vecAC, vecC)
    AR_, ϵR = ar_from_ac_and_c(vecAC, vecC)

    (al = AL_, ar = AR_, ac = vecAC, c = vecC),
    (ac = normAC, c = normC),
    (l = ϵL, r = ϵR),
    (l = infL, r = infR, ac = infAC, c = infC)
end
function tdvpstep(
    O,
    dt,
    AL,
    AC,
    C;
    ishermitian,
    ibc_atol = KrylovDefaults.tol,
    ibc_rtol = KrylovDefaults.tol,
    ibc_krylovdim = KrylovDefaults.krylovdim,
    ibc_maxiter = KrylovDefaults.maxiter,
    eig_tol = KrylovDefaults.tol,
    eig_krylovdim = KrylovDefaults.krylovdim,
    eig_maxiter = KrylovDefaults.maxiter,
)
    AR = permutedims(AL, (2, 1, 3))
    EO = tr(transfer_from_left(I, O, (AL, AC)))
    O_ = reshape(reshape(O, size(O, 1) * size(O, 2), :) - EO * I, size(O))

    (L, infL) = ibcleft(
        O_,
        AL,
        C;
        atol = ibc_atol,
        rtol = ibc_rtol,
        krylovdim = ibc_krylovdim,
        maxiter = ibc_maxiter,
    )
    L = L - tr(transfer_from_left(L, AC)) * I
    R = transpose(L)

    vecAC, infAC = exponentiate(
        dt,
        AC;
        ishermitian = ishermitian,
        tol = eig_tol,
        krylovdim = eig_krylovdim,
        maxiter = eig_maxiter,
    ) do X
        mul_matrix_from_left(X, L) +
        mul_matrix_from_right(X, R) +
        mul_operator_with_left(X, O_, AL) +
        mul_operator_with_right(X, O_, AR)
    end
    infAC.converged == 1 || @warn "AC un-converged"
    normAC = norm(vecAC)
    vecAC /= normAC
    normAC *= exp(dt * EO)

    vecC, infC = exponentiate(
        dt,
        C;
        ishermitian = ishermitian,
        tol = eig_tol,
        krylovdim = eig_krylovdim,
        maxiter = eig_maxiter,
    ) do X
        Y = mul_matrix_from_left(AR, X)
        transfer_from_right(
            I,
            mul_matrix_from_left(Y, L) +
            mul_matrix_from_right(Y, R) +
            mul_operator_with_left(Y, O_, AL) +
            mul_operator_with_right(Y, O_, AR),
            AR,
        )
    end
    infC.converged == 1 || @warn "C un-converged"
    normC = norm(vecC)
    vecC /= normC
    normC *= exp(dt * EO)

    AL_, ϵL = al_from_ac_and_c(vecAC, vecC)

    (al = AL_, ac = vecAC, c = vecC),
    (ac = normAC, c = normC),
    (l = ϵL,),
    (l = infL, ac = infAC, c = infC)
end


text:
using VariableSelector
using Test
using DataFrames

@testset "VariableSelector.jl" begin

    con1 = [1, 2, 3, 4, 5, 6,7 ,8, 9, 10]
    con2 = [.17, .23, .36, .46, .58, .67, .77, .84, .92, 1.2]
    cat1 = ["a", "b", "a", "b", "a", "b", "a", "c", "b", "a"]
    cat2 = ["c", "d", "c", "d", "c", "c", "b", "c", "c", "c"]

    #@test names(VariableSelector.variable_importance(df)) == ["feature", "test", "pvalue"]
    @test VariableSelector.gamma(con1, con2)[1] <= 1
    @test VariableSelector.gamma(cat1, cat2)[1] <= 1

end

max_stars_repo_name: ronbactawar/VariableSelection

max_issues_repo_name: ronbactawar/VariableSelection

max_forks_repo_name: ronbactawar/VariableSelection


text:
#Julia Gadfly Scatterplot

using Gadfly
srand(111)
plot(x=rand(7), y=rand(7))

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
using TypeClasses

# Monoid
# ======

TypeClasses.neutral(::Type{String}) = ""
TypeClasses.combine(s1::String, s2::String) = s1 * s2


text:
testOnAllBackends("Vector Algos") do af
	println("\twhere")
	afArr = array(af, [0.0, 0.1, 0.0, 0.3, 0.0, 0.5, 0.0, 0.7, 0.8, 0.0])
	result = where(afArr)
	@test host(result) == [0x00000001,0x00000003,0x00000005,0x00000007,0x00000008]
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
module MPI

using Compat
import Compat.String

using Elemental: ElComm, ElElementType, ElInt, CommWorld, libEl

function __init__()
    # FixMe! The symbol could probably also be missing for other implementations
    #
    # NOTE! I'm using RTLD_GLOBAL here to avoid the OPEN-MPI error described in
    # https://www.open-mpi.org/faq/?category=troubleshooting#missing-symbols
    if Libdl.dlsym_e(Libdl.dlopen(libEl, Libdl.RTLD_GLOBAL), :MPI_Get_library_version) == C_NULL
        const global MPIImpl = :MPICH2
    else
        versionBuffer = Array(UInt8, 2800)
        len = Cint[0]
        err = ccall((:MPI_Get_library_version, libEl), Cint, (Ptr{UInt8}, Ptr{Cint}), versionBuffer, len)
        versionString = String(versionBuffer[1:len[1]-1])
        if ismatch(r"Open MPI", versionString)
            const global MPIImpl = :OpenMPI
        elseif ismatch(r"MPICH", versionString)
            const global MPIImpl = :MPICH3
        else
            error("don't know which MPI implemetation you are using here")
        end
    end
end

function MPIType(t::DataType)
    if MPIImpl == :OpenMPI
        if t == Float64
            return Libdl.dlsym_e(Libdl.dlopen(libEl), :ompi_mpi_double)
        elseif t == Cint
            return Libdl.dlsym_e(Libdl.dlopen(libEl), :ompi_mpi_int)
        elseif t == Clong
            return Libdl.dlsym_e(Libdl.dlopen(libEl), :ompi_mpi_long_int)
        else
            error("data type not defined yet")
        end
    elseif MPIImpl == :MPICH2 || MPIImpl == :MPICH3
        if t == Float64
            return Cint(0x4c00080b)
        elseif t == Cint
            return Cint(0x4c000405)
        elseif t == Clong
            return Cint(0x4c000807)
        else
            error("data type not defined yet")
        end
    else
        error("MPI implementation not covered yet")
    end
end

function MPIOp(f::Function)
    if MPIImpl == :OpenMPI
        if f == (+)
            return Libdl.dlsym_e(Libdl.dlopen(libEl), :ompi_mpi_op_sum)
        else
            error("operation not defined yet")
        end
    elseif MPIImpl == :MPICH2 || MPIImpl == :MPICH3
        if f == (+)
            return Cint(0x58000003)
        else
            error("operation not defined yet")
        end
    else
        error("MPI implementaion no covered yet")
    end
end

function commRank(comm::ElComm)
    n = Ref{Cint}()
    err = ccall((:MPI_Comm_rank, libEl), Cint, (ElComm, Ref{Cint}), comm, n)
    if err != 0
        error("error value was $err")
    end
    return n[]
end

# FixMe! Should be restricted to support element types
function allreduce{T}(sendbuf::Ref{T}, recvbuf::Ref{T}, count::Integer, op::Function, comm::ElComm = CommWorld)
    err = ccall((:MPI_Allreduce, libEl), Cint,
        (Ref{T}, Ref{T}, Cint, ElComm, ElComm, ElComm),
        sendbuf, recvbuf, count, MPIType(T), MPIOp(op), comm)
    if err != 0
        error("error value was $err")
    end
    return recvbuf
end

allreduce{T}(value::T, op::Function, comm::ElComm = CommWorld) = ElInt(allreduce(Ref(value), Ref{T}(), 1, op, comm)[])

end # module


text:
"Abstract vector space type."
abstract type AbstractVectorSpace end

counter = 0

"""
    VectorSpace(dim::Union{Int, Nothing})

A vector space with optional dimension. Attribute id distinguishes between
different vector spaces of the same dimension. User should not need to set id.
"""
struct VectorSpace <: AbstractVectorSpace
    dim::Union{Int, Nothing}
    id::Int
    function VectorSpace(dim::Union{Int, Nothing})
        global counter
        counter += 1
        new(dim, counter)
    end
end

VectorSpace() = VectorSpace(nothing)

"""
    DualVectorSpace(vectorspace::AbstractVectorSpace)

Dual of a given vector space.
"""
struct DualVectorSpace <: AbstractVectorSpace
    vectorspace::AbstractVectorSpace
end

"""
    RnSpace(dim::Int)

Vector space corresponding to R^n for some n.
"""
struct RnSpace <: AbstractVectorSpace
    dim::Int
end

"""
    dual(V::AbstractVectorSpace)

Returns the dual of V.
"""
dual(V::VectorSpace) = DualVectorSpace(V)
dual(Vstar::DualVectorSpace) = Vstar.vectorspace
dual(R::RnSpace) = R
"""
    dim(V::AbstractVectorSpace)

Returns the dimension V.
"""
dim(V::VectorSpace) = V.dim
dim(Vstar::DualVectorSpace) = Vstar.vectorspace.dim
dim(R::RnSpace) = R.dim

"""
    Base.adjoint(V::AbstractVectorSpace)

Returns the [`dual(V)`](@ref).
"""
Base.adjoint(V::AbstractVectorSpace) = dual(V)

Base.show(io::IO, Vstar::DualVectorSpace) = print(io, Vstar', "*")

Base.show(io::IO, R::RnSpace) = print(io,  "R^",dim(R))

function Base.show(io::IO, V::VectorSpace)
    if dim(V) === nothing
        print(io, "V", "_", V.id)
    else
        print(io, "V", dim(V), "_", V.id)
    end
end

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
# Load X and y variable
using SGDBabysitter
using Pkg
using HDF5
using Plots

data=h5open("test\\mnist.hdf5","r")
xtrain=transpose(read(data,"x_train"))
ytrain=read(data,"t_train")
xvalid=transpose(read(data,"x_valid"))
yvalid=read(data,"t_valid")

n = size(xtrain,1)
xtrain = [ones(n,1) xtrain]
d = size(xtrain,2)

# Choose network structure and randomly initialize weights
include("NeuralNet.jl")
nHidden = [3]
nParams = NeuralNet_nParams(d,nHidden)
w = randn(nParams,1)

# Train with stochastic gradient
maxIter = 10000
j=1
checks=[]
valid=[]
stepSize = 1e-4
W= randn(nParams,1)

for t in 1:maxIter
	global W, j
	# The stochastic gradient update:
	i = rand(1:n)
	f,g = NeuralNet_backprop(W, xtrain[i,:], ytrain[i], nHidden)
	W = W - stepSize*g

	# Every few iterations, plot the data/model:
	if (mod(t-1,round(maxIter/50)) == 0)
		print("Training iteration = $(t-1)",)

		yhat = NeuralNet_predict(W,[ones(size(xvalid,1)) xvalid],nHidden)

		push!(valid,sum((yhat-yvalid).^2)/length(yvalid))
		push!(checks,j)
		j=j+1
	end
end
plot(checks,valid)


text:
export sigma_LorentzBerthelot_mod
function sigma_LorentzBerthelot_mod(sigma::ClapeyronParam,zeta::PairParam)
    sigma = PairParam(sigma)
    σ = sigma.values
    zeta_ = zeta.values
    for i ∈ 1:length(sigma.components), j ∈ 1:length(sigma.components)
        if sigma.ismissingvalues[i,j]
            σ[i,j] = (1-zeta_[i,j])*(σ[i,i] + σ[j,j]) / 2
        end
    end
    return sigma
end

export sigma_LorentzBerthelot
function sigma_LorentzBerthelot(sigma::ClapeyronParam)
    sigma = PairParam(sigma)
    σ = sigma.values
    for i ∈ 1:length(sigma.components), j ∈ 1:length(sigma.components)
        if sigma.ismissingvalues[i,j]
            σ[i,j] = (σ[i,i] + σ[j,j]) / 2
        end
    end
    sigma.ismissingvalues .= false
    return sigma
end

export epsilon_LorentzBerthelot
function epsilon_LorentzBerthelot(epsilon::ClapeyronParam, k::PairParam)
    epsilon = PairParam(epsilon)
    ϵ = epsilon.values
    k_ = k.values
    for i ∈ 1:length(epsilon.components), j ∈ 1:length(epsilon.components)
        if epsilon.ismissingvalues[i,j]
            ϵ[i,j] = (1 - k_[i,j]) * sqrt(ϵ[i,i] * ϵ[j,j])
        end
    end
    epsilon.ismissingvalues .= false
    return epsilon
end

export epsilon_HudsenMcCoubrey
function epsilon_HudsenMcCoubrey(epsilon::ClapeyronParam, sigma::PairParam)
    epsilon = PairParam(epsilon)
    ϵ = epsilon.values
    σ = sigma.values
    for i ∈ 1:length(epsilon.components), j ∈ 1:length(epsilon.components)
        if epsilon.ismissingvalues[i,j]
            ϵ[i,j] = (σ[i,i]^3 * σ[j,j]^3) / σ[i,j]^6 * sqrt(ϵ[i,i] * ϵ[j,j])
        end
    end
    epsilon.ismissingvalues .= false
    return epsilon
end

export lambda_LorentzBerthelot
function lambda_LorentzBerthelot(lambda::ClapeyronParam)
    lambda = PairParam(lambda)
    λ = lambda.values
    for i ∈ 1:length(lambda.components), j ∈ 1:length(lambda.components)
        if lambda.ismissingvalues[i,j]
            λ[i,j] = 3 + sqrt((λ[i,i] - 3) * (λ[j,j] - 3))
        end
    end
    lambda.ismissingvalues .= false
    return lambda
end

export lambda_squarewell
function lambda_squarewell(lambda::ClapeyronParam, sigma::PairParam)
    lambda = PairParam(lambda)
    λ = lambda.values
    σ = sigma.values
    for i ∈ 1:length(lambda.components), j ∈ 1:length(lambda.components)
        if lambda.ismissingvalues[i,j]
            λ[i,j] = (λ[i,i]*σ[i,i] + λ[j,j]*σ[j,j]) / (σ[i,i] + σ[j,j])
        end
    end
    lambda.ismissingvalues .= false
    return lambda
end

"""
    mixing_rule_quad(op, x, p)

returns an efficient implementation of:
sum(x[i]*x[j]*op(p[i],p[j]) for i in @comps for j in @comps)`
where `op` is a function of two arguments that satisfies op(p[i],p[i]) = p[i]

# example
```julia-repl
julia> prop = [200,600];
julia> x0 = [0.3,0.7];
julia> propmix1 = Clapeyron.mixing_rule_quad((x,y)->0.5*(x+y),x0,prop);
julia> propmix2 = sum(x0[i]*x0[j]*(prop[i]+prop[j])*0.5 for i in 1:2 for j in 1:2);
julia> propmix1 ≈ propmix2;
true
```
"""
function mixing_rule_quad(op, x, p)
    N = length(x)
    @boundscheck checkbounds(x, N)
    @inbounds begin
        res1 = zero(eltype(x))
        for i = 1:N
            res1 += p[i] * x[i]^2
            for j = 1:i - 1
                res1 += 2 * x[i] * x[j] * op(p[i], p[j])
            end
        end
    end
    return res1
end


text:
abstract AbstractFormatter

immutable PlaintextFormatter <: AbstractFormatter end
immutable MarkdownFormatter  <: AbstractFormatter end

immutable Format{F <: AbstractFormatter} end

"""
Parsing hook for docstring parsing.

Example:

    import Docile.Formats, Markdown
    const (fmt, md) = (Docile.Formats, Markdown)

    fmt.parsedocs(::fmt.Format{fmt.MarkdownFormatter}, raw, mod, obj) = md.parse(raw)

When registering a package the format is then provided to `PackageData`.

Arguments:

- `raw`: Unparsed docstring extracted from source code.
- `mod`: Module where the docstring and object are defined.
- `obj`: Object that the docstring is associated with.

"""
parsedocs(::Format, raw, mod, obj) = raw