class CC(BaseQuantifier):
  """Classify and Count (CC) quantifier.

  The CC method is the simplest baseline in quantification. It works by
  classifying all unlabeled instances in the test bag using a hard (crisp)
  classifier, and then computing the relative frequency (prevalence) of
  each class based on those predictions.

  Parameters
  ----------
  classifier : estimator object, default = None
    The classifier to be used as the base for quantification.
    If None, an instance of `LogisticRegression()` will be created.

  Attributes
  ----------
  classes_ : ndarray of shape (n_classes,)
    The distinct class labels found during training, sorted ascending
    (as returned by `np.unique`).

  n_classes_ : int
    The number of distinct classes.

  train_prevalence_ : ndarray of shape (n_classes,)
    The prevalence of each class in the training dataset.

  y_prevs_ : ndarray of shape (n_classes,)
    Alias for train_prevalence_ kept for backward compatibility.

  classifier_ : estimator object
    The fitted base classifier trained on the entire dataset.

  Notes
  -----
  The Classify and Count method does not adjust for misclassifications
  made by the base classifier (false positives and false negatives). Therefore,
  its performance is highly dependent on the classification accuracy. 

  The estimated prevalence for class 'c' is given by the formula:

    p_hat(c) = (1 / |X|) * sum( I( f(x) == c ) )

  where f(x) is the crisp prediction of the classifier for instance x,
  and I() is the indicator function.

  References
  ----------
  George Forman. Counting positives accurately despite inaccurate classification.
  In Proceedings of the 16th European Conference on Machine Learning, pages 564-575,
  Porto, Portugal, 2005.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = CC()
  >>> quantifier.fit(X, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """

  def __init__(self, classifier: BaseEstimator = None):
    super().__init__(classifier)

  def fit(self, X: np.ndarray, y: np.ndarray) -> "CC":
    """Adjusts the CC quantifier by fitting the base classifier.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
      Training data.
    y : array-like of shape (n_samples,)
      Labels for the corresponding classes.

    Returns
    -------
    self : object
        Returns the fitted estimator instance itself.
    """
    X, y = check_X_y(X, y, accept_sparse=True)

    self.classes_, counts = np.unique(y, return_counts=True)
    self.n_classes_ = len(self.classes_)
    self.train_prevalence_ = counts / len(y)
    self.y_prevs_ = self.train_prevalence_ # Compatibility purposes
    # lazy validation of the classifier
    base_classifier = self.classifier if self.classifier is not None else LogisticRegression()
    self.classifier_ = clone(base_classifier)
    self.classifier_.fit(X, y) # fit the classifier with all training data

    return self

  def predict(self, X: np.ndarray) -> np.ndarray:
    """Estimate the class prevalences for the test bag X.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
      The test bag with unlabelled instances.

    Returns
    -------
    p_adjusted : ndarray of shape (n_classes,)
      An array with the estimated prevalences for each class,
      normalized to sum up to 1.0.
    """
    check_is_fitted(self)
    X = check_array(X, accept_sparse=True)

    y_pred = self.classifier_.predict(X)

    # fully vectorized counting: since self.classes_ is sorted ascending
    # (np.unique's guarantee), searchsorted maps every prediction to its
    # position in self.classes_ and bincount tallies them in one pass,
    # avoiding the per-call Python dict construction/lookup this replaces
    class_indices = np.searchsorted(self.classes_, y_pred)
    counts = np.bincount(class_indices, minlength=self.n_classes_).astype(float)

    return normalize_prevalence(counts, self.n_classes_)

