import logging
import math
import random
import sys
from collections.abc import MutableSequence, Sequence

from .logging import get_logger

logger = get_logger(__name__, logging.ERROR)

# increase recursion limit for recursive quicksort
sys.setrecursionlimit(3000)

_MIN = 1e-9
_MAX = 1e9


def compute_moiety_conservation_laws(
    stoichiometric_list: Sequence[float],
    num_species: int,
    num_reactions: int,
    max_num_monte_carlo: int = 20,
    rng_seed: None | bool | int = False,
    species_names: Sequence[str] | None = None,
) -> tuple[list[list[int]], list[list[float]]]:
    """Compute moiety conservation laws.

    According to the algorithm proposed by De Martino et al. (2014)
    https://doi.org/10.1371/journal.pone.0100750

    :param stoichiometric_list:
        the stoichiometric matrix as a list (species x reactions,
        column-major ordering)
    :param num_species:
        total number of species in the reaction network
    :param num_reactions:
        total number of reactions in the reaction network
    :param max_num_monte_carlo:
        maximum number of MonteCarlo steps before changing to relaxation
    :param rng_seed:
        Seed for the random number generator. If `False`, the RNG will not be
        re-initialized. Other values will be passed to :func:`random.seed`.
    :param species_names:
        Species names. Optional and only used for logging.
    :returns:
        Integer MCLs as list of lists of indices of involved species and
        list of lists of corresponding coefficients.
    """
    # compute semi-positive conservation laws
    (
        kernel_dim,
        engaged_species,
        int_kernel_dim,
        conserved_moieties,
        cls_species_idxs,
        cls_coefficients,
    ) = _kernel(stoichiometric_list, num_species, num_reactions)
    # if the number of integer MCLs equals total MCLS no MC relaxation
    done = int_kernel_dim == kernel_dim

    if not done:
        # construct interaction matrix
        J, J2, fields = _fill(
            stoichiometric_list, engaged_species, num_species
        )

        # seed random number generator
        if rng_seed is not False:
            random.seed(rng_seed)

        timer = 0
        # maximum number of montecarlo search before starting relaxation
        while not done:
            yes, int_kernel_dim, conserved_moieties = _monte_carlo(
                engaged_species,
                J,
                J2,
                fields,
                conserved_moieties,
                int_kernel_dim,
                cls_species_idxs,
                cls_coefficients,
                num_species,
                max_iter=max_num_monte_carlo,
            )
            # if the number of integer MCLs equals total MCLS then MC done
            done = int_kernel_dim == kernel_dim
            timer = 0 if yes else timer + 1

            if timer == max_num_monte_carlo:
                done = _relax(
                    stoichiometric_list,
                    conserved_moieties,
                    num_reactions,
                    num_species,
                )
                timer = 0
    _reduce(int_kernel_dim, cls_species_idxs, cls_coefficients, num_species)
    _output(
        int_kernel_dim,
        kernel_dim,
        engaged_species,
        cls_species_idxs,
        cls_coefficients,
        species_names,
        verbose=True,
    )

    return cls_species_idxs[:int_kernel_dim], cls_coefficients[:int_kernel_dim]


def _output(
    int_kernel_dim: int,
    kernel_dim: int,
    int_matched: list[int],
    species_indices: list[list[int]],
    species_coefficients: list[list[float]],
    species_names: Sequence[str] | None = None,
    verbose: bool = False,
    log_level: int = logging.DEBUG,
):
    """Log infos on identified conservation laws"""

    def log(*args, **kwargs):
        logger.log(log_level, *args, **kwargs)

    log(
        f"There are {int_kernel_dim} linearly independent conserved "
        f"moieties, engaging {len(int_matched)} state variables."
    )
    if int_kernel_dim == kernel_dim:
        log("They generate all the conservation laws")
    else:
        log(
            f"They don't generate all the conservation laws, "
            f"{kernel_dim - int_kernel_dim} of them are not reducible to "
            "moieties"
        )
    # print all conserved quantities
    if verbose:
        for i, (coefficients, engaged_species_idxs) in enumerate(
            zip(species_coefficients, species_indices, strict=True)
        ):
            if not engaged_species_idxs:
                continue
            log(
                f"Moiety number {i + 1} engages {len(engaged_species_idxs)} "
                "species:"
            )
            for species_idx, coefficient in zip(
                engaged_species_idxs, coefficients, strict=True
            ):
                name = (
                    species_names[species_idx]
                    if species_names
                    else species_idx
                )
                log(f"\t{name}\t{coefficient}")


