#include "fracture.h"

cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp,
                              bool symmetric, cholmod_common *c) {
  const graph_t *g = net->graph;
  uint_t nv;
  uint_t ne;
  uint_t nre;
  uint_t *ev;

  if (use_gp) {
    nv = net->dim;
    ne = net->nep;
    nre = (int_t)net->nep - (int_t)net->num_broken;
    ev = net->evp;
  } else if (dual) {
    nv = g->dnv;
    ne = g->ne;
    nre = net->num_broken;
    ev = g->dev;
  } else {
    nv = g->nv;
    ne = g->ne;
    nre = (int_t)g->ne - (int_t)net->num_broken;
    ev = g->ev;
  }

  uint_t nnz = nre;

  cholmod_triplet *t =
      CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);

  int_t *ri = (int_t *)t->i;
  int_t *ci = (int_t *)t->j;
  double *ai = (double *)t->x;

  t->nnz = nnz;

  uint_t a = 0;

  for (uint_t i = 0; i < ne; i++) {
    if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
      uint_t v1 = ev[2 * i];
      uint_t v2 = ev[2 * i + 1];

      uint_t s1 = v1 < v2 ? v1 : v2;
      uint_t s2 = v1 < v2 ? v2 : v1;

      ri[a] = s2;
      ci[a] = s1;
      ai[a] = 1;

      a++;
    }
  }

  cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
  CHOL_F(free_triplet)(&t, c);

  if (!symmetric) {
    cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
    CHOL_F(free_sparse)(&s, c);
    s = tmp_s;
  }

  return s;
}

cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) {
  const graph_t *g = net->graph;
  uint_t nv = net->dim;
  uint_t ne = net->nep;
  uint_t *ev = net->evp;

  uint_t nnz = nv;

  cholmod_triplet *temp_m =
      CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);

  int_t *rowind = (int_t *)temp_m->i;
  int_t *colind = (int_t *)temp_m->j;
  double *acoo = (double *)temp_m->x;

  temp_m->nnz = nnz;

  for (uint_t i = 0; i < nv; i++) {
    rowind[i] = i;
    colind[i] = i;
    acoo[i] = 0;
  }

  cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);

  for (uint_t i = 0; i < ne; i++) {
    if (!net->fuses[i]) {
      uint_t v1 = ev[2 * i];
      uint_t v2 = ev[2 * i + 1];

      acoo[v1]++;
      acoo[v2]++;
    }
  }

  if (net->voltage_bound && g->boundary != TORUS_BOUND) {
    for (uint_t i = 0; i < net->dim; i++) {
      uint_t v = g->nbi[i];
      for (uint_t j = 0; j < g->vei[v + 1] - g->vei[v]; j++) {
        uint_t e = g->ve[g->vei[v] + j];
        uint_t v0 = g->ev[2 * e];
        uint_t v1 = g->ev[2 * e + 1];

        if (g->bq[v0] || g->bq[v1]) {
          acoo[i]++;
        }
      }
    }
  } else {
    acoo[0]++;
  }

  for (uint_t i = 0; i < nv; i++) {
    if (acoo[i] == 0)
      acoo[i]++;
  }

  // assert(CHOL_F(check_triplet)(temp_m, c));

  cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
  // assert(CHOL_F(check_sparse)(t_out, c));

  double alpha[2] = {1, 0};
  double beta[2] = {-1, 0};
  cholmod_sparse *laplacian =
      CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);

  CHOL_F(free_sparse)(&t_out, c);
  CHOL_F(free_sparse)(&adjacency, c);
  CHOL_F(free_triplet)(&temp_m, c);

  return laplacian;
}