class ACC(BaseCalibratedQuantifier):
  """Adjusted Classify and Count (ACC) quantifier for binary problems.

  ACC is a calibrated quantification method that adjusts the simple Classify 
  and Count (CC) estimates. It uses Out-of-Fold (OOF) predictions to calculate 
  the True Positive Rate (TPR) and False Positive Rate (FPR) of the base classifier. 
  Then, it mathematically corrects the final prevalence using the classical 
  binary adjustment formula.

  Parameters
  ----------
  classifier : estimator object, default=None
    The classifier to be used as the base for quantification.
    If None, an instance of `LogisticRegression()` will be created.
  cv : int, cross-validation generator or an iterable, default = 10
    Determines the cross-validation splitting strategy to generate the 
    Out-of-Fold predictions used for calibration.
  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds (plus
    the final full-data classifier refit). See `BaseCalibratedQuantifier`.
  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs. See
    `BaseCalibratedQuantifier`.

  Attributes
  ----------
  classes_ : ndarray of shape (2,)
    The distinct class labels found during training.

  n_classes_ : int
    The number of distinct classes (expected to be 2).

  train_prevalence_ : ndarray of shape (2,)
    The prevalence of each class in the training dataset.

  tpr_ : float
    The True Positive Rate (Sensitivity) estimated out-of-fold.

  fpr_ : float
    The False Positive Rate (1 - Specificity) estimated out-of-fold.

  classifier_ : estimator object
    The fitted base classifier trained on the entire dataset.

  Notes
  -----
  The ACC method adjusts the observed crisp prevalence (p_pred) using the 
  following analytical equation for the positive class:

    p_adjusted = (p_pred - fpr) / (tpr - fpr)

  The prevalence for the negative class is then derived as 1 - p_corrected.
  If TPR equals FPR, the denominator becomes zero, and the method falls back 
  to the unadjusted p_pred.

  References
  ----------
  George Forman. Counting positives accurately despite inaccurate classification.
  In Proceedings of the 16th European Conference on Machine Learning, pages 564-575,
  Porto, Portugal, 2005.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = ACC()
  >>> quantifier.fit(X, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """

  def __init__(self, classifier: BaseEstimator = None, cv: int = 10,
               n_jobs: int = None, parallel_backend: str = "loky"):
    super().__init__(classifier=classifier, cv=cv, n_jobs=n_jobs, parallel_backend=parallel_backend)

  def _get_oof_method(self) -> str:
    return "predict"

  def _calibrate(self, y_true_oof: np.ndarray, y_pred_oof: np.ndarray):
    """Calculates TPR and FPR using the Out-of-Fold predictions.

    Parameters
    ----------
    y_true_oof : ndarray of shape (n_samples,)
      True labels collected out-of-fold.
    y_pred_oof : ndarray of shape (n_samples,)
      Crisp predictions generated out-of-fold.
    """
    # maps each class using the scikit-learn pattern (asc format)
    neg_label = self.classes_[0]
    pos_label = self.classes_[1]

    is_true_pos = (y_true_oof == pos_label)
    is_true_neg = (y_true_oof == neg_label)

    tp = np.sum(is_true_pos & (y_pred_oof == pos_label))
    fn = np.sum(is_true_pos & (y_pred_oof == neg_label))
    fp = np.sum(is_true_neg & (y_pred_oof == pos_label))
    tn = np.sum(is_true_neg & (y_pred_oof == neg_label))

    self.tpr_ = tp / (tp + fn) if (tp + fn) > 0 else 1.0
    self.fpr_ = fp / (fp + tn) if (fp + tn) > 0 else 0.0

  def _quantify(self, X: np.ndarray) -> np.ndarray:
    """Applies the binary ACC equation to adjust the predicted prevalence.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features)
      Raw data from test bags.

    Returns
    -------
    p_acc : ndarray of shape (2,)
      The adjusted prevalences for [negative_class, positive_class].
    """
    y_pred = self.classifier_.predict(X)
    pos_label = self.classes_[1]

    p_pred = np.mean(y_pred == pos_label)

    denominator = self.tpr_ - self.fpr_
    if denominator == 0:
      p_pos_adjusted = p_pred
    else:
      p_pos_adjusted = (p_pred - self.fpr_) / denominator

    # guardrails: clipping
    p_pos_adjusted = np.clip(p_pos_adjusted, 0.0, 1.0)
    p_neg_adjusted = 1.0 - p_pos_adjusted

    return np.array([p_neg_adjusted, p_pos_adjusted])

  def fit(self, X: np.ndarray, y: np.ndarray) -> "ACC":
    """Fits the ACC quantifier, enforcing the binary-only constraint.

    Raises
    ------
    ValueError
      If `y` contains more than 2 distinct classes.
    """
    if len(np.unique(y)) > 2:
      raise ValueError(
        "ACC method only works for binary quantification. Multiclass "
        "quantification is possible via OVR strategies, but not recommended due to "
        "theoretical issues with that approach."
      )

    return super().fit(X, y)

