mcpnet.h

#ifndef MCPNET_H
#define MCPNET_H

#include <cstddef>
#include <cstdint>

#include "rust/cxx.h"

// =============================================================================
// Adaptive partitioning mutual information
// (mcp::kernel::AdaptivePartitionRankMIKernel2)
// =============================================================================
//
// Compute mutual information between two rank-transformed vectors. Both
// inputs must have the same length and contain integer ranks in
// `[0, count)`. The result is in bits (log base 2).
float  adaptive_mi_kernel_f32(rust::Slice<const ::std::size_t> x_ranks,
                              rust::Slice<const ::std::size_t> y_ranks);
double adaptive_mi_kernel_f64(rust::Slice<const ::std::size_t> x_ranks,
                              rust::Slice<const ::std::size_t> y_ranks);

// =============================================================================
// Rank transform (splash::kernel::Rank<IT, std::size_t, 0, false>)
// =============================================================================
//
// Descending rank (matches MCPNet `mi.cpp`), first rank = 0, output type
// `size_t`. `out_vec` must be at least `in_vec.size()` long.
void rank_f32(rust::Slice<const float>  in_vec,
              rust::Slice<::std::size_t> out_vec);
void rank_f64(rust::Slice<const double> in_vec,
              rust::Slice<::std::size_t> out_vec);

// =============================================================================
// Cumulative histogram (mcp::kernel::IntegralCumulativeHistogramKernel)
// =============================================================================
//
// Builds a cumulative histogram over the integer ranks in `in_vec`. `out_vec`
// must be `in_vec.size() + 1` long and zero-initialized on entry.
void integral_cumulative_histogram(rust::Slice<const ::std::size_t> in_vec,
                                   rust::Slice<::std::size_t>       out_vec);

// =============================================================================
// BSpline weights (mcp::kernel::BSplineWeightsKernel)
// =============================================================================
//
// Computes the BSpline weight matrix for a single sample vector. The
// output vector must be at least `num_bins * num_samples + 1` elements long;
// the last element holds the 1D BSpline entropy of the weights.
void bspline_weights_f32(rust::Slice<const float> in_vec,
                         ::std::int32_t num_bins,
                         ::std::int32_t spline_order,
                         ::std::int32_t num_samples,
                         rust::Slice<float> out_vec);
void bspline_weights_f64(rust::Slice<const double> in_vec,
                         ::std::int32_t num_bins,
                         ::std::int32_t spline_order,
                         ::std::int32_t num_samples,
                         rust::Slice<double> out_vec);

// =============================================================================
// BSpline mutual information (mcp::kernel::BSplineMIKernel)
// =============================================================================
//
// Computes mutual information between two variables from their BSpline
// weight vectors as `H(X) + H(Y) - H(X, Y)`. Each input vector must be
// exactly `num_bins * num_samples + 1` elements long (the layout produced
// by `bspline_weights_*`).
float  bspline_mi_kernel_f32(rust::Slice<const float>  first,
                             rust::Slice<const float>  second,
                             ::std::int32_t num_bins,
                             ::std::int32_t num_samples);
double bspline_mi_kernel_f64(rust::Slice<const double> first,
                             rust::Slice<const double> second,
                             ::std::int32_t num_bins,
                             ::std::int32_t num_samples);

// =============================================================================
// Z-score / standard score (splash::kernel::StandardScore)
// =============================================================================
//
// Population standard score (SampleStats = false), matching MCPNet's
// pearson.cpp pipeline so that pairing with `pearson_kernel_*` yields the
// usual Pearson correlation coefficient. `out_vec` must be at least
// `in_vec.size()` long.
void standard_score_f32(rust::Slice<const float>  in_vec,
                        rust::Slice<float>        out_vec);
void standard_score_f64(rust::Slice<const double> in_vec,
                        rust::Slice<double>       out_vec);

// =============================================================================
// Pearson correlation (mcp::correlation::PearsonKernel)
// =============================================================================
//
// Computes the Pearson correlation coefficient. The two inputs are expected
// to already be standard-scored (see `standard_score_*`). Both inputs must
// have the same length.
float  pearson_kernel_f32(rust::Slice<const float>  first,
                          rust::Slice<const float>  second);
