text:
## Metsa-Oy Forest Products Supply Chain
## Edward J. Xu <edxu96@outlook.com>
## Feb 1st, 2020

function optim_1(alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j, c_j,
  r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik, gamma_jk,
  delta_jk, gamma_paper_k, delta_paper_k, p_fuel)

  mod = Model(with_optimizer(Gurobi.Optimizer, Presolve=0, OutputFlag=0))

  @variable(mod, h_t[1:6]>=0, integer=true)
  @variable(mod, y_i[1:5]>=0)
  @variable(mod, y_j[1:2]>=0)
  @variable(mod, y_paper>=0)
  @variable(mod, z_ik[1:5, 1:4]>=0)
  @variable(mod, z_jk[1:2, 1:4]>=0)
  @variable(mod, z_paper_k[1:4]>=0)

  i_t1 = [[1], [2, 4], [3, 5]]
  i_t2 = [[1], [2, 4], [3, 5]]
  j_t2 = [[1], [], [2]]

  @objective(mod, Max,
    -sum(alpha_t[t] * h_t[t] + beta_t[t] * h_t[t]^2 for t = 1:6) +
    -sum(c_i[i] * y_i[i] for i=1:5) +
    -sum(c_j[j] * y_j[j] for j=1:2) - c_paper * y_paper +
    sum(alpha_t[t] * (
      h_t[t] - sum(a_i[i] * y_i[i] for i in i_t1[t])
      ) for t=1:3) +
    sum(alpha_t[t] * (
      h_t[t] + sum(b_i[i] * y_i[i] for i in i_t2[t-3]) -
      sum(a_j[j] * y_j[j] for j in j_t2[t-3]) -
        a_paper_t[t-3] * y_paper) for t=4:6) +
    sum(e_i[i] * y_i[i] for i=1:5) * p_fuel +
    sum(gamma_ik[i, k] * z_ik[i, k] - delta_ik[i, k] * z_ik[i, k]^2
      for i=1:5, k=1:4) +
    sum(gamma_jk[j, k] * z_jk[j, k] - delta_jk[j, k] * z_jk[j, k]^2
      for j=1:2, k=1:4) +
    sum(gamma_paper_k[k] * z_paper_k[k] -
      delta_paper_k[k] * z_paper_k[k]^2
      for k=1:4)
    )

  @constraint(mod, vec_cons_1[t=1:3],
    h_t[t] >= sum(a_i[i] * y_i[i] for i in i_t1[t]))
  @constraint(mod, vec_cons_2[t=4:6],
    h_t[t] + sum(b_i[i] * y_i[i] for i in i_t2[t-3]) >=
    sum(a_j[j] * y_j[j] for j in j_t2[t-3]) + a_paper_t[t-3] * y_paper)
  @constraint(mod, vec_cons_3[j=1:2],
    y_j[j] >= b_paper_j[j] * y_paper)
  @constraint(mod, vec_cons_4[i=1:5],
    y_i[i] >= sum(z_ik[i, k] for k=1:4))
  @constraint(mod, vec_cons_5[j=1:2],
    y_j[j] - b_paper_j[j] * y_paper >= sum(z_jk[j, k] for k=1:4))
  @constraint(mod, y_paper >= sum(z_paper_k[k] for k=1:4))
  @constraint(mod, sum(y_i[i] for i=1:3) <= r_saw)
  @constraint(mod, sum(y_i[i] for i=4:5) <= r_plywood)
  @constraint(mod, vec_cons_6[j=1:2], y_j[j] <= r_j[j])
  @constraint(mod, y_paper <= r_paper)

  optimize!(mod)

  result_h_t = value_vec(h_t)
  result_y_i = value_vec(y_i)
  result_y_j = value_vec(y_j)
  result_y_paper = value(y_paper)
  result_z_ik = value_mat(z_ik)
  result_z_jk = value_mat(z_jk)
  result_z_paper_k = value_vec(z_paper_k)
  result_obj = objective_value(mod)

  println(result_h_t)
  println(result_y_i)
  println(result_y_j)
  println(result_y_paper)
  println(result_z_ik)
  println(result_z_jk)
  println(result_z_paper_k)
  println(result_obj)

  # return result_h_t, result_y_i, result_y_j, result_y_paper, result_z_ik,
  #   result_z_jk, result_z_paper_k, result_obj
end


text:
module Braking
    export dist, speed

    const EARTH_GRAVITY = 9.81

    function dist(v::Number, μ::Number)::Real
        # convert from km/h to m/s
        v /= 3.6
        # Reaction time: 1
        v + v^2/2/μ/EARTH_GRAVITY
    end

    # Reaction time: 1
    function speed(d::Number, μ::Number)::Real
        μg = μ * EARTH_GRAVITY
        v = sqrt( μg^2 + 2d*μg ) - μg
        # Convert m/s in km/h
        3.6v
    end
end

ext: jl

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module TestAir

using Test
using LinearAlgebra
using Canopy.Air

@time @testset "Canopy.Air" begin
    # TODO:
end

end  # module


text:
@with_kw struct RECFAST{T, AB<:AbstractBackground{T}} <: IonizationIntegrator @deftype T
    bg::AB  # a RECFAST has an associated background evolution
    H0_natural_unit_conversion = ustrip(u"s", unnatural(u"s", 1u"eV^-1"))

    bigH = 100.0e3 / (1e6 * 3.0856775807e16)	 # H₀ in s-1
    C  = 2.99792458e8  # Fundamental constants in SI units
    k_B = 1.380658e-23
    h_P = 6.6260755e-34
    m_e = 9.1093897e-31
    m_H = 1.673575e-27  #	av. H atom
    # note: neglecting deuterium, making an O(e-5) effect
    not4 = 3.9715e0  # mass He/H atom  ("not4" pointed out by Gary Steigman)
    sigma = 6.6524616e-29
    a = 7.565914e-16
    G = 6.6742e-11 	# new value

    Lambda = 8.2245809e0
    Lambda_He = 51.3e0              # new value from Dalgarno
    L_H_ion = 1.096787737e7         # level for H ion. (in m^-1)
    L_H_alpha = 8.225916453e6       # averaged over 2 levels
    L_He1_ion = 1.98310772e7        # from Drake (1993)
    L_He2_ion = 4.389088863e7       # from JPhysChemRefData (1987)
    L_He_2s	= 1.66277434e7          # from Drake (1993)
    L_He_2p	= 1.71134891e7          # from Drake (1993)
    # C	2 photon rates and atomic levels in SI units

    A2P_s = 1.798287e9              # Morton, Wu & Drake (2006)
    A2P_t = 177.58e0                # Lach & Pachuski (2001)
    L_He_2Pt = 1.690871466e7        # Drake & Morton (2007)
    L_He_2St = 1.5985597526e7       # Drake & Morton (2007)
    L_He2St_ion = 3.8454693845e6    # Drake & Morton (2007)
    sigma_He_2Ps = 1.436289e-22     # Hummer & Storey (1998)
    sigma_He_2Pt = 1.484872e-22     # Hummer & Storey (1998)
    # C	Atomic data for HeI

    AGauss1	= -0.14e0               # Amplitude of 1st Gaussian
    AGauss2 = 0.079e0               # Amplitude of 2nd Gaussian
    zGauss1 = 7.28e0                # ln(1+z) of 1st Gaussian
    zGauss2 = 6.73e0                # ln(1+z) of 2nd Gaussian
    wGauss1 = 0.18e0                # Width of 1st Gaussian
    wGauss2 = 0.33e0                # Width of 2nd Gaussian
    # Gaussian fits for extra H physics (fit by Adam Moss, modified by Antony Lewis)

    # the Pequignot, Petitjean & Boisson fitting parameters for Hydrogen
	a_PPB = 4.309
	b_PPB = -0.6166
	c_PPB = 0.6703
	d_PPB = 0.5300
    # the Verner and Ferland type fitting parameters for Helium
    # fixed to match those in the SSS papers, and now correct
	a_VF = 10^(-16.744)
	b_VF = 0.711
	T_0 = 10^(0.477121)	#!3K
	T_1 = 10^(5.114)
    # fitting parameters for HeI triplets
    # (matches Hummer's table with <1% error for 10^2.8 < T/K < 10^4)
	a_trip = 10^(-16.306)
	b_trip = 0.761

    # Set up some constants so they don't have to be calculated later
    Lalpha = 1/L_H_alpha
    Lalpha_He = 1/L_He_2p
    DeltaB = h_P*C*(L_H_ion-L_H_alpha)
    CDB = DeltaB/k_B
    DeltaB_He = h_P*C*(L_He1_ion-L_He_2s)	# 2s, not 2p
    CDB_He = DeltaB_He/k_B
    CB1 = h_P*C*L_H_ion/k_B
    CB1_He1 = h_P*C*L_He1_ion/k_B	# ionization for HeI
    CB1_He2 = h_P*C*L_He2_ion/k_B	# ionization for HeII
    CR = 2π * (m_e/h_P)*(k_B/h_P)
    CK = Lalpha^3/(8π)
    CK_He = Lalpha_He^3/(8π)
    CL = C*h_P/(k_B*Lalpha)
    CL_He = C*h_P/(k_B/L_He_2s)	# comes from det.bal. of 2s-1s
    CT = (8/3)*(sigma/(m_e*C))*a
    Bfact = h_P*C*(L_He_2p-L_He_2s)/k_B

    # Matter departs from radiation when t(Th) > H_frac * t(H)
    H_frac = 1e-3  # choose some safely small number

    # switches
    Hswitch::Int64 = 1
    Heswitch::Int64 = 6

    # Cosmology
    Yp = 0.24
    OmegaB = 0.046  # TODO: should replace during GREAT GENERALIZATION
    HO =  bg.H₀ / H0_natural_unit_conversion
    Tnow = 2.725

    # sort out the helium abundance parameters
    mu_H = 1 / (1 - Yp)			 # Mass per H atom
    mu_T = not4/(not4-(not4-1)*Yp)	 # Mass per atom
    fHe = Yp/(not4*(1 - Yp))		# n_He_tot / n_H_tot

    Nnow = 3 * HO * HO * OmegaB / (8π * G * mu_H * m_H)  # TODO: should replace during GREAT GENERALIZATION
    # fnu = (21/8)*(4/11)^(4/3)
    # (this is explictly for 3 massless neutrinos - change if N_nu.ne.3)  # this is only for H(z) and ∂H/∂z
    # z_eq = (3 * (HO*C)^2 / (8π * G * a * (1+fnu)*Tnow^4))*OmegaT - 1

    fu = (Hswitch == 0) ? 1.14 : 1.125
    b_He = 0.86  # Set the He fudge factor
    tol = 1e-6
end

function recfast_init(𝕣::RECFAST, z)
    if z > 8000.
        x_H0 = 1.
        x_He0 = 1.
        x0 = 1. + 2 * 𝕣.fHe
    elseif z > 3500.
        x_H0 = 1.
        x_He0 = 1.
        rhs = exp( 1.5 * log(𝕣.CR*𝕣.Tnow/(1 + z)) - 𝕣.CB1_He2/(𝕣.Tnow*(1 + z)) ) / 𝕣.Nnow
	    rhs = rhs * 1.  # ratio of g's is 1 for He++ <-> He+
	    x0 = 0.5 * ( sqrt( (rhs - 1 - 𝕣.fHe)^2 + 4 * (1 + 2 * 𝕣.fHe) * rhs) - (rhs - 1 - 𝕣.fHe) )
    elseif z > 2000.
	    x_H0 = 1.
	    rhs = exp( 1.5 * log(𝕣.CR * 𝕣.Tnow / (1 + z)) - 𝕣.CB1_He1/(𝕣.Tnow*(1 + z)) ) / 𝕣.Nnow
	    rhs = 4rhs    # ratio of g's is 4 for He+ <-> He0
	    x_He0 = 0.5 * ( sqrt( (rhs-1)^2 + 4*(1 + 𝕣.fHe)*rhs) - (rhs-1))
	    x0 = x_He0
	    x_He0 = (x0 - 1.)/𝕣.fHe
    else
	    rhs = exp( 1.5 * log(𝕣.CR*𝕣.Tnow/(1 + z)) - 𝕣.CB1/(𝕣.Tnow*(1 + z)) ) / 𝕣.Nnow
	    x_H0 = 0.5 * (sqrt( rhs^2 + 4 * rhs ) - rhs )
	    x_He0 = 0.
	    x0 = x_H0
    end

    return x_H0, x_He0, x0
end