class PCC(BaseQuantifier):
  """Probabilistic Classify and Count (PCC) quantifier.

  The PCC method extends the classical Classify and Count (CC) by using the
  posterior probabilities generated by a classifier instead of its crisp (hard) 
  predictions. The estimated prevalence for each class is computed as the 
  expected value (average) of the predicted probabilities across all instances 
  in the test bag.

  Parameters
  ----------
  classifier : estimator object, default = None
    The classifier to be used as the base for quantification.
    Must implement the `predict_proba` method. If None, an instance of 
    `LogisticRegression()` will be created.

  Attributes
  ----------
  classes_ : ndarray of shape (n_classes,)
    The distinct class labels found during training.

  n_classes_ : int
    The number of distinct classes.

  train_prevalence_ : ndarray of shape (n_classes,)
    The prevalence of each class in the training dataset.

  y_prevs_ : ndarray of shape (n_classes,)
    Alias for train_prevalence_ kept for backward compatibility.

  classifier_ : estimator object
    The fitted base classifier trained on the entire dataset.

  Notes
  -----
  Unlike CC, PCC accounts for the classifier's confidence in its predictions, 
  which often leads to better quantification performance, especially when 
  the classifier is well-calibrated. However, like CC, it is an "unadjusted" 
  method, meaning it does not explicitly correct for systematic classification 
  errors via a confusion matrix.

  The estimated prevalence for class 'c' is given by the formula:

    p_hat(c) = (1 / |X|) * sum( P(c | x) )

  where P(c | x) is the probability that instance x belongs to class c, 
  estimated by the base classifier's `predict_proba` method.

  References
  ----------
  Antonio Bella, Cesar Ferri, José Hernández-Orallo, and María José Ramírez-Quintana.
  Quantification via probability estimators. In 2010 IEEE International Conference on
  Data Mining, pages 737-742, Sydney, Australia, 2010.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = PCC()
  >>> quantifier.fit(X, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """
  def __init__(self, classifier: BaseEstimator = None):
    super().__init__(classifier=classifier)

  def fit(self, X: np.ndarray, y: np.ndarray) -> "PCC":
    """Adjusts the PCC quantifier by fitting the probabilistic base classifier.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
      Training data.
    y : array-like of shape (n_samples,)
      Labels for the corresponding classes.

    Returns
    -------
    self : object
      Returns the fitted estimator instance itself.

    Raises
    ------
    TypeError
      If the provided classifier does not support probability estimation 
      (i.e., lacks a `predict_proba` method).
    """
    X, y = check_X_y(X, y, accept_sparse=True)

    self.classes_, counts = np.unique(y, return_counts=True)
    self.n_classes_ = len(self.classes_)
    self.train_prevalence_ = counts / len(y)
    self.y_prevs_ = self.train_prevalence_ # Compatibility purposes

    base_classifier = self.classifier if self.classifier is not None else LogisticRegression()

    if not hasattr(base_classifier, "predict_proba"):
      raise TypeError(
        f"The classifier {base_classifier.__class__.__name__} does not "
        "support probability estimation. PCC requires 'predict_proba'."
      )

    self.classifier_ = clone(base_classifier)
    self.classifier_.fit(X, y)     

    return self

  def predict(self, X: np.ndarray) -> np.ndarray:
    """Estimate the class prevalences for the test bag X using soft probabilities.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
      The test bag with unlabelled instances.

    Returns
    -------
    p_adjusted : ndarray of shape (n_classes,)
      An array with the estimated prevalences for each class,
      normalized to sum up to 1.0.
    """
    check_is_fitted(self)
    X = check_array(X, accept_sparse=True)
    # obtain the soft probability matrix (n_samples, n_classes)
    y_probas = self.classifier_.predict_proba(X)

    p_pred = np.mean(y_probas, axis=0)
    return normalize_prevalence(p_pred, self.n_classes_)