def _qsort(
    k: int, km: int, order: MutableSequence[int], pivots: Sequence[int]
) -> None:
    """Quicksort

    Recursive implementation of the quicksort algorithm

    :param k:
        number of elements to sort
    :param km:
        current center element
    :param order:
        ordering of the elements
    :param pivots:
        corresponding pivot elements from scaled partial pivoting strategy
    """
    # TODO: Rewrite into an iterative algorithm with pivoting strategy

    if k - km < 1:
        # nothing to sort
        return

    pivot = km + int((k - km) / 2)
    l = 0
    p = k - km - 1
    new_order = [0] * (k - km)
    for i in range(km, k):
        if i != pivot:
            if pivots[order[i]] < pivots[order[pivot]]:
                new_order[l] = order[i]
                l += 1
            else:
                new_order[p] = order[i]
                p -= 1
    new_order[p] = order[pivot]
    order[km:k] = new_order

    # calculate center, then recursive calls on left and right intervals
    centre = p + km
    _qsort(k, centre + 1, order, pivots)
    _qsort(centre, km, order, pivots)


def _kernel(
    stoichiometric_list: Sequence[float], num_species: int, num_reactions: int
) -> tuple[int, list[int], int, list[int], list[list[int]], list[list[float]]]:
    """
    Kernel (left nullspace of :math:`S`) calculation by Gaussian elimination

    To compute the left nullspace of the stoichiometric matrix :math:`S`,
    a Gaussian elimination method with partial scaled pivoting is used to deal
    effectively with a possibly ill-conditioned stoichiometric matrix
    :math:`S`.

    Note that this is the Python reimplementation of the algorithm proposed by
    `De Martino et al. (2014) <https://doi.org/10.1371/journal.pone.0100750>`_
    and thus a direct adaption of the original implementation in C++.

    :param stoichiometric_list:
        the stoichiometric matrix as a list (species x reactions,
        col-major ordering)
    :param num_species:
        total number of species in the reaction network
    :param num_reactions:
        total number of reactions in the reaction network
    :returns:
        kernel dimension, MCLs, integer kernel dimension, integer MCLs and
        indices to species and reactions in the preceding order as a tuple
    """
    matrix: list[list[int]] = [[] for _ in range(num_species)]
    matrix2: list[list[float]] = [[] for _ in range(num_species)]
    i_reaction = 0
    i_species = 0
    for val in stoichiometric_list:
        if val != 0:
            matrix[i_species].append(i_reaction)
            matrix2[i_species].append(val)
        i_species += 1
        if i_species == num_species:
            i_species = 0
            i_reaction += 1
    for i in range(num_species):
        matrix[i].append(num_reactions + i)
        matrix2[i].append(1)

    order: list[int] = list(range(num_species))
    pivots = [
        matrix[i][0] if len(matrix[i]) else _MAX for i in range(num_species)
    ]

    done = False
    while not done:
        _qsort(num_species, 0, order, pivots)
        for j in range(num_species - 1):
            if pivots[order[j + 1]] == pivots[order[j]] != _MAX:
                min1 = _MAX
                if len(matrix[order[j]]) > 1:
                    for i in range(len(matrix[order[j]])):
                        min1 = min(
                            min1,
                            abs(matrix2[order[j]][0] / matrix2[order[j]][i]),
                        )

                min2 = _MAX
                if len(matrix[order[j + 1]]) > 1:
                    for i in range(len(matrix[order[j + 1]])):
                        min2 = min(
                            min2,
                            abs(
                                matrix2[order[j + 1]][0]
                                / matrix2[order[j + 1]][i]
                            ),
                        )

                if min2 > min1:
                    # swap
                    k2 = order[j + 1]
                    order[j + 1] = order[j]
                    order[j] = k2
        done = True

        for j in range(num_species - 1):
            if pivots[order[j + 1]] == pivots[order[j]] != _MAX:
                k1 = order[j + 1]
                k2 = order[j]
                column: list[float] = [0] * (num_species + num_reactions)
                g = matrix2[k2][0] / matrix2[k1][0]
                for i in range(1, len(matrix[k1])):
                    column[matrix[k1][i]] = matrix2[k1][i] * g

                for i in range(1, len(matrix[k2])):
                    column[matrix[k2][i]] -= matrix2[k2][i]

                matrix[k1] = []
                matrix2[k1] = []
                for col_idx, col_val in enumerate(column):
                    if abs(col_val) > _MIN:
                        matrix[k1].append(col_idx)
                        matrix2[k1].append(col_val)

                done = False
                if len(matrix[order[j + 1]]):
                    pivots[order[j + 1]] = matrix[order[j + 1]][0]
                else:
                    pivots[order[j + 1]] = _MAX

    RSolutions = [[] for _ in range(num_species)]
    RSolutions2 = [[] for _ in range(num_species)]
    kernel_dim = 0

    for i in range(num_species):
        done = all(
            matrix[i][j] >= num_reactions for j in range(len(matrix[i]))
        )
        if done and len(matrix[i]):
            for j in range(len(matrix[i])):
                RSolutions[kernel_dim].append(matrix[i][j] - num_reactions)
                RSolutions2[kernel_dim].append(matrix2[i][j])
            kernel_dim += 1
    del matrix, matrix2

    matched = []
    int_matched = []
    cls_species_idxs = [[] for _ in range(num_species)]
    cls_coefficients = [[] for _ in range(num_species)]

    i2 = 0
    for i in range(kernel_dim):
        ok2 = True
        for j in range(len(RSolutions[i])):
            if RSolutions2[i][j] * RSolutions2[i][0] < 0:
                ok2 = False
            if not matched or all(
                cur_matched != RSolutions[i][j] for cur_matched in matched
            ):
                matched.append(RSolutions[i][j])
        if ok2 and len(RSolutions[i]):
            min_value = _MAX
            for j in range(len(RSolutions[i])):
                cls_species_idxs[i2].append(RSolutions[i][j])
                cls_coefficients[i2].append(abs(RSolutions2[i][j]))
                min_value = min(min_value, abs(RSolutions2[i][j]))
                if not int_matched or all(
                    cur_int_matched != cls_species_idxs[i2][j]
                    for cur_int_matched in int_matched
                ):
                    int_matched.append(cls_species_idxs[i2][j])
            for j in range(len(cls_species_idxs[i2])):
                cls_coefficients[i2][j] /= min_value
            i2 += 1
    int_kernel_dim = i2

    assert int_kernel_dim <= kernel_dim
    assert len(cls_species_idxs) == len(cls_coefficients), (
        "Inconsistent number of conserved quantities in coefficients and "
        "species"
    )
    return (
        kernel_dim,
        matched,
        int_kernel_dim,
        int_matched,
        cls_species_idxs,
        cls_coefficients,
    )