double pearson_kernel_f64(rust::Slice<const double> first,
                          rust::Slice<const double> second);

// =============================================================================
// Element-wise combine kernels (mcp::kernel from mcp/transform/combine.hpp)
// =============================================================================
//
// All vector kernels write `count = min(in_a.size(), in_b.size())` output
// elements (or `in.size()` for the unary `scale_kernel`). The output buffer
// must be at least that large.

void add_kernel_f32     (rust::Slice<const float>  in_a,
                         rust::Slice<const float>  in_b,
                         rust::Slice<float>        out);
void add_kernel_f64     (rust::Slice<const double> in_a,
                         rust::Slice<const double> in_b,
                         rust::Slice<double>       out);

void sub_kernel_f32     (rust::Slice<const float>  in_a,
                         rust::Slice<const float>  in_b,
                         rust::Slice<float>        out);
void sub_kernel_f64     (rust::Slice<const double> in_a,
                         rust::Slice<const double> in_b,
                         rust::Slice<double>       out);

void multiply_kernel_f32(rust::Slice<const float>  in_a,
                         rust::Slice<const float>  in_b,
                         rust::Slice<float>        out);
void multiply_kernel_f64(rust::Slice<const double> in_a,
                         rust::Slice<const double> in_b,
                         rust::Slice<double>       out);

void ratio_kernel_f32   (rust::Slice<const float>  in_a,
                         rust::Slice<const float>  in_b,
                         rust::Slice<float>        out);
void ratio_kernel_f64   (rust::Slice<const double> in_a,
                         rust::Slice<const double> in_b,
                         rust::Slice<double>       out);

void max_kernel_f32     (rust::Slice<const float>  in_a,
                         rust::Slice<const float>  in_b,
                         rust::Slice<float>        out);
void max_kernel_f64     (rust::Slice<const double> in_a,
                         rust::Slice<const double> in_b,
                         rust::Slice<double>       out);

void min_kernel_f32     (rust::Slice<const float>  in_a,
                         rust::Slice<const float>  in_b,
                         rust::Slice<float>        out);
void min_kernel_f64     (rust::Slice<const double> in_a,
                         rust::Slice<const double> in_b,
                         rust::Slice<double>       out);

// out = a + coeff * b
void madd_kernel_f32    (rust::Slice<const float>  in_a,
                         rust::Slice<const float>  in_b,
                         double                    coeff,
                         rust::Slice<float>        out);
void madd_kernel_f64    (rust::Slice<const double> in_a,
                         rust::Slice<const double> in_b,
                         double                    coeff,
                         rust::Slice<double>       out);

// out = alpha * a + beta * b
void weighted_mean_kernel_f32(rust::Slice<const float>  in_a,
                              rust::Slice<const float>  in_b,
                              double                    alpha,
                              double                    beta,
                              rust::Slice<float>        out);
void weighted_mean_kernel_f64(rust::Slice<const double> in_a,
                              rust::Slice<const double> in_b,
                              double                    alpha,
                              double                    beta,
                              rust::Slice<double>       out);

// out = coeff * in
void scale_kernel_f32   (rust::Slice<const float>  in,
                         double                    coeff,
                         rust::Slice<float>        out);
void scale_kernel_f64   (rust::Slice<const double> in,
                         double                    coeff,
                         rust::Slice<double>       out);

// =============================================================================
// Stouffer kernels (mcp/transform/stouffer.hpp)
// =============================================================================
//
// The auxiliary (mean, stdev) pairs used in MCPNet are passed here as two
// scalars or two parallel arrays.

float  stouffer_kernel_f32(float  value,
                           float  mean1, float  std1,
                           float  mean2, float  std2);
double stouffer_kernel_f64(double value,
                           double mean1, double std1,
                           double mean2, double std2);

void stouffer_vector_kernel_f32(rust::Slice<const float>  in_vec,
                                float                     aux1_mean,
                                float                     aux1_std,
                                rust::Slice<const float>  aux2_means,
                                rust::Slice<const float>  aux2_stds,
                                rust::Slice<float>        out);