class PACC(BaseCalibratedQuantifier):
  """Probabilistic Adjusted Classify and Count (PACC) quantifier for binary problems.

  PACC is a calibrated quantification method that refines the Probabilistic 
  Classify and Count (PCC) by correcting for the base classifier's systematic 
  probabilistic bias. Using Out-of-Fold (OOF) predictions, it computes the 
  expected predicted probability for the positive class given the true class labels. 
  It then applies an analytical correction formula similar to ACC but entirely 
  based on continuous probability scores.

  Parameters
  ----------
  classifier : estimator object, default = None
    The classifier to be used as the base for quantification.
    Must implement the `predict_proba` method. If None, an instance of 
    `LogisticRegression()` will be created.
  cv : int, cross-validation generator or an iterable, default = 10
    Determines the cross-validation splitting strategy to generate the 
    Out-of-Fold predictions used for calibration.  
  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds (plus
    the final full-data classifier refit). See `BaseCalibratedQuantifier`.
  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs. See
    `BaseCalibratedQuantifier`.

  Attributes
  ----------
  classes_ : ndarray of shape (2,)
    The distinct class labels found during training.

  n_classes_ : int
    The number of distinct classes (expected to be 2).

  train_prevalence_ : ndarray of shape (2,)
    The prevalence of each class in the training dataset.

  y_prevs_ : ndarray of shape (2,)
    Alias for train_prevalence_ kept for backward compatibility.

  mu_pos_pos_ : float
    The mean probability assigned to the positive class for instances that 
    are truly positive, estimated out-of-fold. Analogue to TPR.

  mu_neg_pos_ : float
    The mean probability assigned to the positive class for instances that 
    are truly negative, estimated out-of-fold. Analogue to FPR.

  classifier_ : estimator object
    The fitted base classifier trained on the entire dataset.

  Notes
  -----
  The PACC method adjusts the observed probabilistic prevalence (p_pcc) using 
  the following analytical equation for the positive class:

    p_corrected = (p_pcc - mu_neg_pos) / (mu_pos_pos - mu_neg_pos)

  where 'p_pcc' is the average predicted probability of the positive class 
  in the test bag. The prevalence for the negative class is then derived as 
  1 - p_corrected. If the denominator evaluates to zero, the method falls 
  back to the unadjusted p_pcc.

  References
  ----------
  Antonio Bella, Cesar Ferri, José Hernández-Orallo, and María José Ramírez-Quintana.
  Quantification via probability estimators. In 2010 IEEE International Conference on
  Data Mining, pages 737-742, Sydney, Australia, 2010.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = PACC()
  >>> quantifier.fit(X, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """

  def __init__(self, classifier: BaseEstimator = None, cv: int = 10,
               n_jobs: int = None, parallel_backend: str = "loky"):
    super().__init__(classifier=classifier, cv=cv, n_jobs=n_jobs, parallel_backend=parallel_backend)

  def fit(self, X: np.ndarray, y: np.ndarray) -> "PACC":
    """Fits the PACC quantifier, validating classifier and class-count constraints.

    Raises
    ------
    TypeError
      If the provided classifier does not support probability estimation.
    ValueError
      If `y` contains more than 2 distinct classes.
    """
    base_clf = self.classifier if self.classifier is not None else LogisticRegression()

    if not hasattr(base_clf, "predict_proba"):
      raise TypeError(
        f"The classifier {base_clf.__class__.__name__} does not "
        "support probability estimation. PACC requires 'predict_proba'."
      )

    if len(np.unique(y)) > 2:
      raise ValueError(
        "PACC method only works for binary quantification. Multiclass "
        "quantification is possible via OVR strategies, but not recommended due to "
        "theoretical issues with that approach."
      )

    return super().fit(X, y)

  def _get_oof_method(self) -> str:
    return "predict_proba"

  def _calibrate(self, y_true_oof: np.ndarray, y_pred_oof: np.ndarray):
    """Calculates the expected continuous scores (mu) using Out-of-Fold probabilities.

    Parameters
    ----------
    y_true_oof : ndarray of shape (n_samples,)
      True labels collected out-of-fold.
    y_pred_oof : ndarray of shape (n_samples, 2)
      Continuous probability matrices generated out-of-fold.
    """
    neg_label = self.classes_[0]
    pos_label = self.classes_[1]

    prob_pos_oof = y_pred_oof[:, 1]

    is_truly_pos = (y_true_oof == pos_label)
    is_truly_neg = (y_true_oof == neg_label)

    # mu_neg: probabilities mean attributed to the negative class
    self.mu_pos_pos_ = np.mean(prob_pos_oof[is_truly_pos]) if np.any(is_truly_pos) else 1.0
    self.mu_neg_pos_ = np.mean(prob_pos_oof[is_truly_neg]) if np.any(is_truly_neg) else 0.0

  def _quantify(self, X: np.ndarray) -> np.ndarray:
    """Applies the analytical continuous equation to adjust the probabilistic prevalence.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features)
      Raw data from test bags.

    Returns
    -------
    p_pacc : ndarray of shape (2,)
      The adjusted prevalences for [negative_class, positive_class].
    """
    y_probas = self.classifier_.predict_proba(X)
    p_pcc = np.mean(y_probas[:, 1])

    denominator = self.mu_pos_pos_ - self.mu_neg_pos_

    if denominator == 0:
      p_pos_adjusted = p_pcc
    else:
      p_pos_adjusted = (p_pcc - self.mu_neg_pos_) / denominator

    return np.array([1.0 - p_pos_adjusted, p_pos_adjusted])

class GAC(BaseCalibratedQuantifier, BaseMixtureQuantifier):
  """Generalized Adjusting Confusion Matrix (GAC) Quantifier.

  A distance-minimizing multi-class generalization of the Adjusting Count (AC) 
  algorithm that optimizes target distributions across the discrete labels 
  confusion matrix profile.

  Parameters
  ----------
  classifier : estimator object, default = LogisticRegression
    The underlying base classifier. Defaults to `LogisticRegression()`.

  distance_metric : str, default = 'L2'
    The distance metric minimized.

  cv : int, default = 10
    The number of cross-validation folds.

  use_convex_solver : bool, default = True
    If True, optimizes via CVXPY.

  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds.

  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs.

  References
  ----------
  Aykut Firat. Unified framework for quantification. arXiv preprint arXiv:1606.00868, 2016.
  """

  def __init__(self,
               classifier: BaseEstimator = LogisticRegression(),
               distance_metric: str = "L2", 
               cv: int = 10,
               use_convex_solver: bool = True,
               n_jobs: int = None,
               parallel_backend: str = "loky"):
    BaseCalibratedQuantifier.__init__(self, classifier=classifier, cv=cv,
                                      n_jobs=n_jobs, parallel_backend=parallel_backend)
    BaseMixtureQuantifier.__init__(self, classifier=classifier, distance_metric=distance_metric, 
                                  use_convex_solver=use_convex_solver)

  def _get_oof_method(self) -> str:
    return "predict"

  def _calibrate(self, y_true_oof: np.ndarray, y_pred_oof: np.ndarray):
    # fully vectorized confusion matrix: one-hot encode both true and
    # predicted labels against the sorted classes_ grid, then a single
    # matmul (n_classes, n_samples) @ (n_samples, n_classes) tallies
    # every (pred, true) pair at once, replacing the double Python loop
    # this used to run over classes_ x classes_
    true_idx = np.searchsorted(self.classes_, y_true_oof)
    pred_idx = np.searchsorted(self.classes_, y_pred_oof)
    one_hot_true = np.eye(self.n_classes_)[true_idx]
    one_hot_pred = np.eye(self.n_classes_)[pred_idx]
    confusion_matrix = one_hot_pred.T @ one_hot_true

    _, class_counts = np.unique(y_true_oof, return_counts=True)
    self.conditional_matrix_ = confusion_matrix / class_counts

  def _compute_score(self, X: np.ndarray) -> np.ndarray:
    y_predictions = self.classifier_.predict(X)
    return np.array([np.mean(y_predictions == class_label) for class_label in self.classes_])

  def _quantify(self, X: np.ndarray) -> np.ndarray:
    test_frequencies = self._compute_score(X)
    return self._solve_mixture(test_frequencies)