def _fill(
    stoichiometric_list: Sequence[float],
    matched: Sequence[int],
    num_species: int,
) -> tuple[list[list[int]], list[list[int]], list[int]]:
    """Construct interaction matrix

    Construct the interaction matrix out of the given stoichiometric matrix
    :math:`S`.

    :param stoichiometric_list:
        the stoichiometric matrix given as a flat list
    :param matched:
        found and independent moiety conservation laws (MCL)
    :param num_species:
        number of rows in :math:`S`
    :returns:
        interactions of metabolites and reactions, and matrix of interaction
    """
    dim = len(matched)

    # for each entry in the stoichiometric matrix save interaction
    i_reaction = 0
    i_species = 0
    matrix = [[] for _ in range(dim)]
    matrix2 = [[] for _ in range(dim)]
    for val in stoichiometric_list:
        if val != 0:
            take = dim
            for matched_idx, matched_val in enumerate(matched):
                if i_species == matched_val:
                    take = matched_idx
            if take < dim:
                matrix[take].append(i_reaction)
                matrix2[take].append(val)
        i_species += 1
        if i_species == num_species:
            i_species = 0
            i_reaction += 1

    J = [[] for _ in range(num_species)]
    J2 = [[] for _ in range(num_species)]
    fields = [0] * num_species
    for i in range(dim):
        for j in range(i, dim):
            interactions = 0
            for po in range(len(matrix[i])):
                for pu in range(len(matrix[j])):
                    if matrix[i][po] == matrix[j][pu]:
                        interactions += matrix2[i][po] * matrix2[j][pu]
            if j == i:
                fields[i] = interactions
            elif abs(interactions) > _MIN:
                J[i].append(j)
                J2[i].append(interactions)
                J[j].append(i)
                J2[j].append(interactions)
    return J, J2, fields