function ion_recfast!(f, y, 𝕣::RECFAST, z)

	x_H = y[1]
	x_He = y[2]
	x = x_H + 𝕣.fHe * x_He
	Tmat = y[3]

	n = 𝕣.Nnow * (1+z)^3
	n_He = 𝕣.fHe * 𝕣.Nnow * (1+z)^3
	Trad = 𝕣.Tnow * (1+z)

    a = 1 / (1+z)  # scale factor
    x_a = a2x(a)
	Hz = 𝕣.bg.ℋ(x_a) / a / 𝕣.H0_natural_unit_conversion
	dHdz = (-𝕣.bg.ℋ′(x_a) + 𝕣.bg.ℋ(x_a)) / 𝕣.H0_natural_unit_conversion
	# Hz = 𝕣.HO * sqrt((1+z)^4/(1+𝕣.z_eq)*𝕣.OmegaT + 𝕣.OmegaT*(1+z)^3 + 𝕣.OmegaK*(1+z)^2 + 𝕣.OmegaL)
	# dHdz = (𝕣.HO^2 /2/Hz)*(4*(1+z)^3/(1+𝕣.z_eq)*𝕣.OmegaT + 3*𝕣.OmegaT*(1+z)^2 + 2*𝕣.OmegaK*(1+z))

    # Get the radiative rates using PPQ fit (identical to Hummer's table)
	Rdown=1e-19*𝕣.a_PPB*(Tmat/1e4)^𝕣.b_PPB/(1. + 𝕣.c_PPB*(Tmat/1e4)^𝕣.d_PPB)
	Rup = Rdown * (𝕣.CR*Tmat)^(1.5)*exp(-𝕣.CDB/Tmat)

    # calculate He using a fit to a Verner & Ferland type formula
	sq_0 = sqrt(Tmat/𝕣.T_0)
	sq_1 = sqrt(Tmat/𝕣.T_1)
    # typo here corrected by Wayne Hu and Savita Gahlaut
	Rdown_He = 𝕣.a_VF/(sq_0*(1+sq_0)^(1-𝕣.b_VF))
	Rdown_He = Rdown_He/(1+sq_1)^(1+𝕣.b_VF)
	Rup_He = Rdown_He*(𝕣.CR*Tmat)^(1.5)*exp(-𝕣.CDB_He/Tmat)
	Rup_He = 4. * Rup_He # statistical weights factor for HeI
    # Avoid overflow (pointed out by Jacques Roland)
	if((𝕣.Bfact/Tmat) > 680.)
	  He_Boltz = exp(680.)
	else
	  He_Boltz = exp(𝕣.Bfact/Tmat)
	end

    # now deal with H and its fudges
	if (𝕣.Hswitch == 0)
	    K = 𝕣.CK / Hz # !Peebles coefficient K=lambda_a^3/8piH
	else
        # fit a double Gaussian correction function
        K = 𝕣.CK / Hz*(1.0
            + 𝕣.AGauss1*exp(-((log(1+z)-𝕣.zGauss1)/𝕣.wGauss1)^2)
            + 𝕣.AGauss2*exp(-((log(1+z)-𝕣.zGauss2)/𝕣.wGauss2)^2))
	end

    # add the HeI part, using same T_0 and T_1 values
	Rdown_trip = 𝕣.a_trip/(sq_0*(1+sq_0)^(1-𝕣.b_trip))
	Rdown_trip = Rdown_trip/((1+sq_1)^(1+𝕣.b_trip))
	Rup_trip = Rdown_trip*exp(-𝕣.h_P*𝕣.C*𝕣.L_He2St_ion/(𝕣.k_B*Tmat))
	Rup_trip = Rup_trip*((𝕣.CR*Tmat)^1.5)*(4/3)
    # last factor here is the statistical weight

    # try to avoid "NaN" when x_He gets too small
	if ((x_He < 5.e-9) || (x_He > 0.980))
        Heflag = 0
	else
	    Heflag = 𝕣.Heswitch
	end
	if (Heflag == 0)  # use Peebles coeff. for He
	    K_He = 𝕣.CK_He/Hz
	else # for Heflag>0 		!use Sobolev escape probability
        tauHe_s = 𝕣.A2P_s*𝕣.CK_He*3*n_He*(1-x_He)/Hz
        pHe_s = (1 - exp(-tauHe_s))/tauHe_s
        K_He = 1 / (𝕣.A2P_s*pHe_s*3*n_He*(1-x_He))
        # smoother criterion here from Antony Lewis & Chad Fendt
	    if (((Heflag == 2) || (Heflag >= 5)) && (x_H < 0.9999999))
            # use fitting formula for continuum opacity of H
            # first get the Doppler width parameter
            Doppler = 2*𝕣.k_B*Tmat/(𝕣.m_H*𝕣.not4*𝕣.C*𝕣.C)
            Doppler = 𝕣.C*𝕣.L_He_2p*sqrt(Doppler)
            gamma_2Ps = 3*𝕣.A2P_s*𝕣.fHe*(1-x_He)*𝕣.C*𝕣.C /(
                sqrt(π)*𝕣.sigma_He_2Ps*8π*Doppler*(1-x_H)) /((𝕣.C*𝕣.L_He_2p)^2)
            pb = 0.36 # value from KIV (2007)
            qb = 𝕣.b_He
            # calculate AHcon, the value of A*p_(con,H) for H continuum opacity
            AHcon = 𝕣.A2P_s/(1+pb*(gamma_2Ps^qb))
            K_He = 1/((𝕣.A2P_s*pHe_s+AHcon)*3*n_He*(1-x_He))
	    end
	    if (Heflag >= 3) # include triplet effects
            tauHe_t = 𝕣.A2P_t*n_He*(1. - x_He)*3
            tauHe_t = tauHe_t /(8π*Hz*𝕣.L_He_2Pt^3)
            pHe_t = (1 - exp(-tauHe_t))/tauHe_t
            CL_PSt = 𝕣.h_P*𝕣.C*(𝕣.L_He_2Pt - 𝕣.L_He_2St)/𝕣.k_B
            if ((Heflag == 3) || (Heflag == 5) || (x_H > 0.99999))
                # no H cont. effect
                CfHe_t = 𝕣.A2P_t*pHe_t*exp(-CL_PSt/Tmat)
                CfHe_t = CfHe_t/(Rup_trip+CfHe_t) # "C" factor for triplets
            else # include H cont. effect
                Doppler = 2*𝕣.k_B*Tmat/(𝕣.m_H*𝕣.not4*𝕣.C*𝕣.C)
                Doppler = 𝕣.C*𝕣.L_He_2Pt*sqrt(Doppler)
                gamma_2Pt = (3*𝕣.A2P_t*𝕣.fHe*(1-x_He)*𝕣.C*𝕣.C
                    /(sqrt(π)*𝕣.sigma_He_2Pt*8π*Doppler*(1-x_H))
                    /((𝕣.C*𝕣.L_He_2Pt)^2))
                # use the fitting parameters from KIV (2007) in this case
                pb = 0.66
                qb = 0.9
                AHcon = 𝕣.A2P_t/(1+pb*gamma_2Pt^qb)/3
                CfHe_t = (𝕣.A2P_t*pHe_t+AHcon)*exp(-CL_PSt/Tmat)
                CfHe_t = CfHe_t/(Rup_trip+CfHe_t)  # "C" factor for triplets
            end
	    end
	end

    # Estimates of Thomson scattering time and Hubble time
	timeTh=(1/(𝕣.CT*Trad^4))*(1+x+𝕣.fHe)/x	#!Thomson time
	timeH=2/(3*𝕣.HO*(1+z)^1.5)		#!Hubble time

    # calculate the derivatives
    # turn on H only for x_H<0.99, and use Saha derivative for 0.98<x_H<0.99
    # (clunky, but seems to work)
	if (x_H > 0.99)  # don't change at all
		f[1] = 0.
    # else if ((x_H.gt.0.98d0).and.(Heflag.eq.0)) then	!don't modify
	elseif (x_H > 0.985)  # !use Saha rate for Hydrogen
		f[1] = (x*x_H*n*Rdown - Rup*(1-x_H)*exp(-𝕣.CL/Tmat))/(Hz*(1+z))
        # for interest, calculate the correction factor compared to Saha
        # (without the fudge)
		factor=(1 + K*𝕣.Lambda*n*(1-x_H))/(Hz*(1+z)*(1+K*𝕣.Lambda*n*(1-x)+K*Rup*n*(1-x)))
    else  #!use full rate for H
		f[1] = (((x*x_H*n*Rdown - Rup*(1.0-x_H)*exp(-𝕣.CL/Tmat))
			*(1.0 + K*𝕣.Lambda*n*(1.0-x_H)))
		    /(Hz*(1.0+z)*(1.0/𝕣.fu+K*𝕣.Lambda*n*(1.0-x_H)/𝕣.fu
		    +K*Rup*n*(1.0-x_H))))
	end
    # turn off the He once it is small
	if (x_He < 1e-15)
		f[2] = 0.
	else
		f[2] = (((x*x_He*n*Rdown_He - Rup_He*(1-x_He)*exp(-𝕣.CL_He/Tmat))
            *(1+ K_He*𝕣.Lambda_He*n_He*(1-x_He)*He_Boltz))
            / (Hz*(1+z)
            * (1 + K_He*(𝕣.Lambda_He+Rup_He)*n_He*(1-x_He)*He_Boltz)))
        # Modification to HeI recombination including channel via triplets
	    if (Heflag >= 3)
		    f[2] = f[2] + (x*x_He*n*Rdown_trip
                - (1-x_He)*3*Rup_trip*exp(-𝕣.h_P*𝕣.C*𝕣.L_He_2St/(𝕣.k_B*Tmat))
                ) * CfHe_t/(Hz*(1+z))
	    end
	end

    # follow the matter temperature once it has a chance of diverging

	if (timeTh < 𝕣.H_frac*timeH)
    # f(3)=Tmat/(1.d0+z)	!Tmat follows Trad
    # additional term to smooth transition to Tmat evolution,
    # (suggested by Adam Moss)
		epsilon = Hz*(1+x+𝕣.fHe)/(𝕣.CT*Trad^3*x)
		f[3] = 𝕣.Tnow + epsilon*((1+𝕣.fHe)/(1+𝕣.fHe+x))*(
            (f[1]+𝕣.fHe*f[2])/x) - epsilon* dHdz/Hz + 3*epsilon/(1+z)
	else
		f[3] = 𝕣.CT * (Trad^4) * x / (1+x+𝕣.fHe)* (Tmat-Trad) / (Hz*(1+z)) + 2*Tmat/(1+z)
	end

	return
end

function recfast_xe(𝕣::RECFAST{T};
        Hswitch::Int=1, Heswitch::Int=6, Nz::Int=1000, zinitial=10000., zfinal=0.,
        alg=Tsit5()) where T

    z = zinitial
    n = 𝕣.Nnow * (1 + z)^3
    y = zeros(3)  # array is x_H, x_He, Tmat (Hydrogen ionization, Helium ionization, matter temperature)

    y[3] = 𝕣.Tnow * (1 + z)
    Tmat = y[3]

    x_H0, x_He0, x0 = recfast_init(𝕣, z)
    y[1] = x_H0
    y[2] = x_He0

    out_xe = zeros(T, Nz)
    out_Tmat = zeros(T, Nz)

    for i in 1:Nz
        # calculate the start and end redshift for the interval at each z
        # or just at each z
	    zstart = zinitial + float(i-1)*(zfinal-zinitial)/float(Nz)
	    zend   = zinitial + float(i)*(zfinal-zinitial)/float(Nz)

        # Use Saha to get x_e, using the equation for x_e for ionized helium
        # and for neutral helium.
        # Everyb_trip ionized above z=8000.  First ionization over by z=5000.
        # Assume He all singly ionized down to z=3500, then use He Saha until
        # He is 99% singly ionized, and *then* switch to joint H/He recombination.
	    z = zend
        if (zend > 8000.)
            x_H0 = 1.
            x_He0 = 1.
            x0 = 1 + 2*𝕣.fHe
            y[1] = x_H0
            y[2] = x_He0
            y[3] = 𝕣.Tnow*(1+z)
        elseif (z > 5000.)
            x_H0 = 1.
            x_He0 = 1.
            rhs = exp(1.5 * log(𝕣.CR*𝕣.Tnow/(1+z)) - 𝕣.CB1_He2/(𝕣.Tnow*(1+z)) ) / 𝕣.Nnow
            rhs = rhs*1.  # ratio of g's is 1 for He++ <-> He+
            x0 = 0.5 * (sqrt( (rhs-1-𝕣.fHe)^2 + 4*(1+2𝕣.fHe)*rhs) - (rhs-1-𝕣.fHe) )
            y[1] = x_H0
            y[2] = x_He0
            y[3] = 𝕣.Tnow*(1+z)
	    elseif (z > 3500.)
            x_H0 = 1.
            x_He0 = 1.
            x0 = x_H0 + 𝕣.fHe*x_He0
            y[1] = x_H0
            y[2] = x_He0
            y[3] = 𝕣.Tnow*(1+z)
        elseif (y[2] > 0.99)
            x_H0 = 1.
            rhs = exp(1.5 * log(𝕣.CR*𝕣.Tnow/(1+z)) - 𝕣.CB1_He1/(𝕣.Tnow*(1+z)) ) / 𝕣.Nnow
            rhs = rhs*4.  # ratio of g's is 4 for He+ <-> He0
            x_He0 = 0.5 * ( sqrt( (rhs-1)^2 + 4*(1+𝕣.fHe)*rhs ) - (rhs-1))
            x0 = x_He0
            x_He0 = (x0 - 1) / 𝕣.fHe
            y[1] = x_H0
            y[2] = x_He0
            y[3] = 𝕣.Tnow*(1+z)
	    elseif (y[1] > 0.99)
            rhs = exp(1.5 * log(𝕣.CR*𝕣.Tnow/(1+z)) - 𝕣.CB1/(𝕣.Tnow*(1+z))) / 𝕣.Nnow
            x_H0 = 0.5 * (sqrt(rhs^2+4*rhs) - rhs)

            prob = ODEProblem{true}(ion_recfast!, y, (zstart, zend), 𝕣,
                save_everystep=false, save_start=false)
            sol = solve(prob, alg, reltol=𝕣.tol)
            y .= sol(zend)

            y[1] = x_H0
            x0 = y[1] + 𝕣.fHe*y[2]
	    else
            prob = ODEProblem{true}(ion_recfast!, y, (zstart, zend), 𝕣,
                save_everystep=false, save_start=false)
            sol = solve(prob, alg, reltol=𝕣.tol)
            y .= sol(zend)
            x0 = y[1] + 𝕣.fHe*y[2]
        end

        Trad = 𝕣.Tnow * (1+zend)
        Tmat = y[3]
        x_H = y[1]
        x_He = y[2]
        x = x0

        out_xe[i] = x
        out_Tmat[i] = Tmat

    end

    return out_xe, out_Tmat