void stouffer_vector_kernel_f64(rust::Slice<const double> in_vec,
                                double                    aux1_mean,
                                double                    aux1_std,
                                rust::Slice<const double> aux2_means,
                                rust::Slice<const double> aux2_stds,
                                rust::Slice<double>       out);

// out[k] = (zi[k] + zj[k]) / sqrt(2)
void zscored_stouffer_kernel_f32(rust::Slice<const float>  zi,
                                 rust::Slice<const float>  zj,
                                 rust::Slice<float>        out);
void zscored_stouffer_kernel_f64(rust::Slice<const double> zi,
                                 rust::Slice<const double> zj,
                                 rust::Slice<double>       out);

float  weighted_stouffer_kernel_f32(float  value,
                                    float  mean1, float  std1,
                                    float  mean2, float  std2);
double weighted_stouffer_kernel_f64(double value,
                                    double mean1, double std1,
                                    double mean2, double std2);

// =============================================================================
// Context Likelihood of Relatedness (mcp/transform/clr.hpp)
// =============================================================================

float  clr_kernel_f32(float  value,
                      float  mean1, float  std1,
                      float  mean2, float  std2);
double clr_kernel_f64(double value,
                      double mean1, double std1,
                      double mean2, double std2);

void clr_vector_kernel_f32(rust::Slice<const float>  in_vec,
                           float                     aux1_mean,
                           float                     aux1_std,
                           rust::Slice<const float>  aux2_means,
                           rust::Slice<const float>  aux2_stds,
                           rust::Slice<float>        out);
void clr_vector_kernel_f64(rust::Slice<const double> in_vec,
                           double                    aux1_mean,
                           double                    aux1_std,
                           rust::Slice<const double> aux2_means,
                           rust::Slice<const double> aux2_stds,
                           rust::Slice<double>       out);

// out[k] = sqrt(zi[k]^2 + zj[k]^2)
void zscored_clr_kernel_f32(rust::Slice<const float>  zi,
                            rust::Slice<const float>  zj,
                            rust::Slice<float>        out);
void zscored_clr_kernel_f64(rust::Slice<const double> zi,
                            rust::Slice<const double> zj,
                            rust::Slice<double>       out);

// =============================================================================
// Filter kernels (mcp/filter/threshold.hpp, mcp/filter/dpi.hpp, mcp/filter/mcp.hpp)
// =============================================================================

// ---- threshold.hpp ----------------------------------------------------------
// out[i] = (min <= in[i] <= max) ? in[i] : default_val
void threshold_kernel_f32(rust::Slice<const float>  in_vec,
                          float min_thresh, float max_thresh,
                          float default_val,
                          rust::Slice<float>        out);
void threshold_kernel_f64(rust::Slice<const double> in_vec,
                          double min_thresh, double max_thresh,
                          double default_val,
                          rust::Slice<double>       out);

// out[i] = (in[i] <= min OR in[i] >= max) ? in[i] : default_val
void inverted_threshold_kernel_f32(rust::Slice<const float>  in_vec,
                                   float min_thresh, float max_thresh,
                                   float default_val,
                                   rust::Slice<float>        out);
void inverted_threshold_kernel_f64(rust::Slice<const double> in_vec,
                                   double min_thresh, double max_thresh,
                                   double default_val,
                                   rust::Slice<double>       out);

// out[i] = (min <= aux[i] <= max) ? in[i] : default_val
void threshold2_kernel_f32(rust::Slice<const float>  in_vec,
                           rust::Slice<const float>  aux,
                           float min_thresh, float max_thresh,
                           float default_val,
                           rust::Slice<float>        out);
void threshold2_kernel_f64(rust::Slice<const double> in_vec,
                           rust::Slice<const double> aux,
                           double min_thresh, double max_thresh,
                           double default_val,
                           rust::Slice<double>       out);

// out[i] = (aux[i] <= min OR aux[i] >= max) ? in[i] : default_val
void inverted_threshold2_kernel_f32(rust::Slice<const float>  in_vec,
                                    rust::Slice<const float>  aux,
                                    float min_thresh, float max_thresh,
                                    float default_val,
                                    rust::Slice<float>        out);
void inverted_threshold2_kernel_f64(rust::Slice<const double> in_vec,
                                    rust::Slice<const double> aux,
                                    double min_thresh, double max_thresh,
                                    double default_val,
                                    rust::Slice<double>       out);