def _is_linearly_dependent(
    vector: Sequence[float],
    int_kernel_dim: int,
    cls_species_idxs: Sequence[Sequence[int]],
    cls_coefficients: Sequence[Sequence[float]],
    matched: Sequence[int],
    num_species: int,
) -> bool:
    """Check for linear dependence between MCLs

    Check if the solutions found with Monte Carlo are linearly independent
    with respect to the previous found solution for all MCLs involved

    :param vector:
        found basis
    :param int_kernel_dim:
        number of integer conservative laws
    :param cls_species_idxs:
        NSolutions contains the species involved in the MCL
    :param cls_coefficients:
        NSolutions2 contains the corresponding coefficients in the MCL
    :param matched:
        actual found MCLs
    :param num_species:
        number of rows in :math:`S`
    :returns:
        boolean indicating linear dependence (true) or not (false)
    """
    K = int_kernel_dim + 1
    matrix: list[list[int]] = [[] for _ in range(K)]
    matrix2: list[list[float]] = [[] for _ in range(K)]
    # Populate matrices with species ids and coefficients for CLs
    for i in range(K - 1):
        for j in range(len(cls_species_idxs[i])):
            matrix[i].append(cls_species_idxs[i][j])
            matrix2[i].append(cls_coefficients[i][j])

    order2 = list(range(len(matched)))
    pivots2 = matched[:]
    _qsort(len(matched), 0, order2, pivots2)

    # ensure positivity
    for i in range(len(matched)):
        if vector[order2[i]] > _MIN:
            matrix[K - 1].append(matched[order2[i]])
            matrix2[K - 1].append(vector[order2[i]])

    order = list(range(K))
    pivots = [matrix[i][0] if len(matrix[i]) else _MAX for i in range(K)]

    # check for linear independence of the solution
    ok = False
    while not ok:
        _qsort(K, 0, order, pivots)
        for j in range(K - 1):
            if pivots[order[j + 1]] == pivots[order[j]] != _MAX:
                min1 = _MAX
                if len(matrix[order[j]]) > 1:
                    for i in range(len(matrix[order[j]])):
                        min1 = min(
                            min1,
                            abs(matrix2[order[j]][0] / matrix2[order[j]][i]),
                        )
                min2 = _MAX
                if len(matrix[order[j + 1]]) > 1:
                    for i in range(len(matrix[order[j + 1]])):
                        min2 = min(
                            min2,
                            abs(
                                matrix2[order[j + 1]][0]
                                / matrix2[order[j + 1]][i]
                            ),
                        )
                if min2 > min1:
                    # swap
                    k2 = order[j + 1]
                    order[j + 1] = order[j]
                    order[j] = k2
        ok = True
        for j in range(K - 1):
            if pivots[order[j + 1]] == pivots[order[j]] != _MAX:
                k1 = order[j + 1]
                k2 = order[j]
                column: list[float] = [0] * num_species
                g = matrix2[k2][0] / matrix2[k1][0]
                for i in range(1, len(matrix[k1])):
                    column[matrix[k1][i]] = matrix2[k1][i] * g
                for i in range(1, len(matrix[k2])):
                    column[matrix[k2][i]] -= matrix2[k2][i]

                matrix[k1] = []
                matrix2[k1] = []
                for i in range(num_species):
                    if abs(column[i]) > _MIN:
                        matrix[k1].append(i)
                        matrix2[k1].append(column[i])
                ok = False
                pivots[k1] = matrix[k1][0] if len(matrix[k1]) else _MAX
    K1 = sum(len(matrix[i]) > 0 for i in range(K))
    return K == K1