end

RECFASTredshifts(Nz, zinitial, zfinal) =
    range(zinitial, stop=zfinal, length=Nz+1)[2:end]

function IonizationHistory(𝕣::RECFAST{T}, par::ACP, bg::AB) where
                           {T, ACP<:AbstractCosmoParams, AB<:AbstractBackground}
    x_grid = bg.x_grid

    # GRAFT RECFAST ONTO BOLT. TODO: MEGA-REFACTOR ==============
    Nz = 1000
    Xe_RECFAST, Tmat_RECFAST = recfast_xe(𝕣; Nz=Nz, zinitial=10000., zfinal=0.)
    z_RECFAST = RECFASTredshifts(Nz, 10000., 0.)
    RECFAST_Xₑ_z = spline(reverse(Xe_RECFAST), reverse(z_RECFAST))
    RECFAST_Tmat_z = spline(reverse(Tmat_RECFAST), reverse(z_RECFAST))
    xinitial_RECFAST = z2x(first(z_RECFAST))
    Xₑ_function = x -> (x < xinitial_RECFAST) ?
        first(Xe_RECFAST) : RECFAST_Xₑ_z(x2z(x))
    Trad_function = x -> 𝕣.Tnow * (1 + x2z(x))
    Tmat_function = x -> (x < xinitial_RECFAST) ?
        Trad_function(x) : RECFAST_Tmat_z(x2z(x))
    # =====================================================

    τ, τ′ = τ_functions(x_grid, Xₑ_function, par)
    g̃ = g̃_function(τ, τ′)

    Xₑ_ = spline(Xₑ_function.(x_grid), x_grid)
    τ_ = spline(τ.(x_grid), x_grid)
    g̃_ = spline(g̃.(x_grid), x_grid)
    IT = typeof(Xₑ_)

    Tmat_ = spline(Tmat_function.(x_grid), x_grid)
    Trad_ = spline(Trad_function.(x_grid), x_grid)

    # TO FIX, WHY DOES THIS CONSTRUCTOR REQUIRE {I, IT}???
    return IonizationHistory{T, IT}(
        Xₑ_,
        τ_,
        spline_∂ₓ(τ_, x_grid),
        spline_∂ₓ²(τ_, x_grid),
        g̃_,
        spline_∂ₓ(g̃_, x_grid),
        spline_∂ₓ²(g̃_, x_grid),
        Tmat_,
        Trad_
    )
end


text:
module BaseBenchmarks

using BenchmarkTools
using Compat

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1.0
BenchmarkTools.DEFAULT_PARAMETERS.samples = 10000
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 0.15
BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance = 0.01

const PARAMS_PATH = joinpath(dirname(@__FILE__), "..", "etc", "params.json")
const SUITE = BenchmarkGroup()
const MODULES = Dict("array" => :ArrayBenchmarks,
                     "broadcast" => :BroadcastBenchmarks,
                     "collection" => :CollectionBenchmarks,
                     "dates" => :DatesBenchmarks,
                     "find" => :FindBenchmarks,
                     "io" => :IOBenchmarks,
                     "linalg" => :LinAlgBenchmarks,
                     "micro" => :MicroBenchmarks,
                     "misc" => :MiscellaneousBenchmarks,
                     "union" => :UnionBenchmarks,
                     "parallel" => :ParallelBenchmarks,
                     "problem" => :ProblemBenchmarks,
                     "random" => :RandomBenchmarks,
                     "scalar" => :ScalarBenchmarks,
                     "shootout" => :ShootoutBenchmarks,
                     "simd" => :SIMDBenchmarks,
                     "sort" => :SortBenchmarks,
                     "sparse" => :SparseBenchmarks,
                     "string" => :StringBenchmarks,
                     "tuple" => :TupleBenchmarks)

load!(id::AbstractString; kwargs...) = load!(SUITE, id; kwargs...)

function load!(group::BenchmarkGroup, id::AbstractString; tune::Bool = true)
    modsym = MODULES[id]
    modpath = joinpath(dirname(@__FILE__), id, "$(modsym).jl")
    Core.eval(BaseBenchmarks, :(include($modpath)))
    modsuite = Core.eval(BaseBenchmarks, modsym).SUITE
    group[id] = modsuite
    if tune
        results = BenchmarkTools.load(PARAMS_PATH)[1]
        haskey(results, id) && loadparams!(modsuite, results[id], :evals)
    end
    return group
end

loadall!(; kwargs...) = loadall!(SUITE; kwargs...)

function loadall!(group::BenchmarkGroup; verbose::Bool = true, tune::Bool = true)
    for id in keys(MODULES)
        if verbose
            print("loading group $(repr(id))... ")
            time = @elapsed load!(group, id, tune=false)
            println("done (took $time seconds)")
        else
            load!(group, id, tune=false)
        end
    end
    if tune
        results = BenchmarkTools.load(PARAMS_PATH)[1]
        for (id, suite) in group
            haskey(results, id) && loadparams!(suite, results[id], :evals)
        end
    end
    return group
end

end # module


text:
using Test
using ClimaCore.DataLayouts
using StaticArrays
using ClimaCore.DataLayouts: check_basetype, get_struct, set_struct!

@testset "check_basetype" begin
    @test_throws Exception check_basetype(Real, Float64)
    @test_throws Exception check_basetype(Float64, Real)

    @test isnothing(check_basetype(Float64, Float64))
    @test isnothing(check_basetype(Float64, Complex{Float64}))

    @test_throws Exception check_basetype(Float32, Float64)
    @test_throws Exception check_basetype(Float64, Complex{Float32})

    @test isnothing(check_basetype(Float64, Tuple{}))
    @test isnothing(check_basetype(Tuple{}, Tuple{}))
    @test_throws Exception check_basetype(Tuple{}, Float64)

    @test isnothing(check_basetype(Int, Tuple{Int, Complex{Int}}))
    @test isnothing(check_basetype(Float64, typeof(SA[1.0 2.0; 3.0 4.0])))

    S = typeof((a = ((1.0, 2.0f0), (3.0, 4.0f0)), b = (5.0, 6.0f0)))
    @test isnothing(check_basetype(Tuple{Float64, Float32}, S))

    S = typeof(((), (1.0 + 2.0im, NamedTuple()), 3.0 + 4.0im, ()))
    @test isnothing(check_basetype(Float64, S))
    @test isnothing(check_basetype(Complex{Float64}, S))
end

@testset "get_struct / set_struct!" begin
    array = [1.0, 2.0, 3.0]
    S = Tuple{Complex{Float64}, Float64}
    @test get_struct(array, S) == (1.0 + 2.0im, 3.0)
    set_struct!(array, (4.0 + 2.0im, 6.0))
    @test array == [4.0, 2.0, 6.0]
    @test get_struct(array, S) == (4.0 + 2.0im, 6.0)
end

@testset "IJFH" begin
    Nij = 2 # number of nodal points
    Nh = 2 # number of elements
    S = Tuple{Complex{Float64}, Float64}
    array = rand(Nij, Nij, 3, Nh)
    data = IJFH{S, 2}(array)
    @test getfield(data.:1, :array) == @view(array[:, :, 1:2, :])
    data_slab = slab(data, 1)
    @test data_slab[2, 1] ==
          (Complex(array[2, 1, 1, 1], array[2, 1, 2, 1]), array[2, 1, 3, 1])
    data_slab[2, 1] = (Complex(-1.0, -2.0), -3.0)
    @test array[2, 1, 1, 1] == -1.0
    @test array[2, 1, 2, 1] == -2.0
    @test array[2, 1, 3, 1] == -3.0

    subdata_slab = data_slab.:2
    @test subdata_slab[2, 1] == -3.0
    subdata_slab[2, 1] = -5.0
    @test array[2, 1, 3, 1] == -5.0

    @test sum(data.:1) ≈ Complex(sum(array[:, :, 1, :]), sum(array[:, :, 2, :])) atol =
        10eps()
    @test sum(x -> x[2], data) ≈ sum(array[:, :, 3, :]) atol = 10eps()
end

@testset "IJFH boundscheck" begin
    Nij = 1 # number of nodal points
    Nh = 2 # number of elements
    S = Tuple{Complex{Float64}, Float64}
    array = zeros(Float64, Nij, Nij, 3, 2)
    data = IJFH{S, Nij}(array)

    @test_throws BoundsError slab(data, -1)
    @test_throws BoundsError slab(data, 3)
    @test_throws BoundsError slab(data, 1, -1)
    @test_throws BoundsError slab(data, 1, 3)

    # 2D Slab boundscheck
    sdata = slab(data, 1)
    @test_throws BoundsError sdata[-1, 1]
    @test_throws BoundsError sdata[1, -1]
    @test_throws BoundsError sdata[2, 1]
    @test_throws BoundsError sdata[1, 2]
end

@testset "IJFH type safety" begin
    Nij = 2 # number of nodal points per element
    Nh = 1 # number of elements

    # check that types of the same bitstype throw a conversion error
    SA = (a = 1.0, b = 2.0)
    SB = (c = 1.0, d = 2.0)

    array = zeros(Float64, Nij, Nij, 2, Nh)
    data = IJFH{typeof(SA), Nij}(array)
    data_slab = slab(data, 1)
    data_slab[1, 1] = SA
    @test data_slab[1, 1] isa typeof(SA)
    @test_throws MethodError data_slab[1, 1] = SB
end

@testset "broadcasting between data object + scalars" begin
    FT = Float64
    data1 = ones(FT, 2, 2, 2, 2)
    S = Complex{Float64}
    data1 = IJFH{S, 2}(data1)
    res = data1 .+ 1
    @test res isa IJFH{S}
    @test parent(res) ==
          FT[f == 1 ? 2 : 1 for i in 1:2, j in 1:2, f in 1:2, h in 1:2]

    @test sum(res) == Complex(16.0, 8.0)
    @test sum(Base.Broadcast.broadcasted(+, data1, 1)) == Complex(16.0, 8.0)
end

@testset "broadcasting assignment from scalar" begin
    FT = Float64
    S = Complex{FT}
    data = IJFH{S, 2}(Array{FT}, 3)
    data .= Complex(1.0, 2.0)
    @test parent(data) ==
          FT[f == 1 ? 1 : 2 for i in 1:2, j in 1:2, f in 1:2, h in 1:3]

    data .= 1
    @test parent(data) ==
          FT[f == 1 ? 1 : 0 for i in 1:2, j in 1:2, f in 1:2, h in 1:3]

end

@testset "broadcasting between data objects" begin
    FT = Float64
    data1 = ones(FT, 2, 2, 2, 2)
    data2 = ones(FT, 2, 2, 1, 2)
    S1 = Complex{Float64}
    S2 = Float64
    data1 = IJFH{S1, 2}(data1)
    data2 = IJFH{S2, 2}(data2)
    res = data1 .+ data2
    @test res isa IJFH{S1}
    @test parent(res) ==
          FT[f == 1 ? 2 : 1 for i in 1:2, j in 1:2, f in 1:2, h in 1:2]

    @test sum(res) == Complex(16.0, 8.0)
    @test sum(Base.Broadcast.broadcasted(+, data1, data2)) == Complex(16.0, 8.0)
end