class GPAC(BaseCalibratedQuantifier, BaseMixtureQuantifier):
  """Generalized Probabilistic Adjusting Confusion Matrix (GPAC) Quantifier.

  A distance-minimizing multi-class generalization of the Probabilistic Adjusting 
  Count (PAC) algorithm that optimizes target distributions using soft-probability profiles.

  Parameters
  ----------
  classifier : estimator object, default = LogisticRegression
    The underlying base classifier. Defaults to `LogisticRegression()`.

  distance_metric : str, default = 'L2'
    The distance metric minimized.

  cv : int, default = 10
    The number of cross-validation folds.

  use_convex_solver : bool, default = True
    If True, optimizes via CVXPY.

  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds.

  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs.

  References
  ----------
  Aykut Firat. Unified framework for quantification. arXiv preprint arXiv:1606.00868, 2016.
  """

  def __init__(self,
               classifier: BaseEstimator = LogisticRegression(),
               distance_metric: str = "L2", 
               cv: int = 10,
               use_convex_solver: bool = True,
               n_jobs: int = None,
               parallel_backend: str = "loky"):
    BaseCalibratedQuantifier.__init__(self, classifier=classifier, cv=cv,
                                      n_jobs=n_jobs, parallel_backend=parallel_backend)
    BaseMixtureQuantifier.__init__(self, classifier=classifier, distance_metric=distance_metric, 
                                  use_convex_solver=use_convex_solver)

  def _get_oof_method(self) -> str:
    return "predict_proba"

  def _calibrate(self, y_true_oof: np.ndarray, y_pred_oof: np.ndarray):
    # single matmul replaces the loop over classes_: column l of the
    # result is the sum of predicted-probability rows whose true label is
    # class l, i.e. y_pred_oof.T @ one_hot(true_labels)
    true_idx = np.searchsorted(self.classes_, y_true_oof)
    one_hot_true = np.eye(self.n_classes_)[true_idx]
    probabilistic_matrix = y_pred_oof.T @ one_hot_true

    _, class_counts = np.unique(y_true_oof, return_counts=True)
    self.conditional_matrix_ = probabilistic_matrix / class_counts

  def _compute_score(self, X: np.ndarray) -> np.ndarray:
    return self.classifier_.predict_proba(X).sum(axis=0) / X.shape[0]

  def _quantify(self, X: np.ndarray) -> np.ndarray:
    test_frequencies = self._compute_score(X)
    return self._solve_mixture(test_frequencies)

class FM(BaseCalibratedQuantifier, BaseMixtureQuantifier):
  """Friedman's Method (FM) Quantifier.

  An adjusting prediction mixture model that maps soft classifier probabilities 
  into binary indicator matrices by comparing them against baseline training priors.

  Parameters
  ----------
  classifier : estimator object, default = LogisticRegression
    The underlying base classifier. Defaults to `LogisticRegression()`.

  distance_metric : str, default = 'L2'
    The distance metric minimized.

  cv : int, default = 10
    The number of cross-validation folds.

  use_convex_solver : bool, default = True
    If True, optimizes via CVXPY.

  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds.

  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs.

  References
  ----------
  Jerome H. Friedman. Class counts in future unlabeled samples, 2014.
  Presentation at MIT CSAIL Big Data Event.
  """

  def __init__(self,
               classifier: BaseEstimator = LogisticRegression(),
               distance_metric: str = "L2",
               cv: int = 10,
               use_convex_solver: bool = True,
               n_jobs: int = None,
               parallel_backend: str = "loky"):
    BaseCalibratedQuantifier.__init__(self, classifier=classifier, cv=cv,
                                      n_jobs=n_jobs, parallel_backend=parallel_backend)
    BaseMixtureQuantifier.__init__(self, classifier=classifier, distance_metric=distance_metric, 
                                  use_convex_solver=use_convex_solver)

  def _get_oof_method(self) -> str:
    return "predict_proba"

  def _calibrate(self, y_true_oof: np.ndarray, y_pred_oof: np.ndarray):
    # single matmul replaces the loop over classes_, mirroring GPAC's
    # vectorization but comparing against train_prevalence_ first
    true_idx = np.searchsorted(self.classes_, y_true_oof)
    one_hot_true = np.eye(self.n_classes_)[true_idx]
    above_prior = (y_pred_oof > self.train_prevalence_).astype(float)
    threshold_matrix = above_prior.T @ one_hot_true

    _, class_counts = np.unique(y_true_oof, return_counts=True)
    self.conditional_matrix_ = threshold_matrix / class_counts

  def _compute_score(self, X: np.ndarray) -> np.ndarray:
    return np.sum(self.classifier_.predict_proba(X) > self.train_prevalence_, axis=0) / X.shape[0]

  def _quantify(self, X: np.ndarray) -> np.ndarray:
    test_frequencies = self._compute_score(X)
    return self._solve_mixture(test_frequencies)