def _monte_carlo(
    matched: Sequence[int],
    J: Sequence[Sequence[int]],
    J2: Sequence[Sequence[float]],
    fields: Sequence[float],
    int_matched: MutableSequence[int],
    int_kernel_dim: int,
    cls_species_idxs: MutableSequence[MutableSequence[int]],
    cls_coefficients: MutableSequence[MutableSequence[float]],
    num_species: int,
    initial_temperature: float = 1,
    cool_rate: float = 1e-3,
    max_iter: int = 10,
) -> tuple[bool, int, Sequence[int]]:
    """MonteCarlo simulated annealing for finding integer MCLs

    Finding integer solutions for the MCLs by Monte Carlo, see step (b) in
    the De Martino (2014) paper and Eqs. 11-13 in the publication

    :param matched:
        number of found MCLs
    :param J:
        index of metabolites involved in a MCL
    :param J2:
        coefficients of metabolites involved in a MCL
    :param fields:
        actual number of metabolites involved in a MCL
    :param int_matched:
        actual matched MCLs
    :param int_kernel_dim:
        number of MCLs found in :math:`S`
    :param cls_species_idxs:
        Modified in-place.
    :param cls_coefficients:
        Modified in-place.
    :param initial_temperature:
        initial temperature
    :param cool_rate:
        cooling rate of simulated annealing
    :param max_iter:
        maximum number of MonteCarlo steps before changing to relaxation
    :returns:
        status of MC iteration, number of integer MCLs, number of MCLs,
        metabolites and reaction indices, MCLs and integer MCLs as a tuple

        status indicates if the currently found moiety by the Monte Carlo
        process is linearly dependent (False) or linearly independent (True)
        in case of linear dependence, the current Monte Carlo cycle can be
        considered otherwise the algorithm retries Monte Carlo up to max_iter
    """
    dim = len(matched)
    num = [
        int(2 * random.uniform(0, 1)) if len(J[i]) else 0 for i in range(dim)
    ]
    numtot = sum(num)

    def compute_h():
        H = 0
        for i in range(dim):
            H += fields[i] * num[i] ** 2
            for j in range(len(J[i])):
                H += J2[i][j] * num[i] * num[J[i][j]]
        return H

    H = compute_h()

    count = 0
    howmany = 0
    T1 = initial_temperature
    e = math.exp(-1 / T1)
    while True:
        en = int(random.uniform(0, 1) * dim)
        while not len(J[en]):
            en = int(random.uniform(0, 1) * dim)

        p = -1 if num[en] > 0 and random.uniform(0, 1) < 0.5 else 1
        delta = fields[en] * num[en]
        for i in range(len(J[en])):
            delta += J2[en][i] * num[J[en][i]]
        delta = 2 * p * delta + fields[en]

        if delta < 0 or random.uniform(0, 1) < math.pow(e, delta):
            num[en] += p
            numtot += p
            H += delta

        count += 1
        if count % dim == 0:
            T1 -= cool_rate
            if T1 <= 0:
                T1 = cool_rate
            e = math.exp(-1 / T1)

        if count == dim // cool_rate:
            count = 0
            T1 = initial_temperature
            e = math.exp(-1 / T1)
            en = int(random.uniform(0, 1) * dim)
            while not len(J[en]):
                en = int(random.uniform(0, 1) * dim)
            num = [0] * dim
            num[en] = 1
            numtot = 1

            H = compute_h()
            howmany += 1

        if (H < _MIN and numtot > 0) or (howmany == 10 * max_iter):
            break

    if howmany >= 10 * max_iter:
        return False, int_kernel_dim, int_matched

    # founds MCLS? need to check for linear independence
    if len(int_matched) and not _is_linearly_dependent(
        num,
        int_kernel_dim,
        cls_species_idxs,
        cls_coefficients,
        matched,
        num_species,
    ):
        logger.debug("Found a moiety but it is linearly dependent... next.")
        return False, int_kernel_dim, int_matched

    # reduce by MC procedure
    order2 = list(range(len(matched)))
    pivots2 = matched[:]
    _qsort(len(matched), 0, order2, pivots2)
    for i in range(len(matched)):
        if num[order2[i]] > 0:
            cls_species_idxs[int_kernel_dim].append(matched[order2[i]])
            cls_coefficients[int_kernel_dim].append(num[order2[i]])
    int_kernel_dim += 1
    _reduce(int_kernel_dim, cls_species_idxs, cls_coefficients, num_species)
    min_value = 1000
    for i in range(len(cls_species_idxs[int_kernel_dim - 1])):
        if not len(int_matched) or all(
            cur_int_matched != cls_species_idxs[int_kernel_dim - 1][i]
            for cur_int_matched in int_matched
        ):
            int_matched.append(cls_species_idxs[int_kernel_dim - 1][i])

        min_value = min(min_value, cls_coefficients[int_kernel_dim - 1][i])
    for i in range(len(cls_species_idxs[int_kernel_dim - 1])):
        cls_coefficients[int_kernel_dim - 1][i] /= min_value

    logger.debug(
        f"Found linearly independent moiety, now there are "
        f"{int_kernel_dim} engaging {len(int_matched)} species"
    )

    return True, int_kernel_dim, int_matched