@testset "broadcasting complicated function" begin
    FT = Float64
    S1 = NamedTuple{(:a, :b), Tuple{Complex{Float64}, Float64}}
    data1 = ones(FT, 2, 2, 3, 2)
    S2 = Float64
    data2 = ones(FT, 2, 2, 1, 2)
    data1 = IJFH{S1, 2}(data1)
    data2 = IJFH{S2, 2}(data2)

    f(a1, a2) = a1.a.re * a2 + a1.b
    res = f.(data1, data2)
    @test res isa IJFH{Float64}
    @test parent(res) == FT[2 for i in 1:2, j in 1:2, f in 1:1, h in 1:2]
end


text:
# common utilities for getting and comparing column names

_getnames(x::NamedTuple) = propertynames(x)
_getnames(x::AbstractDataFrame) = _names(x)
_getnames(x::DataFrameRow) = _names(x)
_getnames(x::GroupKey) = parent(x).cols

# this function is needed as == does not allow for comparison between tuples and vectors
function _equal_names(r1, r2)
    n1 = _getnames(r1)
    n2 = _getnames(r2)
    length(n1) == length(n2) || return false
    for (a, b) in zip(n1, n2)
        a == b || return false
    end
    return true
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_stopbutton

"""
    daq_stopbutton(;kwargs...)
    daq_stopbutton(children::Any;kwargs...)
    daq_stopbutton(children_maker::Function;kwargs...)

A StopButton component.
A Stop button component
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The children of the button.
- `id` (String; optional): The ID used to identify this compnent in Dash callbacks
- `buttonText` (String; optional): Text displayed in the button
- `className` (String; optional): Class to apply to the root component element.
- `disabled` (Bool; optional): If true, button cannot be pressesd.
- `label` (optional): Description to be displayed alongside the button.
To control styling, pass an object with label and
style properties.. label has the following type: String | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `label` (String; optional)
- `labelPosition` (a value equal to: 'top', 'bottom'; optional): Where the  label is positioned.
- `n_clicks` (Real; optional): Number of times the button was clicked
- `size` (Real; optional): The size (width) of the stop button in pixels
- `style` (Dict; optional): Style to apply to the root component element.
- `theme` (Dict; optional): Theme configuration to be set by a ThemeProvider
"""
function daq_stopbutton(; kwargs...)
        available_props = Symbol[:children, :id, :buttonText, :className, :disabled, :label, :labelPosition, :n_clicks, :size, :style, :theme]
        wild_props = Symbol[]
        return Component("daq_stopbutton", "StopButton", "dash_daq", available_props, wild_props; kwargs...)
end

daq_stopbutton(children::Any; kwargs...) = daq_stopbutton(;kwargs..., children = children)
daq_stopbutton(children_maker::Function; kwargs...) = daq_stopbutton(children_maker(); kwargs...)


text:
function principal_generator_eichler(I::AlgAssAbsOrdIdl)
  O = left_order(I)
  @assert is_maximal(O)
  A = algebra(O)
  @assert is_eichler(A)

  d = discriminant(O)

  y = integral_coprime_representative(O, I, d)
  J = I*y

  N = normred(J, O)
  @assert denominator(N) == 1 # J should be integral
  N = numerator(N)
  @assert is_principal(N) "Ideal is not principal"

  primes = collect(keys(factor(N).fac))
  valuations = [ valuation(N, p) for p in primes ]
  w = _norm_equation_valuations_only(O, primes, valuations)
  w = evaluate(w)

  Ow = O*O(w)
  ww = _eichler_find_transforming_unit(J, Ow)
  return w*ww*inv(y)
end

# Many functions are in AlgAssRelOrd/Eichler.jl

# Finds at least n units in the order F.maximal_orders[order_num]
function _find_some_units(F::FieldOracle{S, T, U, M}, order_num::Int, n::Int) where { S <: AbsAlgAss{fmpq}, T, U, M }
  O = F.maximal_orders[order_num]
  units = Vector{elem_type(O)}()
  while length(units) < n
    KtoA = get_new_field(F, order_num)
    K = domain(KtoA)
    OK = maximal_order(K)
    UK, mUK = unit_group(OK)
    for j = 1:ngens(UK)
      u = mUK(UK[j])
      if u == -one(OK)
        continue
      end
      push!(units, O(KtoA(elem_in_nf(u, copy = false))))
    end
  end
  return units
end

# Finds x\in A^\times such that I == Jx.
# Requires nr(I) == nr(J) and that nr(I) is coprime to d where d is the
# discriminant of O.
function _eichler_find_transforming_unit(I::AlgAssAbsOrdIdl, J::AlgAssAbsOrdIdl)
  O = left_order(I)
  @assert is_maximal(O)
  # We assume that left_order(J) == O, but testing this would be really expensive

  n = normred(I, O)
  @assert n == normred(J, O)
  @assert denominator(n) == 1
  n = numerator(n)

  fac_n = factor(n)
  if isempty(fac_n)
    return one(algebra(order(O)))
  end

  primes = Vector{fmpz}()
  for (p, e) in fac_n
    for i = 1:e
      push!(primes, p)
    end
  end
  t = _eichler_find_transforming_unit_recurse(I, J, primes)
  return t
end


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# This file contains the ExitCriterion struct's definition and its dispatches.

"""
Defines the exit criterion of each algorithm run.
* maxiter: Max #iterations allowed
* maxtime: Max execution time allowed
* targetaccuracy: Target accuracy to halt algorithm
* loggingfreq: #datapass between logging
"""
struct ExitCriterion
    # A data structure to store exit conditions.

    maxiter::Int              # Max #iterations allowed
    maxtime::Float64          # Max execution time allowed
    targetaccuracy::Float64   # Target accuracy to halt algorithm
    loggingfreq::Int          # #datapass between logging
end

"""
Check if the given exit criterion has been satisfied. Returns true if satisfied else returns false.
"""
function checkexitcondition(
    exitcriterion::ExitCriterion,
    currentiter::Integer,
    elapsedtime,
    measure,
)::Bool
    # A function determine if it's time to halt execution

    if currentiter >= exitcriterion.maxiter
        return true
    elseif elapsedtime >= exitcriterion.maxtime
        return true
    elseif measure <= exitcriterion.targetaccuracy
        return true
    end

    false
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
import Base: @deprecate

# """
#     coef(path::RegularizationPath; select=:AICc)
#
# Returns a p by nλ coefficient array where p is the number of
# coefficients (including any intercept) and nλ is the number of path segments,
# or a selected segment's coefficients.
#
# If model was only initialized but not fit, returns a p vector of zeros.
# Consistent with StatsBase.coef, if the model has an intercept it is included.
#
# # Example:
# ```julia
#   m = fit(LassoPath,X,y)
#   coef(m; select=:CVmin)
# ```
#
# # Keywords
# - `select=:all` returns a p by nλ matrix of coefficients
# - `select=:AIC` selects the AIC minimizing segment
# - `select=:AICc` selects the corrected AIC minimizing segment
# - `select=:BIC` selects the BIC minimizing segment
# - `select=:CVmin` selects the segment that minimizing out-of-sample mean squared
#     error from cross-validation with `nCVfolds` random folds.
# - `select=:CV1se` selects the segment whose average OOS deviance is no more than
#     1 standard error away from the one minimizing out-of-sample mean squared
#     error from cross-validation with `nCVfolds` random folds.
# - `nCVfolds=10` number of cross-validation folds
# - `kwargs...` are passed to [`minAICc(::RegularizationPath)`](@ref) or to
#     [`cross_validate_path(::RegularizationPath)`](@ref)
# """
function StatsBase.coef(path::RegularizationPath, select::Symbol; nCVfolds=10)
    if select == :all
        selector = AllSeg()
        msg = "AllSeg()"
    elseif select == :AIC
        selector = MinAIC()
        msg = "MinAIC()"
    elseif select == :AICc
        selector = MinAICc()
        msg = "MinAICc()"
    elseif select == :BIC
        selector = MinBIC()
        msg = "MinBIC()"
    elseif select == :CVmin
        selector = MinCVmse(path, nCVfolds)
        msg = "MinCVmse(path, nCVfolds)"
    elseif select == :CV1se
        selector = MinCV1se(path, nCVfolds)
        msg = "MinCV1se(path, nCVfolds)"
    else
        error("unknown selector $select")
    end

    Base.depwarn("coef(path; select::Symbol=:$select) is deprecated, use coef(path, $msg) instead", :coef)

    coef(path, selector)
end


text:
const Match = Tuple{Rule, Sub, Int64, Bool} #bool isright
const MatchesBuf = OrderedSet{Match}
# const MatchesBuf = Dict{Rule,Set{Sub}}

import ..genrhsfun
import ..options
import ..@log

using .Schedulers

"""
Exactly like [`addexpr!`](@ref), but instantiate pattern variables
from a substitution, resulting from a pattern matcher run.
"""
function addexprinst_rec!(G::EGraph, pat, sub::Sub, side::Symbol)::EClass
    # e = preprocess(pat)
    # println("========== $pat ===========")

    pat isa EClass && return pat

    if haskey(sub, pat)
        (eclass, lit) = sub[pat]
        pat = (lit != nothing ? lit : eclass)
    elseif pat isa Symbol
        error("unbound pattern variable $pat")
    elseif side == :right && pat isa QuoteNode && pat.value isa Symbol
        pat = pat.value
    end

    # println("pat $pat")
    pat isa EClass && return pat

    if istree(pat)
        args = getargs(pat)
        n = length(args)
        class_ids = Vector{Int64}(undef, length(args))
        for i ∈ 1:n
            # println("child $child")
            @inbounds child = args[i]
            c_eclass = addexprinst_rec!(G, child, sub, side)
            @inbounds class_ids[i] = c_eclass.id
        end
        node = ENode(pat, class_ids)
        return add!(G, node)
    end

    return add!(G, ENode(pat))
end

addexprinst!(g::EGraph, e, sub::Sub, side::Symbol)::EClass =
    addexprinst_rec!(g, preprocess(e), sub, side)

function cached_ids(egraph::EGraph, side)::Vector{Int64}
    # outermost symbol in rule side
    # side = remove_assertions(side)
    # side = unquote_sym(side)
    if istree(side)
        sym = gethead(side)
        get(egraph.symcache, sym, [])
    else
        collect(keys(egraph.M))
    end
end

function eqsat_search!(egraph::EGraph, theory::Vector{Rule},
        scheduler::AbstractScheduler)::MatchesBuf
    matches=MatchesBuf()
    for rule ∈ theory
        # don't apply banned rules
        if shouldskip(scheduler, rule)
            # println("skipping banned rule $rule")
            continue
        end
        if rule.mode ∉ [:symbolic, :dynamic, :equational, :unequal]
            error("unsupported mode in rule ", rule)
        end

        ids = cached_ids(egraph, rule.left)

        for id ∈ ids
            for sub in ematch(egraph, rule.left, id)
                # display(sub); println()
                !isempty(sub) && push!(matches, (rule, sub, id, false))
            end
        end

        if rule.mode == :equational || rule.mode == :unequal
            ids = cached_ids(egraph, rule.right)
            for id ∈ ids
                for sub in ematch(egraph, rule.right, id)
                    # display(sub); println()
                    !isempty(sub) && push!(matches, (rule, sub, id, true))
                end
            end
        end
    end
    return matches
end

function eqsat_apply!(egraph::EGraph, matches::MatchesBuf,
        scheduler::AbstractScheduler, report::Report,
        mod::Module, params::SaturationParams)::Report
    i = 0
    for match ∈ matches
        (rule, sub, id, isright) = match
        i += 1

        if i % 300 == 0
            if params.sizeout > 0 && length(egraph.U) > params.sizeout
                @log "E-GRAPH SIZEOUT"
                report.reason = :sizeout
                return report
            end

            if params.stopwhen()
                @log "Halting requirement satisfied"
                report.reason = :condition
                return report
            end
        end

        writestep!(scheduler, rule)

        # l = remove_assertions(rule.left)
        # println("THE ONE IT MATCHED ON = ", lc)
        # lc = addexprinst!(egraph, l, sub, :left).id
        # println("THE INSTANTIATED ONE = ", lc)

        if rule.mode == :symbolic || rule.mode == :equational || rule.mode == :unequal # symbolic replacement
            if isright
                rc = id
                # println("THE ONE IT MATCHED ON = ", rc)
                # schifo = addexprinst!(egraph, rule.right, sub, :right).id
                # println("THE INSTANTIATED ONE = ", schifo)
                lc = addexprinst!(egraph, remove_assertions(rule.left), sub, :left).id
            else
                lc = id
                # println("THE ONE IT MATCHED ON = ", lc)
                # schifo = addexprinst!(egraph, remove_assertions(rule.left), sub, :left).id
                # println("THE INSTANTIATED ONE = ", schifo)
                rc = addexprinst!(egraph, rule.right, sub, :right).id
            end
            if rule.mode == :unequal
                delete!(matches, match)
                if find(egraph, lc) == find(egraph, rc)
                    @log "Contradiction!" rule
                    report.reason = :contradiction
                    return report
                end
            else
                merge!(egraph, lc, rc)
            end
        elseif rule.mode == :dynamic # execute the right hand!
            lc = id
            (params, f) = getrhsfun(rule, mod)
            actual_params = map(params) do x
                (eclass, literal) = sub[x]
                literal != nothing ? literal : eclass
            end
            r = f(geteclass(egraph, lc), egraph, actual_params...)
            rc = addexpr!(egraph, r)
            merge!(egraph,lc,rc.id)
        else
            error("unsupported rule mode")
        end

        # println(rule)
        # println(sub)
        # println(l); println(r)
        # display(egraph.M); println()
    end
    return report