class X(BaseThresholdQuantifier):
  """Forman's X threshold selection quantifier matching QFY's TSX strategy.

  This method selects the optimal threshold from the dynamic grid that 
  minimizes the absolute distance |TPR - (1 - FPR)|, searching for the 
  intersection point in the ROC curve where TPR + FPR ~= 1.0.

  Parameters
  ----------
  classifier : estimator object, default = None
    The classifier to be used as the base for quantification.
  cv : int, default = 10
    Determines the cross-validation splitting strategy.
  precision : int, default = 3
    The decimal precision used to round the Out-of-Fold probabilities.
  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds. See
    `BaseCalibratedQuantifier`.
  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs.

  References
  ----------
  Forman, G. (2008). Quantifying counts and costs via classification. 
  Data Mining and Knowledge Discovery, 17(2):164-206.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X_, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = X()
  >>> quantifier.fit(X_, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """

  def _quantify(self, X_test: np.ndarray) -> np.ndarray:
    # minimizes the absolute distance |TPR - (1 - FPR)|
    idx = np.argmin(np.abs(self.tpr_by_thresh_ - (1.0 - self.fpr_by_thresh_)))

    t_chosen = self.thresholds_[idx]
    tpr = self.tpr_by_thresh_[idx]
    fpr = self.fpr_by_thresh_[idx]

    pos_probs_test = self.classifier_.predict_proba(X_test)[:, 1]
    p_raw_pos = np.mean(pos_probs_test >= t_chosen)

    p_adj_pos = self._apply_acc_formula(p_raw_pos, tpr, fpr)
    return np.array([1.0 - p_adj_pos, p_adj_pos])

class Max(BaseThresholdQuantifier):
  """Forman's Max threshold selection quantifier matching QFY's TSMax strategy.

  This method selects the threshold that maximizes the separation split 
  between classes, maximizing the absolute difference |TPR - FPR| to achieve 
  the highest denominator stability during shift adjustment.

  Parameters
  ----------
  classifier : estimator object, default = None
    The classifier to be used as the base for quantification.
  cv : int, default = 10
    Determines the cross-validation splitting strategy.
  precision : int, default = 3
    The decimal precision used to round the Out-of-Fold probabilities.
  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds. See
    `BaseCalibratedQuantifier`.
  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs.

  References
  ----------
  Forman, G. (2008). Quantifying counts and costs via classification. 
  Data Mining and Knowledge Discovery, 17(2):164-206.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = Max()
  >>> quantifier.fit(X, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """

  def _quantify(self, X_test: np.ndarray) -> np.ndarray:
    # maximizes the separation split |TPR - FPR|
    idx = np.argmax(np.abs(self.tpr_by_thresh_ - self.fpr_by_thresh_))

    t_chosen = self.thresholds_[idx]
    tpr = self.tpr_by_thresh_[idx]
    fpr = self.fpr_by_thresh_[idx]

    pos_probs_test = self.classifier_.predict_proba(X_test)[:, 1]
    p_raw_pos = np.mean(pos_probs_test >= t_chosen)

    p_adj_pos = self._apply_acc_formula(p_raw_pos, tpr, fpr)
    return np.array([1.0 - p_adj_pos, p_adj_pos])

class T50(BaseThresholdQuantifier):
  """Forman's Threshold 50 quantifier matching QFY's TS50 strategy.

  This method selects the threshold that is closest to achieving a 
  True Positive Rate (TPR) of exactly 0.5, minimizing the cost function 
  |TPR - 0.5| across the candidate grid.

  Parameters
  ----------
  classifier : estimator object, default = None
    The classifier to be used as the base for quantification.
  cv : int, default = 10
    Determines the cross-validation splitting strategy.
  precision : int, default = 3
    The decimal precision used to round the Out-of-Fold probabilities.
  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds. See
    `BaseCalibratedQuantifier`.
  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs.

  References
  ----------
  Forman, G. (2008). Quantifying counts and costs via classification. 
  Data Mining and Knowledge Discovery, 17(2):164-206.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = T50()
  >>> quantifier.fit(X, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """

  def _quantify(self, X_test: np.ndarray) -> np.ndarray:
    # minimizes abs(tpr - 0.5)
    idx = np.argmin(np.abs(self.tpr_by_thresh_ - 0.5))

    t_chosen = self.thresholds_[idx]
    tpr = self.tpr_by_thresh_[idx]
    fpr = self.fpr_by_thresh_[idx]

    pos_probs_test = self.classifier_.predict_proba(X_test)[:, 1]
    p_raw_pos = np.mean(pos_probs_test >= t_chosen)

    p_adj_pos = self._apply_acc_formula(p_raw_pos, tpr, fpr)
    return np.array([1.0 - p_adj_pos, p_adj_pos])