def _relax(
    stoichiometric_list: Sequence[float],
    int_matched: Sequence[int],
    num_reactions: int,
    num_species: int,
    relaxation_max: float = 1e6,
    relaxation_step: float = 1.9,
) -> bool:
    """Relaxation scheme for Monte Carlo final solution

    Checking for completeness using Motzkin's theorem. See Step (c) in
    De Martino (2014) and the Eqs. 14-16 in the corresponding publication

    :param stoichiometric_list:
        stoichiometric matrix :math:`S` as a flat list (column-major ordering)
    :param int_matched:
        number of matched integer CLs
    :param num_reactions:
        number of reactions in reaction network
    :param num_species:
        number of species in reaction network
    :param relaxation_max:
        maximum relaxation step
    :param relaxation_step:
        relaxation step width
    :returns:
        boolean indicating if relaxation has succeeded (``True``) or not
        (``False``)
    """
    K = len(int_matched)
    matrix: list[list[int]] = [[] for _ in range(K)]
    matrix2: list[list[float]] = [[] for _ in range(K)]
    i_reaction = 0
    i_species = 0
    for val in stoichiometric_list:
        if val != 0:
            take = K
            if K > 0:
                for i in range(K):
                    if i_species == int_matched[i]:
                        take = i
            if take < K:
                matrix[take].append(i_reaction)
                matrix2[take].append(val)
        i_species += 1
        if i_species == num_species:
            i_species = 0
            i_reaction += 1

    # reducing the stoichiometric matrix of conserved moieties to row echelon
    #  form by Gaussian elimination
    order = list(range(K))
    pivots = [matrix[i][0] if len(matrix[i]) else _MAX for i in range(K)]
    done = False
    while not done:
        _qsort(K, 0, order, pivots)
        for j in range(K - 1):
            if pivots[order[j + 1]] == pivots[order[j]] != _MAX:
                min1 = _MAX
                if len(matrix[order[j]]) > 1:
                    for i in range(len(matrix[order[j]])):
                        min1 = min(
                            min1,
                            abs(matrix2[order[j]][0] / matrix2[order[j]][i]),
                        )
                min2 = _MAX
                if len(matrix[order[j + 1]]) > 1:
                    for i in range(len(matrix[order[j + 1]])):
                        min2 = min(
                            min2,
                            abs(
                                matrix2[order[j + 1]][0]
                                / matrix2[order[j + 1]][i]
                            ),
                        )
                if min2 > min1:
                    # swap
                    k2 = order[j + 1]
                    order[j + 1] = order[j]
                    order[j] = k2
        done = True
        for j in range(K - 1):
            if pivots[order[j + 1]] == pivots[order[j]] != _MAX:
                k1 = order[j + 1]
                k2 = order[j]
                column: list[float] = [0] * num_reactions
                g = matrix2[k2][0] / matrix2[k1][0]
                for i in range(1, len(matrix[k1])):
                    column[matrix[k1][i]] = matrix2[k1][i] * g
                for i in range(1, len(matrix[k2])):
                    column[matrix[k2][i]] -= matrix2[k2][i]

                matrix[k1] = []
                matrix2[k1] = []
                for col_idx, col_val in enumerate(column):
                    if abs(col_val) > _MIN:
                        matrix[k1].append(col_idx)
                        matrix2[k1].append(col_val)
                done = False
                if len(matrix[order[j + 1]]):
                    pivots[order[j + 1]] = matrix[order[j + 1]][0]
                else:
                    pivots[order[j + 1]] = _MAX

    # normalize
    for matrix2_i in matrix2:
        if len(matrix2_i):
            norm = matrix2_i[0]
            for j in range(len(matrix2_i)):
                matrix2_i[j] /= norm

    for k1 in reversed(range(K - 1)):
        k = order[k1]
        if len(matrix[k]) <= 1:
            continue

        i = 0
        while i < len(matrix[k]):
            for i_species in range(k1 + 1, K):
                j = order[i_species]
                if not len(matrix[j]) or matrix[j][0] != matrix[k][i]:
                    continue

                # subtract rows
                # matrix2[k] = matrix2[k] - matrix2[j] * matrix2[k][i]
                row_k: list[float] = [0] * num_reactions
                for a in range(len(matrix[k])):
                    row_k[matrix[k][a]] = matrix2[k][a]
                for a in range(len(matrix[j])):
                    row_k[matrix[j][a]] -= matrix2[j][a] * matrix2[k][i]
                # filter
                matrix[k] = [
                    row_idx
                    for row_idx, row_val in enumerate(row_k)
                    if row_val != 0
                ]
                matrix2[k] = [row_val for row_val in row_k if row_val != 0]

                if len(matrix[k]) <= i:
                    break
            i += 1

    indip = [K + 1] * num_reactions
    for i in range(K):
        if len(matrix[i]):
            indip[matrix[i][0]] = i
    M1 = 0
    for i in range(num_reactions):
        if indip[i] == K + 1:
            indip[i] = K + M1
            M1 += 1

    matrixAus = [[] for _ in range(M1)]
    matrixAus2 = [[] for _ in range(M1)]
    i_reaction = 0
    for i in range(num_reactions):
        if indip[i] >= K:
            matrixAus[i_reaction].append(i)
            matrixAus2[i_reaction].append(1)
            i_reaction += 1
        else:
            t = indip[i]
            if len(matrix[t]) > 1:
                for k in range(1, len(matrix[t])):
                    idx = indip[matrix[t][k]] - K
                    matrixAus[idx].append(i)
                    matrixAus2[idx].append(-matrix2[t][k])
    del matrix

    N1 = num_species - K
    matrix_aus = [[] for _ in range(N1)]
    matrix_aus2 = [[] for _ in range(N1)]
    k1 = 0
    i_reaction = 0
    i_species = 0
    for val in stoichiometric_list:
        take = 1
        for i in range(len(int_matched)):
            if i_species == int_matched[i]:
                take -= 1
        if val != 0 and take == 1:
            matrix_aus[k1].append(i_reaction)
            matrix_aus2[k1].append(val)
        i_species += 1
        k1 += take
        if i_species == num_species:
            i_species = 0
            k1 = 0
            i_reaction += 1

    matrixb = [[] for _ in range(N1)]
    matrixb2 = [[] for _ in range(N1)]
    for i in range(M1):
        for j in range(N1):
            if len(matrix_aus[j]) * len(matrixAus[i]):
                prod = 0
                for ib in range(len(matrixAus[i])):
                    for jb in range(len(matrix_aus[j])):
                        if matrixAus[i][ib] == matrix_aus[j][jb]:
                            prod += matrixAus2[i][ib] * matrix_aus2[j][jb]
                if abs(prod) > _MIN:
                    matrixb[j].append(i)
                    matrixb2[j].append(prod)
    del matrixAus, matrixAus2, matrix_aus, matrix_aus2

    var = [_MIN] * M1
    time = 0
    cmin_idx = 0
    while True:
        cmin = 1000
        for j in range(N1):
            constr = 0
            if len(matrixb[j]):
                for i in range(len(matrixb[j])):
                    constr += matrixb2[j][i] * var[matrixb[j][i]]
                if constr < cmin:
                    cmin_idx = j
                    cmin = constr
        if cmin >= 0:
            # constraints satisfied
            break

        # Motzkin relaxation
        alpha = -relaxation_step * cmin
        fact = sum(val**2 for val in matrixb2[cmin_idx])
        alpha /= fact
        alpha = max(1e-9 * _MIN, alpha)
        for j in range(len(matrixb[cmin_idx])):
            var[matrixb[cmin_idx][j]] += alpha * matrixb2[cmin_idx][j]

        time += 1
        if time >= relaxation_max:
            # timeout
            break

    return done


