#include "header_cpp.h"


int char2int(char a);
long timeseed();
double ran3(long *idum);
float ran4(long *idum);
char **shaffle_vertical(char **ali, int nal, int alilen, int N);

int main(int argc, char *argv[]) {

//  cout<<"***"<<endl;

  ifstream hmmaln(argv[1], ios::in);
  if(!hmmaln) {
    cout << "cannot open the hmm alignment file " << endl;
    exit (0 );
  }

  int seqnum = 0;
  int alilen = 0;

  char seqname[100];
  char seq[1000];
  int N;
 
  while(hmmaln>>seqname>>seq) {
    N = strlen(seq); seqnum ++;
    if(seqnum==1) alilen = N;
    if(alilen!=N) {
	cout <<"Not all the sequences are of the same length"<< endl;
	exit(0);
    }
 }

// cout<<"***"<<endl;

 //cout <<seqnum<< "\t"<< alilen<< endl;
 
 int i =0;
 //char sequence[seqnum][alilen];
 char **sequence = new char *[seqnum];
 for(i=0;i<seqnum;i++) sequence[i] = new char [alilen+1];
 char seqnames[seqnum][100];

 hmmaln.clear();
 hmmaln.seekg( 0 );

 i = 0;
 while(hmmaln>>seqnames[i]>>sequence[i]) {
   //cout <<i << endl;
   i++;if(i==seqnum) break;
 }

 //for(i=0;i<seqnum;i++) cout<<sequence[i]<<endl;

 char **shaffle;

 //shaffle = shaffle_vertical(sequence, seqnum, alilen, 2);

 if(argc<=2) {
    shaffle = shaffle_vertical(sequence, seqnum, alilen, 1);
 }
 else {

    int n = atoi(argv[2]);
    shaffle = shaffle_vertical(sequence, seqnum, alilen, n);
 }



}


char **shaffle_vertical(char **ali, int nal, int alilen, int N)
{

/* shaffle a sequence alignment ali[0..nal-1][0..alilen-1] vertically 
   for N times. 
*/

  int i, j, k, l;

  long sseed = -1;

  char **shaffle_ali = new char *[N*nal];
  for(i=0;i<N*nal;i++) shaffle_ali[i] = new char [alilen];

  char **original_ali = new char *[nal];
  for(i=0;i<nal;i++) original_ali[i] = new char [alilen];

  int random_num;

  //cout << "N="<<N<<endl;
  for(i=0;i<N;i++) {
	for(j=0;j<nal;j++) for(k=0;k<alilen;k++) original_ali[j][k]=ali[j][k];

   	/*for(j=0;j<nal;j++) {
	  for(k=0;k<alilen;k++) {
		cout << original_ali[j][k];
	  }
	  cout << endl;
	}*/

   	//cout << "i=" << i<< endl;

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

	  for(k=0;k<nal;k++)  {
	     random_num = (int) (ran4(&sseed)*(nal-k));
   	     //cout << random_num << " " ;
	     shaffle_ali[i*nal+k][j] = original_ali[random_num][j];
	     //cout << i<< " "<< i*nal+k << " " << j<<" " <<  shaffle_ali[i*nal+k][j] <<endl;
	     for(l=random_num+1;l<(nal-k);l++) {
		original_ali[l-1][j] = original_ali[l][j];
	     }
	  }
        }
  }

  for(i=0;i<N*nal;i++) {
	cout << "Ran"<< i << "\t";
	for(j=0;j<alilen;j++) {
	   cout << shaffle_ali[i][j];
	}
	cout << endl;
  }

  return shaffle_ali;

}


int char2int(char a) {
  switch(a) {
    case 'c': case 'C': return 0; break;
    case 'm': case 'M': return 1; break;
    case 'f': case 'F': return 2;break;
    case 'i': case 'I': return 3; break;
    case 'l': case 'L': return 4; break;
    case 'v': case 'V': return 5; break;
    case 'w': case 'W': return 6; break;
    case 'y': case 'Y': return 7; break;
    case 'a': case 'A': return 8; break;
    case 'g': case 'G': return 9; break;
    case 't': case 'T': return 10; break;
    case 's': case 'S': return 11; break;
    case 'q': case 'Q': return 12; break;
    case 'n': case 'N': return 13; break;
    case 'e': case 'E': return 14; break;
    case 'd': case 'D': return 15; break;
    case 'h': case 'H': return 16; break;
    case 'r': case 'R': return 17; break;
    case 'k': case 'K': return 18; break;
    case 'p': case 'P': return 19; break;
    case 'x': case 'X': case '.': case '-': return 20; break;
  }
}

long timeseed(){

/* used to generate a negative long int as a seed for ran1 and ran2 from current time */
long seed;
time_t ttime;
time(&ttime);
seed=921303200-ttime;
if(seed>0)seed=-seed;if(seed==0){fprintf(stderr, "Something is wrong with time");exit(0);}
//fprintf(stderr,"The seed is %ld\n",seed);
return seed;
}


#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)

double ran3(long *idum)
/* returns uniform random deviate between 0.0 and 1.0
idum - is a negative integer, do not alter it for successive calls*/

{
        static int inext,inextp;
        static long ma[56];
        static int iff=0;
        long mj,mk;
        int i,ii,k;

        if (*idum < 0 || iff == 0) {
                iff=1;
                mj=MSEED-(*idum < 0 ? -*idum : *idum);
                mj %= MBIG;
                ma[55]=mj;
                mk=1;
                for (i=1;i<=54;i++) {
                        ii=(21*i) % 55;
                        ma[ii]=mk;
                        mk=mj-mk;
                        if (mk < MZ) mk += MBIG;
                        mj=ma[ii];
                }
                for (k=1;k<=4;k++)
                        for (i=1;i<=55;i++) {
                                ma[i] -= ma[1+(i+30) % 55];
                                if (ma[i] < MZ) ma[i] += MBIG;
                        }
                inext=0;
                inextp=31;
                *idum=1;
        }
        if (++inext == 56) inext=1;
        if (++inextp == 56) inextp=1;
        mj=ma[inext]-ma[inextp];
        if (mj < MZ) mj += MBIG;
        ma[inext]=mj;
        return mj*FAC;
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC


float ran4(long *idum)
{
        void psdes(unsigned long *lword, unsigned long *irword);
        unsigned long irword,itemp,lword;
        static long idums = 0;
#if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
        static unsigned long jflone = 0x00004080;
        static unsigned long jflmsk = 0xffff007f;
#else
        static unsigned long jflone = 0x3f800000;
        static unsigned long jflmsk = 0x007fffff;
#endif

        if (*idum < 0) {
                idums = -(*idum);
                *idum=1;
        }
        irword=(*idum);
        lword=idums;
        psdes(&lword,&irword);
        itemp=jflone | (jflmsk & irword);
        ++(*idum);
        return (*(float *)&itemp)-1.0;
}

#define NITER 4

void psdes(unsigned long *lword, unsigned long *irword)
{
        unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
        static unsigned long c1[NITER]={
                0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
        static unsigned long c2[NITER]={
                0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};

        for (i=0;i<NITER;i++) {
                ia=(iswap=(*irword)) ^ c1[i];
                itmpl = ia & 0xffff;
                itmph = ia >> 16;
                ib=itmpl*itmpl+ ~(itmph*itmph);
                *irword=(*lword) ^ (((ia = (ib >> 16) |
                        ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
                *lword=iswap;
        }
}
#undef NITER