end

"""
Core algorithm of the library: the equality saturation step.
"""
function eqsat_step!(egraph::EGraph, theory::Vector{Rule}, mod::Module,
        scheduler::AbstractScheduler, match_hist::MatchesBuf, params::SaturationParams)

    report = Report(egraph)
    instcache = Dict{Rule, Dict{Sub, Int64}}()

    readstep!(scheduler)

    search_stats = @timed eqsat_search!(egraph, theory, scheduler)
    matches = search_stats.value
    report.search_stats = report.search_stats + discard_value(search_stats)

    # mmm = unique(matches)
    # mmm = symdiff(match_hist, matches)

    matches = setdiff!(matches, match_hist)
    # for (rule, subset) ∈ matches
    #     setdiff!(subset, match_hist[rule])
    # end

    # println("============ WRITE PHASE ============")
    # println("\n $(length(matches)) $(length(match_hist))")
    # if length(matches) > matchlimit
    #     matches = matches[1:matchlimit]
    # #     # mmm = Set(collect(mmm)[1:matchlimit])
    # #     #
    # #     # report.reason = :matchlimit
    # #     # @goto quit_rebuild
    # end
    # println(" diff length $(length(matches))")

    apply_stats = @timed eqsat_apply!(egraph, matches,
        scheduler, report, mod, params)
    report = apply_stats.value
    report.apply_stats = report.apply_stats + discard_value(apply_stats)

    # don't add inequalities to history
    union!(match_hist, matches)

    # display(egraph.parents); println()
    # display(egraph.M); println()
    @label quit_rebuild
    report.saturated = isempty(egraph.dirty) && report.reason == nothing
    rebuild_stats = @timed rebuild!(egraph)
    report.rebuild_stats = report.rebuild_stats + discard_value(rebuild_stats)

    # TODO produce proofs with match_hist
    total_time!(report)

    return report, egraph
end

"""
Given an [`EGraph`](@ref) and a collection of rewrite rules,
execute the equality saturation algorithm.
"""
saturate!(egraph::EGraph, theory::Vector{Rule}; mod=@__MODULE__) =
    saturate!(egraph, theory, SaturationParams(); mod=mod)

function saturate!(egraph::EGraph, theory::Vector{Rule}, params::SaturationParams;
    mod=@__MODULE__,)

    curr_iter = 0

    # init the scheduler
    sched = params.scheduler(egraph, theory, params.schedulerparams...)

    match_hist = MatchesBuf()

    tot_report = Report()

    # GC.enable(false)
    while true
        curr_iter+=1

        options[:printiter] && @info("iteration ", curr_iter)

        report, egraph = eqsat_step!(egraph, theory, mod, sched, match_hist, params)

        tot_report = tot_report + report

        # report.reason == :matchlimit && break
        report.reason == :condition && break
        report.reason == :contradiction && break
        cansaturate(sched) && report.saturated && (tot_report.saturated = true; break)
        curr_iter >= params.timeout && (tot_report.reason = :timeout; break)
        params.sizeout > 0 && length(egraph.U) > params.sizeout && (tot_report.reason = :sizeout; break)
        params.stopwhen() && (tot_report.reason = :condition; break)
    end
    # println(match_hist)

    # display(egraph.M); println()
    tot_report.iterations = curr_iter
    @log tot_report
    # GC.enable(true)

    return tot_report
end


text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(TriComplex(0, 0.0))
end

@test begin
    iszerodivisor(TriComplex(0,1,1,0))
end

@test begin
    !iszerodivisor(TriComplex(1,2,3,4))
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(TriComplex(a))
end

@test begin
    a = 1
    !isreal(TriComplex(a, a))
end

