/* sparse.gp */

/*
zCs = [C,E]~
- C: a t_VECSMALL of indices
- E: a t_VECSMALL of coefficients
representing sum_i E[i]*b_{C[i]} (b_j standard basis).

SMS format: https://hpac.imag.fr/

branch aurel-amt or pari-A2026

A KzMs (complex of sparse matrices) is a t_VEC with two components:
- a t_VECSMALL of dimensions d_0,...,d_n
- a t_VEC of zMs D_1,...,D_n representing differentials
  D_i : Z^{d_{i-1}} <- Z^{d_i}
  s.t. D_i*D_{i+1} = 0 for all 1<=i<n.
*/

install("zMs_to_ZM","GL","smattomat");
install("KzMs_simplify","G");

tosmat(M,ncols=-1) =
{
  if(type(M)=="t_MAT",
    my(m);
    m = vector(#M);
    for(j=1,#m,
      my(col=M[,j],C=List(),E=List());
      for(i=1,#col,
        if(col[i],
          listput(~C,i);
          listput(~E,col[i]);
        );
      );
      m[j] = [Vecsmall(C),Vecsmall(E)]
    );
    m
  ,type(M)=="t_VEC",
    M = vecsort(M,[2,1]); \\sort by column then row
    my(j=1,C=List(),E=List(),m=List());
    for(k=1,#M,
      my([i,jj,c] = M[k]);
      if(j!=jj,
        listput(~m,[Vecsmall(C),Vecsmall(E)]);
        j++;
        while(j<jj,
          listput(~m,[Vecsmall([]),Vecsmall([])]);
          j++;
        );
        C = List();
        E = List();
      );
      listput(~C,i);
      listput(~E,c);
    );
    listput(~m,[Vecsmall(C),Vecsmall(E)]);
    if(ncols>=0,
      while(j<ncols,
        listput(~m,[Vecsmall([]),Vecsmall([])]);
        j++;
      );
    );
    Vec(m)
  ,type(M)=="t_LIST", \\map
    M = apply(concat,Col(Mat(M)))~;
    tosmat(M,ncols)
  );
};

smattovec(M) =
{
  my(L = List());
  for(j=1,#M,
    my([C,E] = M[j]);
    for(k=1,#E,
      my(i,c); 
      i = C[k];
      c = E[k];
      listput(~L,[i,j,c]);
    );
  );
  Vec(L)
};

smattomap(M) =
{
  my(m = Map());
  for(j=1,#M,
    my([C,E] = M[j]);
    for(k=1,#E,
      my(i,c); 
      i = C[k];
      c = E[k];
      mapput(~m, [i,j], c);
    );
  );
  m
};

smattranspose(M) =
{
  my(L = List());
  for(j=1,#M,
    my([C,E] = M[j]);
    for(k=1,#E,
      my(i,c); 
      i = C[k];
      c = E[k];
      listput(~L,[j,i,c]);
    );
  );
  tosmat(Vec(L))
};

/* import from SMS format */
smatimport(file) =
{
  my(v = readstr(file),v2,nrows,ncols,data,m=0);
  v2 = vector(#v-1);
  data = strsplit(v[1]," ");
  if(data[3]!="M", error("wrong format")); \\TODO better error message
  nrows = eval(data[1]);
  ncols = eval(data[2]);
  for(i=2,#v,
    data = strsplit(v[i]," ");
    my(c = apply(eval,data));
    if(c[1]==0, m=i; break);
    v2[i-1] = c;
  );
  if(!m, error("wrong file format"));
  v2 = v2[1..m-2];
  [nrows,ncols,tosmat(v2,ncols)]
};

/* export to SMS format */
smatexport(file,M,nrows) =
{
  my(ncols = #M, V);
  write(file, nrows, " ", ncols, " M");
  V = smattovec(M);
  V = vecsort(V, [1,2]);
  foreach(V,v,
    write(file, v[1], " ", v[2], " ", v[3]);
  );
  write(file, "0 0 0");
};

densecomplex(C) =
{
  my([dims,diffs] = C);
  [dims, [smattomat(diffs[i],dims[i]) | i <- [1..#diffs]]];
};

/*
  No need to compute integer kernel:
  tor(Z/B) = tor(C/B) because Z is saturated (it is a kernel)
  Z/B -> C/B is injective, and surjective on torsion:
  n*c = 0 mod B => n*c is in B => n*c in Z => c in Z.
*/
densehomology(C,a,b,verb=0) =
{
  my(dims,diffs,t,i,rprec,m,r,tor,betti,snf,H);
  H = List();
  [dims,diffs] = C;
  if(verb,printf("dim %d: %d cells\n", a-1, dims[1]));
  rprec = matrank(diffs[1]);
  for(k=a,b-1,
    i = 1+k-a;
    if(verb,printf("\ndim %d: %d cells\n", k, dims[i+1]));
    if(verb,printf("homology H_%d: ",k));
    t = getabstime();
    if(diffs[i+1]==0,
      snf=vector(dims[i+1])
    ,\\else
      snf = matsnf(diffs[i+1],4);
    );
    t = getabstime()-t;
    tor = [d | d <- snf, d!=0];
    m = #snf - #tor;
    r = dims[i+1] - m;
    betti = dims[i+1] - (r+rprec);
    if(verb,
      if(betti==0 && #tor==0, print1("0"));
      if(betti==1, print1("Z"), betti>1, printf("Z^%d", betti));
      if(betti>0 && #tor>0, print1(" + "));
      if(#tor>0, print1(tor));
      print1(" ; #tor=", vecprod(tor));
      print("\n[",strtime(t),"]");
    );
    listput(~H, [betti,tor]);
    rprec = r;
  );
  if(verb,printf("\ndim %d: %d cells\n", b, dims[b-a+2]));
  Vec(H);
};

zCs_permute(C,prow) =
{
  my(c,e,psort);
  [c,e] = C;
  c = vector(#c,i,prow[c[i]]);
  psort = vecsort(c,,1);
  [Vecsmall(vector(#c,i,c[psort[i]])), Vecsmall(vector(#e,i,e[psort[i]]))]~
};

zMs_permute(M,prow,pcol) =
{
  vector(#M,i,zCs_permute(M[pcol[i]],prow));
};

permutecomplex(C,Lperm) =
{
  my([dims,diffs]=C,prow,pcol);
  for(i=1,#diffs,
    if(dims[i]>0,
      prow = Lperm[i]^(-1);
      pcol = Lperm[i+1];
      diffs[i] = zMs_permute(diffs[i],prow,pcol);
    );
  );
  [dims,diffs]
};

randperm(n) =
{
  my(v = [1..n],j,tmp);
  for(i=1,n,
    j = random([i,n]);
    tmp = v[i];
    v[i] = v[j];
    v[j] = tmp;
  );
  Vecsmall(v)
};

zCs_nbcoefs(C) = #C[1];
zCs_nbinvcoefs(C) = #[c | c <- C[2], abs(c)==1];
zMs_nbcoefs(M) =
{
  sum(i=1,#M,zCs_nbcoefs(M[i]));
};
zMs_nbinvcoefs(M) =
{
  sum(i=1,#M,zCs_nbinvcoefs(M[i]));
};
K_nbcoefs(C) =
{
  sum(i=1,#C[2],zMs_nbcoefs(C[2][i]));
};
K_nbinvcoefs(C) =
{
  sum(i=1,#C[2],zMs_nbinvcoefs(C[2][i]));
};

eps = 0.01;
homology(C,a,b,verb=0) =
{
  my(C2=0,dC,Lperm,count=0,t,H);
  if(verb,
    print1("starting dimensions: ", C[1]);
    print1(" nb coefs: ", K_nbcoefs(C));
    print(" nb inv coefs: ", K_nbinvcoefs(C));
  );
  t = getabstime();
  while(C2==0 || vecsum(Vec(C[1]))<=(1-eps)*vecsum(Vec(C2[1])), 
    C2 = C;
    C = KzMs_simplify(C2);
    if(verb,
      print1("dimensions: ", C[1]);
      print1(" nb coefs: ", K_nbcoefs(C));
      print(" nb inv coefs: ", K_nbinvcoefs(C));
    );
  );
  C2=0;
  if(verb,print("now shuffling."));
  while(count<1, 
    C2 = C;
    Lperm = vector(#C2[1],i,randperm(C2[1][i]));
    C2 = permutecomplex(C2,Lperm);
    C = KzMs_simplify(C2);
    if(verb,
      print1("dimensions: ", C[1]);
      print1(" nb coefs: ", K_nbcoefs(C));
      print(" nb inv coefs: ", K_nbinvcoefs(C));
    );
    if(C2!=0 && vecsum(Vec(C[1]))>(1-eps)*vecsum(Vec(C2[1])), count++, count=0);
  );
  t = getabstime()-t;
  if(verb,
    print("time simplification: ", strtime(t));
    print("dimensions: ", C[1]);
    print("nb coefs: ", vector(#C[2],i,zMs_nbcoefs(C[2][i])));
    print("nb inv coefs: ", vector(#C[2],i,zMs_nbinvcoefs(C[2][i])));
    print1("nb coefs: ", K_nbcoefs(C));
    print(" nb inv coefs: ", K_nbinvcoefs(C));
  );
  dC = densecomplex(C);
  H = densehomology(dC,a,b);
  if(verb,print("done."));
  H
};