def _reduce(
    int_kernel_dim: int,
    cls_species_idxs: MutableSequence[MutableSequence[int]],
    cls_coefficients: MutableSequence[MutableSequence[float]],
    num_species: int,
) -> None:
    """Reducing the solution which has been found by the Monte Carlo process

    In case of superpositions of independent MCLs one can reduce by
    iteratively subtracting the other independent MCLs, taking care
    to maintain then non-negativity constraint, see Eq. 13 in De Martino (2014)

    :param int_kernel_dim:
        number of found MCLs
    :param cls_species_idxs:
        Species indices involved in each of the conservation laws.
        Modified in-place.
    :param cls_coefficients:
        Coefficients for each of the species involved in each of the
        conservation laws. Modified in-place.
    :param num_species:
        number of species / rows in :math:`S`
    """
    K = int_kernel_dim
    order = list(range(K))
    pivots = [-len(cls_species_idxs[i]) for i in range(K)]

    done = False
    while not done:
        _qsort(K, 0, order, pivots)
        done = True
        for i in range(K - 1):
            k1 = order[i]
            for j in range(i + 1, K):
                k2 = order[j]
                column: list[float] = [0] * num_species
                for species_idx, coefficient in zip(
                    cls_species_idxs[k1], cls_coefficients[k1], strict=True
                ):
                    column[species_idx] = coefficient
                ok1 = True
                for species_idx, coefficient in zip(
                    cls_species_idxs[k2], cls_coefficients[k2], strict=True
                ):
                    column[species_idx] -= coefficient
                    if column[species_idx] < -_MIN:
                        ok1 = False
                        break
                if not ok1:
                    continue

                done = False
                cls_species_idxs[k1] = []
                cls_coefficients[k1] = []
                for col_idx, col_val in enumerate(column):
                    if abs(col_val) > _MIN:
                        cls_species_idxs[k1].append(col_idx)
                        cls_coefficients[k1].append(col_val)
                pivots[k1] = -len(cls_species_idxs[k1])
