/* Generate H for  LDPC*/

/* Based on sparse.c by Matt Davey <mcdavey@mrao.cam.ac.uk> with his permission*/

#include <math.h>

#include "mex.h"



/* Input Arguments: parameters*/

#define	M_IN	prhs[0]          /* number of parity checks */

#define	N_IN	prhs[1]          /* blocklength */ 

#define	T_IN	prhs[2]          /* mean column weight */

#define	Q_IN	prhs[3]          /* GF base  */

#define	SEED_IN	prhs[4]          /* seed for random generator */





/* Output Arguments: matrices*/

#define	H_OUT	plhs[0]





void mexFunction(

                 int nlhs,       mxArray *plhs[],

                 int nrhs, const mxArray *prhs[]

		 )

{

  short	**M_list, **N_list, *M_target;

  double *pp,*sr,*s,*ss;

  int N,M,q,i,j,k,nzmax,*irs,*jcs,l,

	  tr,tm,tc,done,redo,tmp,regime,tr_max,t_max, m_low;

  float t;

  long seed;

  void adjust(int *, int *, int *, int *);

  unsigned int K2,M2;

  char c;

  mxArray *arg_in[2], *arg_out[1]; /* to call rand generator of Matlab*/



  /* Check for proper number of arguments */

  if (nrhs != 5) {

    mexErrMsgTxt("GENERATE requires five input arguments.");

  } else if (nlhs > 1) {

    mexErrMsgTxt("GENERATE requires one output argument.");

  }



  pp = mxGetPr(N_IN);    N = (int) (*pp);

  pp = mxGetPr(M_IN);    M = (int) (*pp);

  pp = mxGetPr(T_IN);    t = (float) (*pp);

  pp = mxGetPr(Q_IN);    q = (int) (*pp);

  pp = mxGetPr(SEED_IN); seed = (int) (*pp);

  



  arg_in[0] = mxCreateString("state");

  arg_in[1] = mxCreateDoubleMatrix(1, 1, mxREAL);

  s = mxGetPr(arg_in[1]);

  s[0] = seed; /* this will be used to call rand*/



  /* initialize random generator */

  mexCallMATLAB(0, NULL, 2, arg_in, "rand"); /* rand('state',seed) */

  s[0] = 1; /* use s to store "1"*/





/* I have no idea about the details of the following - please

ask the author - :).  igor*/











/* Generate some sparse matrices for error-correcting codes.  Supply

 * the following parameters:

 * 	N: Blocklength

 * 	M: Number of parity checks

 * 	t: Mean column weight.

 *

 * If t<3, then we generate weight 2 columns systematically in the

 * form of blocks of identity matrices to reduce probability of

 * getting short cycle lengths.  Can't have more than M weight 2

 * columns, though.

 *

 * Having generated any weight 2 columns, we fill the rest.  Calculate

 * number of columns to fill with floor(t) and number with ceiling(t).

 * Find mean row weight and calculate number of rows to fill with

 * floor(r) and ceiling(r).

 *

 * We fill rows as follows:

 * (Regime 0): Fill so that rows contain <= tr ones until M-tm rows

 * 	    contain tr ones.

 * (Regime 1): Fill so that rows contain <= (tr+1) ones until tm rows

 * 	    contain (tr+1) ones.

 * (Regime 2): Fill remaining rows containing < tr ones until done.

 */



 

  /* N=6112;M=4512;t=2.3;*/



     

  t_max=(int)ceil(t);

  M_list=(short **)mxMalloc(N*sizeof(short *));

  M_target=(short *)mxMalloc(M*sizeof(short *));

  N_list=(short **)mxMalloc(M*sizeof(short *));

  for(i=0;i<N;i++){

    M_list[i]=(short *)mxMalloc((t_max+1)*sizeof(short));

  }

  i=0;

  /* Do we have any weight 2 columns?

   *

   * If so, do this first.  Remember that M mightn't have a large

   * power of two as a divisor, so might need to find some M'<M to use

   * as unit length.

   */

  K2=0;

  if(t<3){

    K2=ceil((double)N*(3-t));

	if(K2>M)

		mexErrMsgTxt("GENERATE: Can't have more than M weight 2 columns.");

    j=2;

    done=0;

    for(i=0;!done;i++){

      M2=floor((double)M/(double)j);

      if((M2*(j-1))>=K2) done=1;

      j*=2;

    }

    M2*=(j/4);

  }

  /* 

   * i contains number of identity blocks we'll need.... */

  tr=((short)floor((double)(t*N)/(double)M));

  /* Now we want `tr' to be final minimum row weight, `tr_max' to be

   * final maximum row weight, `tm' to be number of rows which will

   * have weight greater than `tr'.  'tc' will be a running count of

   * how many rows we still have to fill up to weight `tr'.  Once we

   * hit this many, we can start overfilling rows.

   */

  if (i>tr){

    tr_max=i;

    /* If identity blocks make overheavy rows, we need to calculate

     * the minimum row weight.

     */

    done=0;

    k=1;

    j=floor((double)t*N)-2*K2; /* Number of ones left to distribute */

    for(i=0;!done;i++){

      /* (M-2*M2) rows will be empty after identity blocks */

      j-=((M-2*M2)+(2*M2*(k-1))/k);

      if(j<0) {

	done=1;

      }

      else {

	k*=2;

      }

    }

    tr=i-1;

    tm=M+j;

  }

  else {

    /* This is easier! */

    tr_max=tr+1;

    tm=(int)floor((((double)t*N)/(double)M-tr)*M +0.5);

  }

  tc=M-tm;

  for(i=0;i<M;i++){

    N_list[i]=(short *)mxMalloc((tr_max+1)*sizeof(short));  

  }

  for(i=0;i<M;i++){

    N_list[i][0]=0;

  }



  regime=0;

  /* Generate weight 2 columns.  First create two identity matrices on

   * top of each other, then two 1/2 size ones in the lower rows of

   * the matrix, and so on.

   *

   * j:   length of current identity block

   * k:   base of this block

   * i:   current column position

   */

  j=M2;

  k=0;

  i=0;

  while(i<K2){

    for(;(i-k)<j && i<K2;i++){

      M_list[i][0]=2;

      M_list[i][1]=i;

      M_list[i][2]=i+j;

      N_list[i][0]++;

      if(N_list[i][0]==tr) adjust(&tm,&tr,&tc,&regime);

      N_list[i][N_list[i][0]]=i;

      N_list[i+j][0]++;

      if(N_list[i+j][0]==tr) adjust(&tm,&tr,&tc,&regime);

      N_list[i+j][N_list[i+j][0]]=i;

    }

    k=i;

    j/=2;

  }

 

 /* Now fill the unsystematic columns, ensuring weight per row as even as poss. */

  i=K2;

  if(K2==0){

    /* Fill low weight columns */

    for(i=0;i<(int)(N*(t_max-t)+0.5);i++){

      for(k=1;k<=(int)floor(t);k++){

	done=0;

	do {

      mexCallMATLAB(1, arg_out,1 , &arg_in[1], "rand"); /* ss = rand(1) */

	  ss = mxGetPr(arg_out[0]);

	  j=(short)floor(M*ss[0]);

      mexCallMATLAB(1, arg_out,1 , &arg_in[1], "rand");

	  ss = mxGetPr(arg_out[0]);

	  if((ss[0])<(1-(double)N_list[j][0]/(double)tr)) {

	    done=1;

	    for(l=1;l<k;l++) if(j==M_list[i][l]) done=0;

	  }

	} while(!done);

	N_list[j][0]++;

	N_list[j][N_list[j][0]]=i;

	if(N_list[j][0]==tr) adjust(&tm,&tr,&tc,&regime);

	M_list[i][k]=j;

      }

      M_list[i][0]=k-1;

    }

  }

  redo=1;

  for(;i<N;i++){

    fprintf(stderr,"%d\r",i);

    for(k=1;k<=t_max;k++){

      done=0;

      do {

	/* find the lowest weight rows, and fill one of them */

	if(redo){

	  l=tr_max;

	  for(j=0;j<M;j++) if (N_list[j][0]<l) l=N_list[j][0];

	  m_low=0;

	  for(j=0;j<M;j++) if (N_list[j][0]==l) {M_target[m_low]=j; m_low++;}

	}

	mexCallMATLAB(1, arg_out,1 , &arg_in[1], "rand");

	ss = mxGetPr(arg_out[0]);

	j=M_target[tmp=(short)floor(m_low*ss[0])];

	/*	if(ss[0]<(1-(double)N_list[j][0]/(double)tr)) {*/

	  done=1;

	  for(l=1;l<k;l++) if(j==M_list[i][l]) done=0;

	  if(done==1){

	    if(m_low==1) redo=1;

	    else {

	      for(;tmp<(m_low-1);tmp++) M_target[tmp]=M_target[tmp+1];

	      m_low--;

	      redo=0;

	    }

	  }

	  /*	}*/

	  

      } while(!done);

      N_list[j][0]++;

      N_list[j][N_list[j][0]]=i;

      if(N_list[j][0]==tr) adjust(&tm,&tr,&tc,&regime);

      M_list[i][k]=j;

    }

    M_list[i][0]=k-1;

  }

  tr=((short)ceil((double)(3*N-K2)/(double)M));

  



  for(i=0;i<M;i++){

    mxFree(N_list[i]);

  }

  mxFree(N_list);





/* done generating H matrix - positions only */





    /* Allocate space for sparse matrix */

    nzmax=0; for(j=0 ; j<N ; j++) nzmax +=M_list[j][0];

	mexPrintf("%d \n",nzmax);

    /* NOTE: The maximum number of non-zero elements cannot be less

       than the number of columns in the matrix. */

    if (N>nzmax){

	nzmax=N;

    }

    plhs[0] = mxCreateSparse(M,N,nzmax,mxREAL);

    sr  = mxGetPr(plhs[0]);

    irs = mxGetIr(plhs[0]);  /* row */

    jcs = mxGetJc(plhs[0]);  /* column */

    

    /* Copy nonzeros */

	/* elements of GFq are substitudet randomply */

    k = 0; 

    for (j=0; (j<N ); j++) {

	jcs[j] = k;

	for (i=1; (i<=M_list[j][0] ); i++) {

		mexCallMATLAB(1, arg_out,1 , &arg_in[1], "rand"); /* ss = rand(q) */

		ss = mxGetPr(arg_out[0]);

		sr[k] = floor(1 + ss[0]*(q-1));

		irs[k] = M_list[j][i];

		k++;	    

	}

    }

    jcs[N] = k;







  for(i=0;i<N;i++){

    mxFree(M_list[i]);

  }

  mxFree(M_list);



  return;

}







void adjust(int *tm, int *tr, int *tc, int *regime){

  switch(*regime){

  case 0:

    (*tc)--;

    if((*tc)==0){

      *regime=1;

      (*tr)++;

    }

    break;

  case 1:

    (*tm)--;

    if((*tm)==0){

      *regime=2;

      (*tr)--;

    }

    break;

  }

}