// out[i] = ((negate ^ (mask[i] != 0)) ? in[i] : default_val
// `mask` is a u8 slice where 0 means false and any non-zero means true.
void mask_kernel_f32(rust::Slice<const float>     in_vec,
                     rust::Slice<const ::std::uint8_t> mask,
                     bool negate, float default_val,
                     rust::Slice<float>           out);
void mask_kernel_f64(rust::Slice<const double>    in_vec,
                     rust::Slice<const ::std::uint8_t> mask,
                     bool negate, double default_val,
                     rust::Slice<double>          out);

// ---- dpi.hpp ----------------------------------------------------------------
// Returns true when the edge (x, z) should be deleted under DPI: there is some
// y for which mi_xy and mi_yz both meet the i_xz tolerance bound.
// EQUALITY_EXCLUSION is set to true (the default in the C++ code).
bool dpi_kernel_f32(::std::size_t x, ::std::size_t z,
                    rust::Slice<const float>  row_x,
                    rust::Slice<const float>  row_z,
                    double tolerance);
bool dpi_kernel_f64(::std::size_t x, ::std::size_t z,
                    rust::Slice<const double> row_x,
                    rust::Slice<const double> row_z,
                    double tolerance);

// As above but with a per-node TF mask (TFs[i] > 0 marks node `i` as a TF).
bool dpi_tf_kernel_f32(::std::size_t x, ::std::size_t z,
                       rust::Slice<const float>  row_x,
                       rust::Slice<const float>  row_z,
                       rust::Slice<const float>  tfs,
                       double tolerance);
bool dpi_tf_kernel_f64(::std::size_t x, ::std::size_t z,
                       rust::Slice<const double> row_x,
                       rust::Slice<const double> row_z,
                       rust::Slice<const double> tfs,
                       double tolerance);

// ---- mcp.hpp ----------------------------------------------------------------
// max_y(min(row_x[y], row_z[y])) excluding y in {x, z}; if `tfs` is non-empty,
// it is used as a per-node mask (only y with tfs[y] > 0 are considered).
float  mcp2_maxmin_kernel_f32(::std::size_t x, ::std::size_t z,
                              rust::Slice<const float>  row_x,
                              rust::Slice<const float>  row_z,
                              rust::Slice<const float>  tfs);
double mcp2_maxmin_kernel_f64(::std::size_t x, ::std::size_t z,
                              rust::Slice<const double> row_x,
                              rust::Slice<const double> row_z,
                              rust::Slice<const double> tfs);

// row_x[z] / max_y(min(row_x[y], row_z[y]))
float  mcp_ratio_kernel_f32(::std::size_t x, ::std::size_t z,
                            rust::Slice<const float>  row_x,
                            rust::Slice<const float>  row_z,
                            rust::Slice<const float>  tfs);
double mcp_ratio_kernel_f64(::std::size_t x, ::std::size_t z,
                            rust::Slice<const double> row_x,
                            rust::Slice<const double> row_z,
                            rust::Slice<const double> tfs);

// 1 - ratio (with optional clamping into [0, 1]).
float  mcp_tolerance_kernel_f32(::std::size_t x, ::std::size_t z,
                                rust::Slice<const float>  row_x,
                                rust::Slice<const float>  row_z,
                                rust::Slice<const float>  tfs,
                                bool clamped);
double mcp_tolerance_kernel_f64(::std::size_t x, ::std::size_t z,
                                rust::Slice<const double> row_x,
                                rust::Slice<const double> row_z,
                                rust::Slice<const double> tfs,
                                bool clamped);

// row_x[z] - max_y(min(row_x[y], row_z[y]))
float  mcp_diff_kernel_f32(::std::size_t x, ::std::size_t z,
                           rust::Slice<const float>  row_x,
                           rust::Slice<const float>  row_z,
                           rust::Slice<const float>  tfs);
double mcp_diff_kernel_f64(::std::size_t x, ::std::size_t z,
                           rust::Slice<const double> row_x,
                           rust::Slice<const double> row_z,
                           rust::Slice<const double> tfs);

#endif // MCPNET_H