@test begin
    z = random(TriComplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(TriComplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, TriComplex(1, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1, H: 2, G: 3, GH: 4, F: 5, FH: 6, FG: 7, FGH: 8]"
    l == r
end

@test begin
    x = random(TriComplex{BigFloat})
    y = random(TriComplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(TriComplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(TriComplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(BiComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    barH(barH(x)) == x
end

@test begin
    x = random(TriComplex{BigInt})
    barG(barG(x)) == x
end

@test begin
    x = random(TriComplex{BigInt})
    barF(barF(x)) == x
end

@test begin
    x = random(TriComplex{BigInt})
    barG(barH(x)) == barH(barG(x))
end

@test begin
    x = random(TriComplex{BigInt})
    barF(barH(x)) == barH(barF(x))
end

@test begin
    x = random(TriComplex{BigInt})
    barG(barF(x)) == barF(barG(x))
end

@test begin
    x = random(TriComplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barH(x * y)
    r = barH(x) * barH(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barH(x + y)
    r = barH(x) + barH(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barH(x - y)
    r = barH(x) - barH(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barG(x * y)
    r = barG(x) * barG(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barG(x + y)
    r = barG(x) + barG(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barG(x - y)
    r = barG(x) - barG(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barF(x * y)
    r = barF(x) * barF(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barF(x + y)
    r = barF(x) + barF(y)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = barF(x - y)
    r = barF(x) - barF(y)
    l == r
end

@test begin
    x = random(TriComplex{Rational{BigInt}})
    y = random(TriComplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(TriComplex{Rational{BigInt}})
    y = random(TriComplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(TriComplex{Rational{BigInt}})
    y = random(TriComplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(TriComplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    norm(x) > big(0)
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = norm(x * y)
    r = norm(x) * norm(y)
    l == r
end

@test_throws ErrorException begin
    inv(TriComplex(0, 0.0))
end

@test_throws ErrorException begin
    random(TriComplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(TriComplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(TriComplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(TriComplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(TriComplex{BigInt})
    z == z + zero(TriComplex{BigInt})
end

@test begin
    z = random(TriComplex{BigInt})
    z == z * one(TriComplex{BigInt})
end

@test begin
    z = random(TriComplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(TriComplex{BigInt})
    z == z * one(z)
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = commutator(x, y)
    r = zero(TriComplex{BigInt})
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = alternatorL(x, y)
    r = zero(TriComplex{BigInt})
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = alternatorR(x, y)
    r = zero(TriComplex{BigInt})
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    l = flexator(x, y)
    r = zero(TriComplex{BigInt})
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = associator(x, y, z)
    r = zero(TriComplex{BigInt})
    l == r
end

@test begin
    x = random(TriComplex{BigInt})
    y = random(TriComplex{BigInt})
    z = random(TriComplex{BigInt})
    l = jacobiator(x, y, z)
    r = 3 * x * y * z
    l == r
end

@test begin
    l = TriComplex(Cmplex(1,2))
    r = TriComplex(1,2)
    l == r
end

@test begin
    l = TriComplex(BiComplex(1,2))
    r = TriComplex(1,2,0)
    l == r
end

@test begin
    l = TriComplex(1, 2, 3.0)
    r = TriComplex(1.0, 2.0, 3.0)
    l == r
end

@test begin
    l = TriComplex(1, 2, 3, 4.0)
    r = TriComplex(1.0, 2.0, 3.0, 4.0)
    l == r
end

@test begin
    l = TriComplex(1, 2, 3, 4, 5.0)
    r = TriComplex(1.0, 2.0, 3.0, 4.0, 5.0)
    l == r
end

@test begin
    l = TriComplex(1, 2, 3, 4, 5, 6.0)
    r = TriComplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
    l == r
end

@test begin
    l = TriComplex(1, 2, 3, 4, 5, 6, 7.0)
    r = TriComplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0)
    l == r
end

@test begin
    l = TriComplex(1, 2, 3, 4, 5, 6, 7, 8.0)
    r = TriComplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
    l == r
end


text:
using Base: @pure, @propagate_inbounds

@inline function fllog2(x::T) where {T <: Integer}
    return (sizeof(T) * 8 - 1) - leading_zeros(x)
end
@inline function cllog2(x::T) where {T <: Integer}
    return fllog2(x - 1) + 1
end

@inline function fld2(x::T) where {T <: Integer}
    return x >>> true
end
@inline function cld2(x::T) where {T <: Integer}
    return (x + true) >>> true
end

@pure nbits(::Type{T}) where {T} = sizeof(T) * 8
@pure log2nbits(::Type{T}) where {T} = flrlog2(nbits(T))

@inline function flpow1m(x::T) where {T <: Integer}
    return (1 << x) - 1
end

@inline undefs(T::Type, dims::Vararg{Any, N}) where {N} = Array{T, N}(undef, dims...)

zero!(arr) = fill!(arr, zero(eltype(arr)))

le_nest(f, x, ns) = _le_nest(f, x, ns...)
_le_nest(f, x, n) = f(n)
function _le_nest(f, x, n, ns...)
    return quote
        if $x <= $n
            $(f(n))
        else
            $(_le_nest(f, x, ns...))
        end
    end
end

eq_nest(f, x, ns) = _eq_nest(f, x, ns...)
_eq_nest(f, x, n) = f(n)
function _eq_nest(f, x, n, ns...)
    return quote
        if $x == $n
            $(f(n))
        else
            $(_eq_nest(f, x, ns...))
        end
    end
end

@memoize arch_id() = (sha2_256(string(cpuinfo()))...,)


text:
#Contains MPI related things that the library is missing
import MPI

#TODO document

# The RSend constand isn't defined on windows
const MPI_RSEND = if Sys.iswindows()
                      (:MPI_RSEND, MPI.libmpi)
                  else
                      MPI.MPI_RSEND
                  end

function MPI_Rsend(buf::MPI.MPIBuffertype{T}, count::Integer,
                dest::Integer, tag::Integer, comm::MPI.Comm) where T
    ccall(MPI_RSEND, Nothing,
        (Ptr{T}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint},
           Ref{Cint}),
        buf, count, MPI.mpitype(T), dest, tag, comm.val, 0)
end

function MPI_Rsend(buf::AbstractArray{T}, dest::Integer, tag::Integer, comm::MPI.Comm) where T
    MPI_Rsend(buf, length(buf), dest, tag, comm)
end

function MPI_Rsend(obj::T, dest::Integer, tag::Integer, comm::MPI.Comm) where T
    MPI_Rsend([obj], dest, tag, comm)
end


text:
struct PolynomialBasis
    polynomial::AbstractVector{T} where T <: AbstractPolynomialLike
    x::AbstractVector{T} where T <: AbstractPolynomialLike
end

"""
# Notes
`with_bias=true` provides polynomials with the degree of from 0 to `d`.
`with_bias=false` provides polynomials with the degree of `d`.
# Variables
n ∈ N: length of array, i.e., x ∈ Rⁿ
d ∈ N: degree
"""
function PolynomialBasis(n::Int, d::Int; with_bias=true)
    @assert n >= 1 && d >= 0
    _n = with_bias ? n+1 : n
    exponents = multiexponents(_n, d)
    @polyvar x[1:n]
    _x = with_bias ? [x..., 1] : x
    polynomial = exponents |> Map(exponent -> prod(_x.^exponent)) |> collect
    PolynomialBasis(polynomial, x)
end

function (a::PolynomialBasis)(x::AbstractArray)
    a.polynomial |> Map(poly -> poly(a.x => x)) |> collect
end


text:
module BiobakeryUtils

export
    import_abundance_tables,
    import_abundance_table,
    clean_abundance_tables,
    taxfilter,
    taxfilter!,
    rm_strat!,
    permanova,
    humann2_regroup,
    humann2_rename,
    humann2_barplot,
    humann2_barplots,
    qvalue!

using DataFrames
using Statistics
using CSV
using RCall
using Microbiome

include("general.jl")
include("metaphlan.jl")
include("humann2.jl")

end


text:
module LIBSVM_

export LinearSVC, SVC
export NuSVC, NuSVR
export EpsilonSVR
export OneClassSVM

import MLJModelInterface
import MLJModelInterface: Table, Continuous, Count, Finite, OrderedFactor,
                          Multiclass

const MMI = MLJModelInterface

import ..LIBSVM
using Statistics

"""
    LinearSVC(; kwargs...)

Linear support vector machine classifier using LIBLINEAR: https://www.csie.ntu.edu.tw/~cjlin/liblinear/

See also SVC, NuSVC
"""
mutable struct LinearSVC <: MMI.Deterministic
    solver::LIBSVM.Linearsolver.LINEARSOLVER
    weights::Union{Dict, Nothing}
    tolerance::Float64
    cost::Float64
    p::Float64
    bias::Float64
end

function LinearSVC(
    ;solver::LIBSVM.Linearsolver.LINEARSOLVER = LIBSVM.Linearsolver.L2R_L2LOSS_SVC_DUAL
    ,weights::Union{Dict, Nothing} = nothing
    ,tolerance::Float64 = Inf
    ,cost::Float64 = 1.0
    ,p::Float64 = 0.1
    ,bias::Float64= -1.0)

    model = LinearSVC(
        solver
        ,weights
        ,tolerance
        ,cost
        ,p
        ,bias
    )

    message = MMI.clean!(model)   #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

"""
    SVC(; kwargs...)

Kernel support vector machine classifier using LIBSVM: https://www.csie.ntu.edu.tw/~cjlin/libsvm/
If `gamma==-1.0` then  `gamma = 1/nfeatures is used in
fitting.
If `gamma==0.0` then a `gamma = 1/(var(X) * nfeatures)` is
used in fitting

See also LinearSVC, NuSVC
"""
mutable struct SVC <: MMI.Deterministic
    kernel::LIBSVM.Kernel.KERNEL
    gamma::Float64
    weights::Union{Dict, Nothing}
    cost::Float64
    cachesize::Float64
    degree::Int32
    coef0::Float64
    tolerance::Float64
    shrinking::Bool
    probability::Bool
end

function SVC(
    ;kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis
    ,gamma::Float64 = 0.0
    ,weights::Union{Dict, Nothing} = nothing
    ,cost::Float64 = 1.0
    ,cachesize::Float64=200.0
    ,degree::Int32 = Int32(3)
    ,coef0::Float64 = 0.0
    ,tolerance::Float64 = .001
    ,shrinking::Bool = true
    ,probability::Bool = false)

    model = SVC(
        kernel
        ,gamma
        ,weights
        ,cost
        ,cachesize
        ,degree
        ,coef0
        ,tolerance
        ,shrinking
        ,probability
    )

    message = MMI.clean!(model)   #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

"""
    NuSVC(; kwargs...)

Kernel support vector machine classifier using LIBSVM: https://www.csie.ntu.edu.tw/~cjlin/libsvm/

If `gamma==-1.0` then  `gamma = 1/nfeatures is used in
fitting.
If `gamma==0.0` then a `gamma = 1/(var(X) * nfeatures)` is
used in fitting

See also LinearSVC, SVC
"""
mutable struct NuSVC <: MMI.Deterministic
    kernel::LIBSVM.Kernel.KERNEL
    gamma::Float64
    weights::Union{Dict, Nothing}
    nu::Float64
    cost::Float64
    cachesize::Float64
    degree::Int32
    coef0::Float64
    tolerance::Float64
    shrinking::Bool
end

function NuSVC(
    ;kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis
    ,gamma::Float64 = 0.0
    ,weights::Union{Dict, Nothing} = nothing
    ,nu::Float64 = 0.5
    ,cost::Float64 = 1.0
    ,cachesize::Float64 = 200.0
    ,degree::Int32 = Int32(3)
    ,coef0::Float64 = 0.
    ,tolerance::Float64 = .001
    ,shrinking::Bool = true)

    model = NuSVC(
        kernel
        ,gamma
        ,weights
        ,nu
        ,cost
        ,cachesize
        ,degree
        ,coef0
        ,tolerance
        ,shrinking
    )

    message = MMI.clean!(model)   #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

mutable struct OneClassSVM <: MMI.Unsupervised
    kernel::LIBSVM.Kernel.KERNEL
    gamma::Float64
    nu::Float64
    cost::Float64
    cachesize::Float64
    degree::Int32
    coef0::Float64
    tolerance::Float64
    shrinking::Bool
end

function OneClassSVM(
    ;kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis
    ,gamma::Float64 = 0.0
    ,nu::Float64 = 0.1
    ,cost::Float64 = 1.0
    ,cachesize::Float64 = 200.0
    ,degree::Int32 = Int32(3)
    ,coef0::Float64 = 0.0
    ,tolerance::Float64 = .001
    ,shrinking::Bool = true)

    model = OneClassSVM(
        kernel
        ,gamma
        ,nu
        ,cost
        ,cachesize
        ,degree
        ,coef0
        ,tolerance
        ,shrinking
    )

    message = MMI.clean!(model)   #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

"""
    NuSVR(; kwargs...)

Kernel support vector machine regressor using LIBSVM: https://www.csie.ntu.edu.tw/~cjlin/libsvm/

If `gamma==-1.0` then  `gamma = 1/nfeatures is used in
fitting.
If `gamma==0.0` then a `gamma = 1/(var(X) * nfeatures)` is
used in fitting

See also EpsilonSVR
"""
mutable struct NuSVR <: MMI.Deterministic
    kernel::LIBSVM.Kernel.KERNEL
    gamma::Float64
    nu::Float64
    cost::Float64
    cachesize::Float64
    degree::Int32
    coef0::Float64
    tolerance::Float64
    shrinking::Bool
end

function NuSVR(
    ;kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis
    ,gamma::Float64 = 0.0
    ,nu::Float64 = 0.5
    ,cost::Float64 = 1.0
    ,cachesize::Float64 = 200.0
    ,degree::Int32 = Int32(3)
    ,coef0::Float64 = 0.
    ,tolerance::Float64 = .001
    ,shrinking::Bool = true)

    model = NuSVR(
        kernel
        ,gamma
        ,nu
        ,cost
        ,cachesize
        ,degree
        ,coef0
        ,tolerance
        ,shrinking
    )

    message = MMI.clean!(model)   #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

"""
    EpsilonSVR(; kwargs...)

Kernel support vector machine regressor using LIBSVM: https://www.csie.ntu.edu.tw/~cjlin/libsvm/

If `gamma==-1.0` then  `gamma = 1/nfeatures is used in
fitting.
If `gamma==0.0` then a `gamma = 1/(var(X) * nfeatures)` is
used in fitting

See also NuSVR
"""
mutable struct EpsilonSVR <: MMI.Deterministic
    kernel::LIBSVM.Kernel.KERNEL
    gamma::Float64
    epsilon::Float64
    cost::Float64
    cachesize::Float64
    degree::Int32
    coef0::Float64
    tolerance::Float64
    shrinking::Bool
end

function EpsilonSVR(
    ;kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis
    ,gamma::Float64 = 0.0
    ,epsilon::Float64 = 0.1
    ,cost::Float64 = 1.0
    ,cachesize::Float64 = 200.0
    ,degree::Int32 = Int32(3)
    ,coef0::Float64 = 0.
    ,tolerance::Float64 = .001
    ,shrinking::Bool = true)

    model = EpsilonSVR(
        kernel
        ,gamma
        ,epsilon
        ,cost
        ,cachesize
        ,degree
        ,coef0
        ,tolerance
        ,shrinking
    )

    message = MMI.clean!(model)   #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

const SVM = Union{LinearSVC, SVC, NuSVC, NuSVR, EpsilonSVR, OneClassSVM} # all SVM models defined here

"""
    map_model_type(model::SVM)

Helper function to map the model to the correct LIBSVM model type needed for function dispatch.
"""
function map_model_type(model::SVM)
    if isa(model, LinearSVC)
        return LIBSVM.LinearSVC
    elseif isa(model, SVC)
        return LIBSVM.SVC
    elseif isa(model, NuSVC)
        return LIBSVM.NuSVC
    elseif isa(model, NuSVR)
        return LIBSVM.NuSVR
    elseif isa(model, EpsilonSVR)
        return LIBSVM.EpsilonSVR
    elseif isa(model, OneClassSVM)
        return LIBSVM.OneClassSVM
    else
        error("Got unsupported model type: $(typeof(model))")
    end
end

"""
    get_svm_parameters(model::Union{SVC, NuSVC, NuSVR, EpsilonSVR, OneClassSVM})

Helper function to get the parameters from the SVM model struct.
"""
function get_svm_parameters(model::Union{SVC, NuSVC, NuSVR, EpsilonSVR, OneClassSVM})
    #Build arguments for calling svmtrain
    params = Tuple{Symbol, Any}[]
    push!(params, (:svmtype, map_model_type(model))) # get LIBSVM model type
    for fn in fieldnames(typeof(model))
        push!(params, (fn, getfield(model, fn)))
    end

    return params
end

function MMI.fit(model::LinearSVC, verbosity::Int, X, y)

    Xmatrix = MMI.matrix(X)' # notice the transpose
    y_plain = MMI.int(y)
    decode  = MMI.decoder(y[1]) # for predict method

    cache = nothing

    result = LIBSVM.LIBLINEAR.linear_train(y_plain, Xmatrix,
        weights = model.weights, solver_type = Int32(model.solver),
        C = model.cost, p = model.p, bias = model.bias,
        eps = model.tolerance, verbose = ifelse(verbosity > 1, true, false)
    )

    fitresult = (result, decode)
    report = nothing

    return fitresult, cache, report
end

function MMI.fit(model::Union{SVC, NuSVC}, verbosity::Int, X, y)

    Xmatrix = MMI.matrix(X)' # notice the transpose
    y_plain = MMI.int(y)
    decode  = MMI.decoder(y[1]) # for predict method

    cache = nothing

    model = deepcopy(model)
    model.gamma == -1.0 && (model.gamma = 1.0/size(Xmatrix, 1))
    model.gamma == 0.0 && (model.gamma = 1.0/(var(Xmatrix) * size(Xmatrix, 1)) )
    result = LIBSVM.svmtrain(Xmatrix, y_plain;
        get_svm_parameters(model)...,
        verbose = ifelse(verbosity > 1, true, false)
    )

    fitresult = (result, decode)
    report = (gamma=model.gamma,)

    return fitresult, cache, report
end

function MMI.fit(model::Union{NuSVR, EpsilonSVR}, verbosity::Int, X, y)

    Xmatrix = MMI.matrix(X)' # notice the transpose

    cache = nothing

    model = deepcopy(model)
    model.gamma == -1.0 && (model.gamma = 1.0/size(Xmatrix, 1))
    model.gamma == 0.0 && (model.gamma = 1.0/(var(Xmatrix) * size(Xmatrix, 1)) )
    fitresult = LIBSVM.svmtrain(Xmatrix, y;
        get_svm_parameters(model)...,
        verbose = ifelse(verbosity > 1, true, false)
    )

    report = (gamma=model.gamma,)

    return fitresult, cache, report
end

function MMI.fit(model::OneClassSVM, verbosity::Int, X)

    Xmatrix = MMI.matrix(X)' # notice the transpose

    cache = nothing

    model = deepcopy(model)
    model.gamma == -1.0 && (model.gamma = 1.0/size(Xmatrix, 1))
    model.gamma == 0.0 && (model.gamma = 1.0/(var(Xmatrix) * size(Xmatrix, 1)) )
    fitresult = LIBSVM.svmtrain(Xmatrix;
        get_svm_parameters(model)...,
        verbose = ifelse(verbosity > 1, true, false)
    )

    report = (gamma=model.gamma,)

    return fitresult, cache, report
end

function MMI.predict(model::LinearSVC, fitresult, Xnew)
    result, decode = fitresult
    (p,d) = LIBSVM.LIBLINEAR.linear_predict(result, MMI.matrix(Xnew)')
    return decode(p)
end

function MMI.predict(model::Union{SVC, NuSVC}, fitresult, Xnew)
    result, decode = fitresult
    (p,d) = LIBSVM.svmpredict(result, MMI.matrix(Xnew)')
    return decode(p)
end

function MMI.predict(model::Union{NuSVR, EpsilonSVR}, fitresult, Xnew)
    (p,d) = LIBSVM.svmpredict(fitresult, MMI.matrix(Xnew)')
    return p
end

function MMI.transform(model::OneClassSVM, fitresult, Xnew)
    (p,d) = LIBSVM.svmpredict(fitresult, MMI.matrix(Xnew)')
    return MMI.categorical(p)
end

# metadata
MMI.load_path(::Type{<:LinearSVC}) = "MLJModels.LIBSVM_.LinearSVC"
MMI.load_path(::Type{<:SVC}) = "MLJModels.LIBSVM_.SVC"
MMI.load_path(::Type{<:NuSVC}) = "MLJModels.LIBSVM_.NuSVC"
MMI.load_path(::Type{<:NuSVR}) = "MLJModels.LIBSVM_.NuSVR"
MMI.load_path(::Type{<:EpsilonSVR}) = "MLJModels.LIBSVM_.EpsilonSVR"
MMI.load_path(::Type{<:OneClassSVM}) = "MLJModels.LIBSVM_.OneClassSVM"

MMI.package_name(::Type{<:SVM}) = "LIBSVM"
MMI.package_uuid(::Type{<:SVM}) = "b1bec4e5-fd48-53fe-b0cb-9723c09d164b"
MMI.is_pure_julia(::Type{<:SVM}) = false
MMI.package_url(::Type{<:SVM}) = "https://github.com/mpastell/LIBSVM.jl"
MMI.input_scitype(::Type{<:SVM}) = Table(Continuous)
MMI.target_scitype(::Type{<:Union{LinearSVC, SVC, NuSVC}}) = AbstractVector{<:Finite}
MMI.target_scitype(::Type{<:Union{NuSVR, EpsilonSVR}}) = AbstractVector{Continuous}
MMI.output_scitype(::Type{<:OneClassSVM}) = AbstractVector{<:Finite{2}} # Bool (true means inlier)

end # module


text:
# Autogenerated wrapper script for Ipopt_jll for i686-w64-mingw32-libgfortran3-cxx11
export amplexe, libipopt, libipoptamplinterface

using ASL_jll
using MUMPS_seq_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Ipopt")
JLLWrappers.@declare_executable_product(amplexe)
JLLWrappers.@declare_library_product(libipopt, "libipopt-3.dll")
JLLWrappers.@declare_library_product(libipoptamplinterface, "libipoptamplinterface-3.dll")
function __init__()
    JLLWrappers.@generate_init_header(ASL_jll, MUMPS_seq_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        amplexe,
        "bin\\ipopt.exe",
    )

    JLLWrappers.@init_library_product(
        libipopt,
        "bin\\libipopt-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libipoptamplinterface,
        "bin\\libipoptamplinterface-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
if !(@isdefined peaks)
    include("peaks.jl")
end

####################
## main
####################

# Parameters
const N=31
xmat, ymat, P = peaks(N)
xvec = vec(xmat[1,:])
yvec = vec(ymat[:,1])

# find min
minval, ind = findmin(P)
cind = CartesianIndices(P)[ind]
row_min, col_min = cind[1], cind[2]

# find max
maxval, ind = findmax(P)
cind = CartesianIndices(P)[ind]
row_max, col_max = cind[1], cind[2]

####################

####################
## plot
####################

# directory where figures are printed
figdir="./fig"
if !isdir(figdir); mkdir(figdir); end

using Plots
pyplot()
#gr()

# see http://docs.juliaplots.org/latest/colors/#colorschemes
clibrary(:misc)

# Two-dimensional contour
contour(xvec, yvec, P, c=(:rainbow), fill=false, tickfont=12,
        xlims=(-3.,3.), ylims=(-3.,3.), clims=(-6.,8.),
        xlabel="X", ylabel="Y", axis_ratio=:equal)
# max and min
scatter!([xvec[col_min]], [yvec[row_min]], ms=10., color=:yellow, lab="min")
scatter!([xvec[col_max]], [yvec[row_max]], ms=10., color=:magenta, lab="max",
         legendfont=12,legend=:bottomleft, size=(800,600))
# save
savefig(joinpath(figdir,"contour_2d.png"))

# Three dimensional surface plot
# plot(..., linetype=:surface) , or surface(...)
# plot(xvec,yvec,P,c=(:rainbow),linetype=:surface,fillalpha=0.9, tickfont=12,
surface(xvec,yvec,P,c=(:rainbow),fillalpha=0.9, tickfont=12,
        xlims=(-3.,3.), ylims=(-3.,3.), zlims=(-8.,10.), clims=(-6.,8.),
        xlabel="X", ylabel="Y", zlabel="Z", size=(800,600), colorbar=:best)
# max and min
scatter!([xvec[col_min]], [yvec[row_min]], [minval], ms=8., color=:yellow, lab="min")
scatter!([xvec[col_max]], [yvec[row_max]], [maxval], ms=8., color=:magenta, lab="max",
         legendfont=12, legend=:bottomleft)
# save
savefig(joinpath(figdir,"surface_3d.png"))
#######################


text:
using SimpleWorkflows
using Test

@testset "SimpleWorkflow.jl" begin
    include("workflow.jl")
end


text:
__precompile__()
module RCall

using Requires
using Dates
using Libdl
using Random
using REPL
if VERSION ≤ v"1.1.1"
   using Missings
end
using CategoricalArrays
using DataFrames
using StatsModels

import DataStructures: OrderedDict

import Base: eltype, convert, isascii,
    names, length, size, getindex, setindex!,
    show, showerror, write
import Base.Iterators: iterate, IteratorSize, IteratorEltype, Pairs, pairs

export RObject,
   Sxp, NilSxp, StrSxp, CharSxp, LglSxp, IntSxp, RealSxp, CplxSxp,
   ListSxp, VecSxp, EnvSxp, LangSxp, ClosSxp, S4Sxp,
   getattrib, setattrib!, getnames, setnames!, getclass, setclass!, attributes,
   globalEnv,
   isnull, isna, anyna,
   robject, rcopy, rparse, rprint, reval, rcall, rlang,
   rimport, @rimport, @rlibrary, @rput, @rget, @var_str, @R_str

const depfile = joinpath(dirname(@__FILE__),"..","deps","deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("RCall not properly installed. Please run Pkg.build(\"RCall\")")
end

include("setup.jl")
include("types.jl")
include("Const.jl")
include("methods.jl")
include("convert/base.jl")
include("convert/missing.jl")
include("convert/categorical.jl")
include("convert/datetime.jl")
include("convert/dataframe.jl")
include("convert/formula.jl")

include("convert/default.jl")
include("eventloop.jl")
include("eval.jl")
include("language.jl")
include("io.jl")
include("callback.jl")
include("namespaces.jl")
include("render.jl")
include("macros.jl")
include("operators.jl")
include("RPrompt.jl")
include("ijulia.jl")
include("deprecated.jl")

end # module


text:
using CGP
using Logging
using PyCall
using ArgParse
using Printf
using Random

@pyimport gym
@pyimport pybullet
@pyimport pybullet_envs
CGP.Config.init("cfg/base.yaml")
CGP.Config.init("cfg/classic.yaml")

function play_env(c::Chromosome, env, render::Bool=false)
    env[:seed](0)
    if render
        env[:render](mode="human")
    end
    ob = env[:reset]()
    total_reward = 0.0
    done = false
    reward = 0.0

    while ~done
        action = process(c, ob./5.0)
        try
            ob, reward, done, _ = env[:step](action)
            if render
                env[:render](mode="human")
            end
        catch
            done = true
        end
        total_reward += reward
    end

    total_reward
end

function get_args()
    s = ArgParseSettings()

    @add_arg_table(
        s,
        "--seed", arg_type = Int, default = 0,
        "--log", arg_type = String, default = "gym.log",
        "--id", arg_type = String, default = "MountainCarContinuous-v0",
        "--ea", arg_type = String, default = "oneplus",
        "--chromosome", arg_type = String, default = "CGPChromo",
        "--cfg", arg_type = String,
    )

    parse_args(CGP.Config.add_arg_settings!(s))
end

if ~isinteractive()
    args = get_args()

    if args["cfg"] != nothing
        CGP.Config.init(args["cfg"])
    end

    CGP.Config.init(Dict([k=>args[k] for k in setdiff(
        keys(args), ["seed", "log", "id", "ea", "chromosome", "cfg", "graph"])]...))

    Random.seed!(args["seed"])
    global_logger(SimpleLogger(open(args["log"], "a+")))
    ea = eval(Meta.parse(args["ea"]))
    ctype = eval(Meta.parse(args["chromosome"]))
    env = gym.make(args["id"])
    nin = length(env[:observation_space][:low])
    nout = length(env[:action_space][:low])

    fit = x->play_env(x, env, false)
    maxfit, best = ea(nin, nout, fit; seed=args["seed"], ctype=ctype, id=args["id"])

    @info(@sprintf("E%0.6f", -maxfit))

    best_ind = ctype(best, nin, nout)
    play_env(best_ind, env, true)
end


text:
include("functions_collapse.jl")

function pairsnp_square(nuc_bit_array_A, nuc_bit_V)

    diff_list_A = get_list_of_differences(nuc_bit_array_A, nuc_bit_V)

    snp_count_array = zeros(Int, size(nuc_bit_array_A, 2), size(nuc_bit_array_A, 2))

    # For parallelising:
    n_Chunks = Threads.nthreads()
    if n_Chunks > size(nuc_bit_array_A, 2) / 2
        n_Chunks::Int = floor(size(nuc_bit_array_A, 2) / 2)
    end
    ranges = get_ranges(size(nuc_bit_array_A, 2), n_Chunks)

    Threads.@threads for range in ranges
        for a in range[1]:range[2]
            for b in (a + 1):size(nuc_bit_array_A, 2)

                snp_counter = 0

                for r in diff_list_A[a]
                    @inbounds x = nuc_bit_array_A[r,a]
                    @inbounds y = nuc_bit_array_A[r,b]
                    different = (x & y) < 16

                    if different
                        snp_counter+=1
                    end
                end

                for r in diff_list_A[b]
                    if r in diff_list_A[a]
                        continue
                    end

                    @inbounds x = nuc_bit_array_A[r,a]
                    @inbounds y = nuc_bit_array_A[r,b]
                    different = (x & y) < 16

                    if different
                        snp_counter+=1
                    end
                end

                snp_count_array[a,b] = snp_counter
                snp_count_array[b,a] = snp_counter
            end
        end
    end

    return snp_count_array
end

function pairsnp_rectangle(nuc_bit_array_A, nuc_bit_array_B, nuc_bit_V)

    diff_list_A = get_list_of_differences(nuc_bit_array_A, nuc_bit_V)
    diff_list_B = get_list_of_differences(nuc_bit_array_B, nuc_bit_V)

    snp_count_array = zeros(Int, size(nuc_bit_array_B, 2), size(nuc_bit_array_A, 2))

    Threads.@threads for a in 1:size(nuc_bit_array_A, 2)
        for b in 1:size(nuc_bit_array_B, 2)

            snp_counter = 0

            for r in diff_list_A[a]
                @inbounds x = nuc_bit_array_A[r,a]
                @inbounds y = nuc_bit_array_B[r,b]
                different = (x & y) < 16

                if different
                    snp_counter+=1
                end
            end

            for r in diff_list_B[b]
                if r in diff_list_A[a]
                    continue
                end

                @inbounds x = nuc_bit_array_A[r,a]
                @inbounds y = nuc_bit_array_B[r,b]
                different = (x & y) < 16

                if different
                    snp_counter+=1
                end
            end

            snp_count_array[b,a] = snp_counter
        end
    end

    return snp_count_array
end

function write_pairsnp_square(outfile, int_array, seq_name_vector)

    open(outfile, "w") do io
        @inbounds for i in 1:size(int_array, 2)
            write(io, seq_name_vector[i])
            write(io, ",")

            write(io, join(int_array[:,i], ",") * '\n')
        end
    end
end

function write_pairsnp_rectangle(outfile, int_array, seq_name_vector_A, seq_name_vector_B)

    open(outfile, "w") do io
        write(io, ",")
        write(io, join(seq_name_vector_B, ",") * "\n")

        @inbounds for i in 1:size(int_array, 2)
            write(io, seq_name_vector_A[i])
            write(io, ",")
            write(io, join(int_array[:,i], ",") * '\n')
        end
    end
end

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
#=
In this file the core computations for creating a recurrence matrix
are defined (via multiple dispatch).

The low level interface is contained in the function
`recurrence_matrix`, and this is where any specialization should happen.
=#

################################################################################
# Distance matrix
################################################################################
"""
    distancematrix(x [, y = x], metric = "euclidean")

Create a matrix with the distances between each pair of points of the
time series `x` and `y` using `metric`.

The time series `x` and `y` can be `Dataset`s or vectors or matrices with data points
in rows.
The data point dimensions (or number of columns) must be the same for `x` and `y`.
The returned value is a `n×m` matrix, with `n` being the length (or number of rows)
of `x`, and `m` the length of `y`.

The metric can be identified by a string, or any of the `Metric`s defined in
the [`Distances` package](https://github.com/JuliaStats/Distances.jl).
The list of strings available to define the metric are:

* `"max"` or `"inf"` for the maximum or L∞ norm
  (`Chebyshev()` in the `Distances` package).
* `"euclidean"` for the L2 or Euclidean norm, used by default
  (`Euclidean()` in `Distances`).
* `"manhattan"`, `"cityblock"`, `"taxicab"` or `"min"` for the Manhattan or L1 norm
  (`Cityblock()` in `Distances`).
"""
distancematrix(x, metric::Union{Metric,String}=DEFAULT_METRIC) = distancematrix(x, x, metric)

# For 1-dimensional arrays (vectors), the distance does not depend on the metric
distancematrix(x::Vector, y::Vector, metric=DEFAULT_METRIC) = abs.(x .- y')

# If the metric is supplied as a string, get the corresponding Metric from Distances
distancematrix(x, y, metric::String) = distancematrix(x, y, getmetric(metric))

const MAXDIM = 9
function distancematrix(x::Tx, y::Ty, metric::Metric=DEFAULT_METRIC) where
         {Tx<:Union{AbstractMatrix, Dataset}} where {Ty<:Union{AbstractMatrix, Dataset}}
    sx, sy = size(x), size(y)
    if sx[2] != sy[2]
        error("the dimensions of `x` and `y` data points must be the equal")
    end
    if sx[2] ≥ MAXDIM && typeof(metric) == Euclidean # Blas optimization
        return _distancematrix(Matrix(x), Matrix(y), metric)
    elseif Tx <: Matrix && Ty <: Matrix && metric == Chebyshev()
        return _distancematrix(x, y, metric)
    else
        return _distancematrix(Dataset(x), Dataset(y), metric)
    end
end

# Core function for Matrices (wrapper of `pairwise` from the Distances package)
_distancematrix(x::AbstractMatrix, y::AbstractMatrix, metric::Metric) =
pairwise(metric, x', y')
# Core function for Datasets
function _distancematrix(x::Dataset{S,Tx}, y::Dataset{S,Ty},
    metric::Metric) where {S, Tx, Ty}

    x = x.data
    y = y.data
    d = zeros(promote_type(Tx,Ty), length(x), length(y))
    for j in 1:length(y)
        for i in 1:length(x)
            @inbounds d[i,j] = evaluate(metric, x[i], y[j])
        end
    end
    return d
end

################################################################################
# AbstractRecurrenceMatrix type definitions and documentation strings
################################################################################
abstract type AbstractRecurrenceMatrix end
const ARM = AbstractRecurrenceMatrix
struct RecurrenceMatrix <: AbstractRecurrenceMatrix
    data::SparseMatrixCSC{Bool,Int}
end
struct CrossRecurrenceMatrix <: AbstractRecurrenceMatrix
    data::SparseMatrixCSC{Bool,Int}
end
struct JointRecurrenceMatrix <: AbstractRecurrenceMatrix
    data::SparseMatrixCSC{Bool,Int}
end

function Base.summary(R::AbstractRecurrenceMatrix)
    N = nnz(R.data)
    return "$(nameof(typeof(R))) of size $(size(R.data)) with $N entries"
end
Base.show(io::IO, R::AbstractRecurrenceMatrix) = println(io, summary(R))

# Propagate used functions:
begin
    extentions = [
        (:Base, (:getindex, :size, :length, :view, :iterate)),
        (:LinearAlgebra, (:diag, :triu, :tril, :issymmetric)),
        (:SparseArrays, (:nnz, :rowvals, :nzrange, :nonzeros))
    ]
    for (M, fs) in extentions
        for f in fs
            @eval $M.$(f)(x::ARM, args...) = $(f)(x.data, args...)
        end
    end
end

for operator in [:(==), :(!=)]
    @eval Base.$operator(x::ARM, y::ARM) = $operator(x.data, y.data)
end

LinearAlgebra.issymmetric(::RecurrenceMatrix) = true
# column values in sparse matrix (parallel to rowvals)
function colvals(x::SparseMatrixCSC)
    cv = zeros(Int,nnz(x))
    @inbounds for c=1:size(x,2)
        cv[nzrange(x,c)] .= c
    end
    cv
end
colvals(x::ARM) = colvals(x.data)

"""
    RecurrenceMatrix(x, ε; kwargs...)

Create a recurrence matrix from trajectory `x`.
Objects of type `<:AbstractRecurrenceMatrix` are displayed as a [`recurrenceplot`](@ref).

## Description

The recurrence matrix is a numeric representation of a "recurrence plot" [1, 2],
in the form of a sparse square matrix of Boolean values.

`x` must be a `Vector` or a `Dataset` or a `Matrix` with data points in rows
(possibly representing and embedded phase space; see [`embed`](@ref)).
If `d(x[i], x[j]) ≤ ε` (with `d` the distance function),
then the cell `(i, j)` of the matrix will have a `true`
value. The criteria to evaluate distances between data points are defined
by the following keyword arguments:

* `scale=1` : a function of the distance matrix (see [`distancematrix`](@ref)),
  or a fixed number, used to scale the value of `ε`. Typical choices are
  `maximum` or `mean`, such that the threshold `ε` is defined as a ratio of the
  maximum or the mean distance between data points, respectively (using
  `mean` or `maximum` calls specialized versions that are faster than the naive
  approach).  Use `1` to keep the distances unscaled (default).
* `fixedrate::Bool=false` : a flag that indicates if `ε` should be
  taken as a target fixed recurrence rate (see [`recurrencerate`](@ref)).
  If `fixedrate` is set to `true`, `ε` must be a value between 0 and 1,
  and `scale` is ignored.
* `metric="euclidean"` : metric of the distances, either `Metric` or a string,
   as in [`distancematrix`](@ref).

See also: [`CrossRecurrenceMatrix`](@ref), [`JointRecurrenceMatrix`](@ref) and
use [`recurrenceplot`](@ref) to turn the result of these functions into a plottable format.

## References
[1] : N. Marwan *et al.*, "Recurrence plots for the analysis of complex systems",
*Phys. Reports 438*(5-6), 237-329 (2007).

[2] : N. Marwan & C.L. Webber, "Mathematical and computational foundations of
recurrence quantifications", in: Webber, C.L. & N. Marwan (eds.), *Recurrence
Quantification Analysis. Theory and Best Practices*, Springer, pp. 3-43 (2015).
"""
function RecurrenceMatrix(x, ε; metric = DEFAULT_METRIC, kwargs...)
    m = getmetric(metric)
    s = resolve_scale(x, m, ε; kwargs...)
    m = recurrence_matrix(x, m, s; kwargs...)
    return RecurrenceMatrix(m)
end

"""
    CrossRecurrenceMatrix(x, y, ε; kwargs...)

Create a cross recurrence matrix from trajectories `x` and `y`.

The cross recurrence matrix is a bivariate extension of the recurrence matrix.
For the time series `x`, `y`, of length `n` and `m`, respectively, it is a
sparse `n×m` matrix of Boolean values, such that if `d(x[i], y[j]) ≤ ε`,
then the cell `(i, j)` of the matrix will have a `true` value.

See [`RecurrenceMatrix`](@ref) for details, references and keywords.
See also: [`JointRecurrenceMatrix`](@ref).
"""
function CrossRecurrenceMatrix(x, y, ε; metric = DEFAULT_METRIC, kwargs...)
    m = getmetric(metric)
    s = resolve_scale(x, y, m, ε; kwargs...)
    m = recurrence_matrix(x, y, m, s)
    return CrossRecurrenceMatrix(m)
end

"""
    JointRecurrenceMatrix(x, y, ε; kwargs...)

Create a joint recurrence matrix from `x` and `y`.

The joint recurrence matrix considers the recurrences of the trajectories
of `x` and `y` separately, and looks for points where both recur
simultaneously. It is calculated by the element-wise multiplication
of the recurrence matrices of `x` and `y`. If `x` and `y` are of different
length, the recurrences are only calculated until the length of the shortest one.

See [`RecurrenceMatrix`](@ref) for details, references and keywords.
See also: [`CrossRecurrenceMatrix`](@ref).
"""
function JointRecurrenceMatrix(x, y, ε; kwargs...)
    n = min(size(x,1), size(y,1))
    if n == size(x,1) && n == size(y,1)
        rm1 = RecurrenceMatrix(x, ε; kwargs...)
        rm2 = RecurrenceMatrix(y, ε; kwargs...)
    else
        rm1 = RecurrenceMatrix(x[1:n,:], ε; kwargs...)
        rm2 = RecurrenceMatrix(y[1:n,:], ε; kwargs...)
    end
    return JointRecurrenceMatrix(rm1.data .* rm2.data)
end

################################################################################
# Scaling / fixed rate
################################################################################
# here args... is (x, y, metric, ε) or just (x, metric, ε)
function resolve_scale(args...; scale=1, fixedrate=false)
    ε = args[end]
    # Check fixed recurrence rate - ε must be within (0, 1)
    if fixedrate
        sfun = (m) -> quantile(vec(m), ε)
        return resolve_scale(Base.front(args)..., 1.0; scale=sfun, fixedrate=false)
    else
        scale_value = _computescale(scale, Base.front(args)...)
        return ε*scale_value
    end
end

# If `scale` is a function, compute the numeric value of the scale based on the
# distance matrix; otherwise return the value of `scale` itself
_computescale(scale::Real, args...) = scale
_computescale(scale::Function, x, metric::Metric) =
_computescale(scale, x, x, metric)

# generic method that uses `distancematrix`
function _computescale(scale::Function, x, y, metric)
    if x===y
        distances = zeros(Int(length(x)*(length(x)-1)/2))
        c = 0
        @inbounds for i in 1:length(x)-1, j=(i+1):length(x)
            distances[c+=1] = evaluate(metric, x[i], y[j])
        end
    else
        distances = distancematrix(x, y, metric)
    end
    return scale(distances)
end
# specific methods to avoid `distancematrix`
function _computescale(scale::typeof(maximum), x, y, metric::Metric)
    maxvalue = zero(eltype(x))
    if x===y
        @inbounds for i in 1:length(x)-1, j=(i+1):length(x)
            newvalue = evaluate(metric, x[i], y[j])
            (newvalue > maxvalue) && (maxvalue = newvalue)
        end
    else
        @inbounds for xi in x, yj in y
            newvalue = evaluate(metric, xi, yj)
            (newvalue > maxvalue) && (maxvalue = newvalue)
        end
    end
    return maxvalue
end
function _computescale(scale::typeof(mean), x, y, metric::Metric)
    meanvalue = 0.0
    if x===y
        @inbounds for i in 1:length(x)-1, j=(i+1):length(x)
            meanvalue += evaluate(metric, x[i], y[j])
        end
        denominator = length(x)*(length(x)-1)/2
    else
        @inbounds for xi in x, yj in y
            meanvalue += evaluate(metric, xi, yj)
        end
        denominator = Float64(length(x)*length(y))
    end
    return meanvalue/denominator
end

################################################################################
# recurrence_matrix - Low level interface
################################################################################
# TODO: increase the efficiency here by not computing everything:
recurrence_matrix(x, metric::Metric, ε::Real) =
recurrence_matrix(x, x, metric, ε)

# Convert Matrices to Datasets (better performance in all cases)
function recurrence_matrix(x::AbstractMatrix, y::AbstractMatrix,
                           metric::Metric, ε)
    return recurrence_matrix(Dataset(x), Dataset(y), metric, ε)
end

# Core function
function recurrence_matrix(xx::Dataset, yy::Dataset, metric::Metric, ε)
    x = xx.data
    y = yy.data
    rowvals = Vector{Int}()
    colvals = Vector{Int}()
    for j in 1:length(y)
        nzcol = 0
        for i in 1:length(x)
            @inbounds if evaluate(metric, x[i], y[j]) ≤ ε
                push!(rowvals, i)
                nzcol += 1
            end
        end
        append!(colvals, fill(j, (nzcol,)))
    end
    nzvals = fill(true, (length(rowvals),))
    return sparse(rowvals, colvals, nzvals, length(x), length(y))
end

# Vector version can be more specialized (and metric is irrelevant)
function recurrence_matrix(x::AbstractVector, y::AbstractVector, metric, ε)
    rowvals = Vector{Int}()
    colvals = Vector{Int}()
    for j in 1:length(y)
        nzcol = 0
        for i in 1:length(x)
            if @inbounds abs(x[i] - y[j]) ≤ ε
                push!(rowvals, i)
                nzcol += 1
            end
        end
        append!(colvals, fill(j, (nzcol,)))
    end
    nzvals = fill(true, (length(rowvals),))
    return sparse(rowvals, colvals, nzvals, length(x), length(y))
end


text:
@testset "logisticsoftmax" begin
    N, d = 100, 1
    K = 3
    k = 2.0 * SqExponentialKernel() ∘ ScaleTransform(10.0)
    X, f1 = generate_f(N, d, k)
    X, f2 = generate_f(N, d, k; X=X)
    X, f3 = generate_f(N, d, k; X=X)
    f = [f1, f2, f3]
    y = getindex.(findmax(hcat(f1, f2, f3); dims=2)[2], 2)[:]
    # scatter(X,[f1,f2,f3])
    # scatter!(X,y)

    floattypes = [Float64]
    tests_likelihood(
        LogisticSoftMaxLikelihood(K),
        MultiClassLikelihood,
        Dict(
            "VGP" => Dict("AVI" => true, "QVI" => false, "MCVI" => true),
            "SVGP" => Dict("AVI" => true, "QVI" => false, "MCVI" => true),
            "OSVGP" => Dict("AVI" => true, "QVI" => false, "MCVI" => false),
            "MCGP" => Dict("Gibbs" => true, "HMC" => false),
        ),
        floattypes,
        "MultiClass",
        K,
        X,
        f,
        y,
        k,
    )
end


text:
using Documenter
using MarketData

makedocs(
  format = Documenter.HTML(),
  sitename = "MarketData.jl",
  pages = [
    "index.md",
    "getting_started.md",
    "test_data.md",
    "company_financial_series.md",
    "downloads.md",
  ],
  modules = [MarketData],
)

deploydocs(
  repo         = "github.com/JuliaQuant/MarketData.jl.git",
  devbranch    = "master",
  target       = "build",
  deps         = nothing,  # we use the `format = :html`, without `mkdocs`
  make         = nothing,  # we use the `format = :html`, without `mkdocs`
  push_preview = true,
)