class MedianSweep(BaseThresholdQuantifier):
  """Forman's Median Sweep (MS) quantifier faithful to QFY's implementation.

  Median Sweep evaluates individual ACC adjustments across all candidates in the 
  dynamic grid. It enforces a strict filter where thresholds with a denominator
  (TPR - FPR) smaller than `delta_min` are ignored.
  If no thresholds satisfy the filter, it falls back to the estimate with the 
  largest available separation. Otherwise, it returns the robust median of all 
  valid predictions.

  Parameters
  ----------
  classifier : estimator object, default=None
    The classifier to be used as the base for quantification.
  cv : int, default=10
    Determines the cross-validation splitting strategy.
  precision : int, default=3
    The decimal precision used to round the Out-of-Fold probabilities.
  delta_min : float, default=0.25
    The minimum required threshold separation (TPR - FPR) to accept an 
    individual ACC adaptation, avoiding denominator instability.
  n_jobs : int, default = None
    Number of jobs to run in parallel while fitting the `cv` folds. See
    `BaseCalibratedQuantifier`.
  parallel_backend : str, default = "loky"
    `joblib.Parallel` backend used for the CV/final-fit jobs.

  Notes
  -----
  The fallback path (used only when **no** threshold satisfies
  `delta > delta_min`) picks the threshold with the largest `delta` among
  *all* candidates: as long as no threshold is ever accepted as valid,
  every candidate is a fallback contender, so the largest-`delta`
  candidate found while sweeping the whole grid is exactly the one a
  sequential scan would have kept until the end.

  References
  ----------
  Forman, G. (2008). Quantifying counts and costs via classification. 
  Data Mining and Knowledge Discovery, 17(2):164-206.

  Examples
  --------
  >>> from sklearn.datasets import make_classification
  >>> X, y = make_classification(n_samples=1000, n_classes=2, random_state=42)
  >>> quantifier = MedianSweep()
  >>> quantifier.fit(X, y)
  >>> X_test, _ = make_classification(n_samples=200, n_classes=2, random_state=7)
  >>> prevalences = quantifier.predict(X_test)
  >>> print(prevalences)
  """

  def __init__(self, classifier: BaseEstimator = None, cv: int = 10, precision: int = 3,
               delta_min: float = 0.25, n_jobs: int = None, parallel_backend: str = "loky"):
    super().__init__(classifier=classifier, cv=cv, precision=precision,
                     n_jobs=n_jobs, parallel_backend=parallel_backend)
    self.delta_min = delta_min

  def _quantify(self, X_test: np.ndarray) -> np.ndarray:
    pos_probs_test = self.classifier_.predict_proba(X_test)[:, 1]
    delta = self.tpr_by_thresh_ - self.fpr_by_thresh_

    # vectorized ACC adjustment for every threshold at once: broadcasting
    # (n_test, n_thresholds) comparisons, then averaging each column gives
    # the raw positive rate at every threshold in one shot
    p_raw_pos_by_thresh = np.mean(
      pos_probs_test[:, np.newaxis] >= self.thresholds_[np.newaxis, :], axis=0
    )

    valid_mask = delta > self.delta_min

    if np.any(valid_mask):
      p_adj_valid = (p_raw_pos_by_thresh[valid_mask] - self.fpr_by_thresh_[valid_mask]) / delta[valid_mask]
      p_adj_pos = np.median(np.clip(p_adj_valid, 0.0, 1.0))
    else:
      # no threshold cleared delta_min: fall back to the single threshold
      # with the largest available separation across the whole grid
      # (equivalent to the sequential scan's behavior in this scenario,
      # see the class-level Notes section)
      idx = np.argmax(delta)
      if delta[idx] == 0:
        p_max_fallback = self.tpr_by_thresh_[idx]
      else:
        p_max_fallback = (p_raw_pos_by_thresh[idx] - self.fpr_by_thresh_[idx]) / delta[idx]
      p_adj_pos = np.clip(p_max_fallback, 0.0, 1.0)

    return np.array([1.0 - p_adj_pos, p_adj_pos])

