#include "fracture.h"

double *net_currents(const net_t *net, const double *voltages,
                     cholmod_common *c) {
  uint_t ne = net->graph->ne;
  uint_t dim = net->graph->nv;
  cholmod_sparse *voltcurmat = net->voltcurmat;

  cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);

  double *tmp_x = x->x;
  x->x = (void *)voltages;

  cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c);

  double alpha[2] = {1, 0};
  double beta[2] = {0, 0};
  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);

  double *currents = (double *)malloc(ne * sizeof(double));

  for (int i = 0; i < ne; i++) {
    if (net->fuses[i]) {
      currents[i] = 0;
    } else {
      currents[i] = ((double *)y->x)[i];
    }
  }

  if (net->graph->boundary == TORUS_BOUND) {
    for (uint_t i = 0; i < net->graph->bi[1]; i++) {
      uint_t e = net->graph->b[i];
      uint_t v1 = net->graph->ev[2 * e];
      uint_t v2 = net->graph->ev[2 * e + 1];
      double v1y = net->graph->vx[2 * v1 + 1];
      double v2y = net->graph->vx[2 * v2 + 1];

      if (v1y > v2y) {
        currents[e] += 1;
      } else {
        currents[e] -= 1;
      }
    }
  }

  x->x = tmp_x;
  CHOL_F(free_dense)(&x, c);
  CHOL_F(free_dense)(&y, c);

  return currents;
}