class HDx(BaseMixtureQuantifier):
  """Hellinger Distance x (HDx) quantifier.

  HDx is a non-parametric feature-space mixture model that operates directly on 
  categorical or discretized continuous features without training an underlying 
  classifier. It projects each feature column independently, constructs a global 
  marginal conditional probability matrix during training, and minimizes the 
  Hellinger Distance to estimate the test class prevalences.

  Parameters
  ----------
  use_convex_solver : bool, default=True
    If True, attempts to solve the statistical mixture distribution using `cvxpy`.
    If False, falls back to the Golden Section Search numerical solver.

  Attributes
  ----------
  classes_ : ndarray of shape (n_classes,)
    The distinct class labels found during the training phase.

  n_classes_ : int
    The total number of unique classes.

  train_prevalence_ : ndarray of shape (n_classes,)
    The baseline prevalence proportion of each class observed in the training data.

  feature_spaces_ : list of ndarray
    A list of length `n_features`, where each element contains the unique sorted 
    values observed for that specific feature column during training.

  conditional_matrix_ : ndarray of shape (n_total_unique_bins, n_classes)
    The stacked conditional probability matrix built during the `fit` phase.
    Represents the marginal distribution profiles for each class.

  References
  ----------
  Víctor González-Castro, Rocío Alaiz-Rodríguez, and Enrique Alegre. Class distribution
  estimation based on the Hellinger distance. Information Sciences, 218(1):146-164, 2013
  """
  def __init__(self, use_convex_solver: bool = True):
    # HDx operates on features directly (classifier=None) and strictly uses Hellinger Distance ("HD")
    super().__init__(classifier=None,
                     distance_metric="HD",
                     use_convex_solver=use_convex_solver)
    self.feature_spaces_ = None

  def fit(self, X: np.ndarray, y: np.ndarray) -> 'HDx':
    """Fits the HDx mixture model by building the marginal conditional matrix.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features)
      The training feature matrix.
    y : ndarray of shape (n_samples,)
      The target class labels.

    Returns
    -------
    self : object
      Returns the instance itself.
    """
    X, y = check_X_y(X, y, accept_sparse=False)

    self.classes_, class_counts = np.unique(y, return_counts=True)
    self.n_classes_ = len(self.classes_)
    self.train_prevalence_ = class_counts / len(y)

    if self.n_classes_ < 2:
      raise ValueError("HDx requires at least 2 distinct classes to fit.")

    n_features = X.shape[1]
    class_idx = np.searchsorted(self.classes_, y)

    # map and store the unique token space for each individual feature column
    self.feature_spaces_ = [np.unique(X[:, j]) for j in range(n_features)]

    # build the conditional matrix (CM): the loop over features is
    # unavoidable (each column has a different number of unique values,
    # so the blocks can't be stacked into a single rectangular operation),
    # but within each feature the value x class crosstab is fully
    # vectorized via a combined-index bincount instead of the previous
    # nested (class x unique_value) Python loop
    conditional_blocks = []
    for j in range(n_features):
      unique_values = self.feature_spaces_[j]
      val_idx = np.searchsorted(unique_values, X[:, j])

      combined_idx = val_idx * self.n_classes_ + class_idx
      counts_flat = np.bincount(combined_idx, minlength=len(unique_values) * self.n_classes_)
      crosstab_counts = counts_flat.reshape(len(unique_values), self.n_classes_)

      # normalize counts by each class size to form conditional probabilities
      conditional_blocks.append(crosstab_counts / class_counts)

    # vertically stack all independent column representations into a single global system matrix
    self.conditional_matrix_ = np.vstack(conditional_blocks)

    return self

  def _compute_score(self, X: np.ndarray) -> np.ndarray:
    """Extracts the empirical marginal test frequencies across all features.

    Calculates the relative sample distribution frequency over the saved 
    training feature spaces for the incoming test batch.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features)
      The testing feature matrix.

    Returns
    -------
    test_frequencies : ndarray of shape (n_total_unique_bins,)
      The stacked vector representing empirical marginal frequencies of the test batch.
    """
    n_samples = X.shape[0]
    n_features = X.shape[1]
    frequencies_list = []

    # compute marginal frequencies matching the exact bins established
    # during training; the per-feature loop remains for the same reason
    # as in fit (ragged bin counts across columns), but the inner
    # per-unique-value counting is now a single searchsorted + bincount
    # pass instead of one np.count_nonzero comparison per unique value.
    # Test values absent from the trained feature space (i.e. not an
    # exact match to any trained unique value) contribute 0, exactly
    # matching the original np.count_nonzero(X[:, j] == val) semantics.
    for j in range(n_features):
      unique_values = self.feature_spaces_[j]
      col = X[:, j]

      idx = np.searchsorted(unique_values, col)
      idx_clipped = np.clip(idx, 0, len(unique_values) - 1)
      exact_match = unique_values[idx_clipped] == col

      counts = np.bincount(idx_clipped[exact_match], minlength=len(unique_values))
      frequencies_list.append(counts / n_samples)

    # combine the individual frequency vectors into the stacked global vector
    return np.hstack(frequencies_list)