//*****************************************************************************************//
// The program simulates the performance of RM(r,m)^2-turbo code using the trellis-based MAP
// decoding  algorithm for linear block codes. The minimal trellis was drawn using Massey 
// Construction. Five decoding iteration is set and BER is written in file:BTC_General.dat
// User can input RM code parameters, i.e., r and m of RM(r,m) code. For example, RM(3,5) code
// is RM code with code length = 32 and information length = 26.  This program works up to 
// RM(4,6) or code length = 64 and information length = 57. 
// This program is written by Ms. Bo Yin
//*****************************************************************************************//

#include <iostream.h>
#include <stdlib.h>
#include <iomanip.h>
#include <fstream.h>
#include <math.h>
#include "ECHELON64.h"
#include "RM_64.h"
#include "Encoder_RM64.h"

class Node {
 public:
  Node();
  int getTime(void);
  double getStage(void);
  void setStage(double s);
  void setTime(int t);
  void setForwardmatric(double fm);
  void setBackwardmatric(double bm); 
  double getForwardmatric(void);
  double getBackwardmatric(void);
 private: 
  int time;
  double stage;
  double forwardmatric;
  double backwardmatric;
};

Node::Node() {time=0; stage=0.0; forwardmatric = 0.0; backwardmatric = 0.0;}

int Node::getTime(void) {return time;}
double Node::getStage(void) {return stage;}
double Node::getForwardmatric(void) {return forwardmatric;}
double Node::getBackwardmatric(void) {return backwardmatric;}
void Node::setStage(double s) {stage = s;}
void Node::setTime(int t) {time = t;}
void Node::setForwardmatric(double fm) {forwardmatric = fm;}
void Node::setBackwardmatric(double bm) {backwardmatric = bm;} 

class Branch {
public:
  Branch();
  int getOriginaltime(void);
  double getOriginal(void);
  double getDestination(void);
  int getBit(void); // the bit for each branch
  int getOriginalnode(void);// get the index of each original node
  int getDestinationnode(void);// get the index of each destination node
  void setOriginaltime(int ot);
  void setOriginal(double o);
  void setDestination(double d);
  void setBit(int b);//set the bit for each branch
  void setOriginalnode(int on);//set the index of each original node 
  void setDestinationnode(int dn);//set the index of each destination node
private:
  int originaltime;
  double original;
  int originalnode; 
  int destinationnode;
  double destination; 
  int bit;
};

Branch::Branch() {originaltime=bit=0; original=destination=0.0;}
int Branch::getOriginaltime(void) {return originaltime;}
double Branch::getOriginal(void) {return original;}
double Branch::getDestination(void) {return destination;}
int Branch::getBit(void) {return bit;}
int Branch::getOriginalnode(void) {return originalnode;}
int Branch::getDestinationnode(void) {return destinationnode;}

void Branch::setOriginaltime(int ot) {originaltime = ot;}
void Branch::setOriginal(double o) {original = o;}
void Branch::setDestination(double d) {destination = d;}
void Branch::setBit(int b) {bit = b;}
void Branch::setOriginalnode(int on) {originalnode = on;}
void Branch::setDestinationnode(int dn) {destinationnode = dn;}

/*******************************************************/
/*                    Main part                        */
/*******************************************************/

int findbit(int originaltime,double original,int *index, int Gx[][65], int RM_N, int RM_K);
void newAddXYtoX(int *X,int *Y,int leng, int RM_N);
double findstate(int *index,int length, int Gx[][65], int RM_N, int RM_K);
void fillid(int *in,int *out,int s,int a, int RM_N);
fstream tout("BTC_General.dat", ios::out);

int main()
{

  int lastbranch, lastnode;
  int numoneindex = 0;
  int numindex = 0;
  int startnodeindex[65];
  int startbranchindex[65];
  int i,j;
  int countnode = 0 ,numbranch = 0;
  int r, m, k, n;
  int Left[58];
  int original_G[58][65], G[58][65];
  const int *G1ptr, *G2ptr;
  int statecomplexfile[65], *stateptr;
  int nodenumber = 0; 
  int *L;
  
  cout << "Enter (r, m) : ";
  cin >> r >> m;
  RM_64 mycode(m, r);// the definition of an object for RM code original G_matrix
  G1ptr = mycode.get_G();
  
  k = mycode.get_k();
  cout << "Information length=" << k << endl;
  tout << "Information length=" << k << endl;
  n = mycode.get_n();
  cout << "Code length =" << n << endl;
  tout << "Code length =" << n << endl << endl;
  for ( i = 0; i < k; i++)
    {
      for ( j = 0; j < n; j++)
	original_G[i][j] = G1ptr[65*i + j];
    }
 
  
  ECHELON64 echelon_G( original_G, k, n);//the definition of an object for Echelon G_matrix
  G2ptr = echelon_G.get_echelon();
  for ( i = 0; i < k; i++)
    { for ( j = 0; j < n; j++)
      G[i][j] = G2ptr[65*i + j];}
 
  echelon_G.display_MSGM(); 
  stateptr = echelon_G.statecomplex();
  for (i = 0; i <= n; i++)
    statecomplexfile[i] = *(stateptr + i);
  for ( i=0; i<=n; i++)
    for (i = 0; i < n; i++)
      countnode += ( int) pow (2, statecomplexfile[i]);
  
  Node node[3000];
  Branch branch[6000];
  for(i=0; i < n + 1; i++)
    {
      for(j=0; j < pow(2,statecomplexfile[i]); j++)
	{
	  node[nodenumber].setTime(i);
	  nodenumber++;
	}
    }
  
  startnodeindex[0] = 0;
  startnodeindex[1] = 1;
  for(i=2; i < n + 1; i++)
    startnodeindex[i] = startnodeindex[i-1] + (int)pow(2,statecomplexfile[i-1]);
  
  startbranchindex[0] = 0;
 for(i=1; i < n; i++)
   {
     if(statecomplexfile[i-1] < statecomplexfile[i])
       startbranchindex[i] = startbranchindex[i-1] + (int)pow(2,statecomplexfile[i]);
     else if(statecomplexfile[i-1] == statecomplexfile[i])
       startbranchindex[i] = startbranchindex[i-1] + 2*(int)pow(2,statecomplexfile[i]);
     else if (statecomplexfile[i-1] > statecomplexfile[i])
       startbranchindex[i] = startbranchindex[i-1] + (int)pow(2,statecomplexfile[i-1]);
   }
 startbranchindex[n] = 0;

 int id[2638][10];
 int dim,s;


 for(i = 0; i < 2638; i++) 
   for(j = 0; j < 10; j++)
     id[i][j] = -1;
 
 id[0][0] = -2;

 nodenumber = 0;
 numbranch  = 0;

 
 
 L =  echelon_G.get_Left();
 for ( i = 0; i < k; i++)
   Left[i] = *(L +i);
 

/*******************************/
/*      Trellis part           */
/*******************************/
 
 for(i=0; i < n; i++)
   {
     for(dim=0; dim < k; dim++)
	{
	  if(i == Left[dim])
	    {s = dim; break;}
	  else {s = -1;}
	}
	
     if(s >= 0)
       {
	 int nn=0;
	 for(int kk=0;kk < pow(2,statecomplexfile[i]); kk++)
	   {
	     if(kk == 0)
	       {
		 for(j=0; j < 2; j++)
		   {
		     branch[numbranch].setOriginaltime(i);
		     branch[numbranch].setOriginal(findstate(id[nodenumber],n-i, G,n,k));
		     branch[numbranch].setOriginalnode(nodenumber);
		   
		     if(j==0)
		       {
			 branch[numbranch].setDestinationnode(startnodeindex[i+1]+2*kk+j);
			 fillid(id[branch[numbranch].getOriginalnode()],id[branch[numbranch].getDestinationnode()],s,-1,n);
			 branch[numbranch].setDestination(findstate(id[startnodeindex[i+1]+2*kk+j],n-i-1,G,n,k));
			 branch[numbranch].setBit(findbit(i,branch[numbranch].getOriginal(),id[startnodeindex[i+1]+2*kk+j], G, n, k));
			 node[startnodeindex[i+1]+2*kk+j].setStage(branch[numbranch].getDestination());
			 nn++;
		       }
		     else
		       {
			 branch[numbranch].setDestinationnode(startnodeindex[i+1]+2*kk+j);
			 fillid(id[branch[numbranch].getOriginalnode()],id[branch[numbranch].getDestinationnode()],s,1,n);
			 branch[numbranch].setDestination(findstate(id[startnodeindex[i+1]+2*kk+j],n-i-1,G,n,k));
			 branch[numbranch].setBit(findbit(i,branch[numbranch].getOriginal(),id[startnodeindex[i+1]+2*kk+j], G, n,k));
			 node[startnodeindex[i+1]+2*kk+j].setStage(branch[numbranch].getDestination());
			 nn++;
		       }
		     numbranch++;
		   }// loop j
		 nodenumber++;
	       } // if(kk==0) statement
	     else
	       {
		 int l,b, arrayvalue[65];
		 double value;
		 char found;
		 for(j=0; j < 2; j++)
		   {
		     branch[numbranch].setOriginaltime(i);
		     branch[numbranch].setOriginal(findstate(id[nodenumber],n-i,G,n,k));
		     branch[numbranch].setOriginalnode(nodenumber);
		     if(j==0)
		       {
			 
			 value = branch[numbranch].getOriginal();
			 for(l=n-i-1; l >= 0; l--)
			   {
			     arrayvalue[l] = (int)fmod(value,2);
			     value = value/2;
			   }
			 
			 value = 0.0;
			 for(l=n-1-i; l >= 1; l--)
			   {
			     if (arrayvalue[l]==1)
			       value +=  pow(2,n-1-i-l);
			   }
		       }
		     else
		       {
			 int idtemp[65];
			 for(int p = 0; p < n; p++)
			   idtemp[p] = -1;
			 fillid(id[branch[numbranch].getOriginalnode()],idtemp,s,1,n);
			 value = findstate(idtemp,n-i-1,G,n,k);
			 
		       }
		    
		     for(b = 0; b < 2 * pow(2, statecomplexfile[i]); b++)
		       {
			double temp;
			temp = branch[startbranchindex[i] + b].getDestination();
			 if(value == temp)
			   {
			     branch[numbranch].setDestination(temp);
			     branch[numbranch].setDestinationnode(branch[startbranchindex[i]+b].getDestinationnode());
			     branch[numbranch].setBit(branch[startbranchindex[i]+b].getBit()*-1);
			     found = 'y';
			     break;
			   }
			 else
			   {
			     found = 'n';
			   }
		       } // loop b
		     
		     
		     if(found == 'n')
		       {
			 branch[numbranch].setDestinationnode(startnodeindex[i+1]+nn);
			 if(j==0)
			   {
			     fillid(id[branch[numbranch].getOriginalnode()],id[branch[numbranch].getDestinationnode()],s,-1,n);
			   }
			 else
			   {
			     fillid(id[branch[numbranch].getOriginalnode()],id[branch[numbranch].getDestinationnode()],s,1,n);
			   } 
			 branch[numbranch].setDestination(findstate(id[startnodeindex[i+1]+nn],n-i-1,G,n,k));
			 branch[numbranch].setBit(findbit(i,branch[numbranch].getOriginal(),id[startnodeindex[i+1]+nn], G,n,k));
			 node[startnodeindex[i+1]+nn].setStage(branch[numbranch].getDestination());
			 nn++;
		       }// if(found == 'n')
		     
		     numbranch++;
		     
		   } //loop j   
		 nodenumber++;
	       } // from else statement of (kk > 0)
	     // nodenumber++;
	   }// loop kk  
       }
     else
       {
	 int nn=0;
	 for(int kk=0; kk < pow(2,statecomplexfile[i]); kk++)
	   {
	     branch[numbranch].setOriginaltime(i);
	     branch[numbranch].setOriginal(findstate(id[nodenumber],n-i,G,n,k));
	     branch[numbranch].setOriginalnode(nodenumber);
	     if(kk==0)
	       {
		 branch[numbranch].setDestinationnode(startnodeindex[i+1]+kk);
		 fillid(id[branch[numbranch].getOriginalnode()],id[branch[numbranch].getDestinationnode()],s,-1,n);
		 branch[numbranch].setDestination(findstate(id[startnodeindex[i+1]+kk],n-i-1,G,n,k));
		 branch[numbranch].setBit(findbit(i,branch[numbranch].getOriginal(),id[startnodeindex[i+1]+kk],G,n,k));
		 node[startnodeindex[i+1]+kk].setStage(branch[numbranch].getDestination());
		 nn++;
	       }
	     else
	       {
		double value;
		int l,b, arrayvalue[65];
		char found;
		 value = branch[numbranch].getOriginal();
		 for(l=n-i-1; l >= 0; l--)
		   {
		     arrayvalue[l] = (int)fmod(value, 2);
		     value = value / 2;
		   }
		      
		 value = 0.0;
		 for(l = n-1-i; l >= 1; l--)
		   {
		     if (arrayvalue[l] == 1)
		       value += pow(2, n-1-i-l);
		   }
		 
		 for(b = 0; b < pow(2, statecomplexfile[i]); b++)
		   {
		     double temp;
		     temp = branch[startbranchindex[i]+b].getDestination();
		     if(value == temp)
		       {
			 branch[numbranch].setDestination(temp);
			 branch[numbranch].setDestinationnode(branch[startbranchindex[i]+b].getDestinationnode());
			 branch[numbranch].setBit(branch[startbranchindex[i]+b].getBit()*-1);
			 found = 'y';
			 break;
		       }
		     else
		       {
			 found = 'n';
		       }
		   }
		 
		 if(found == 'n')
		   {
		     branch[numbranch].setDestinationnode(startnodeindex[i+1]+nn);
		     fillid(id[branch[numbranch].getOriginalnode()],id[branch[numbranch].getDestinationnode()],s,-1,n);
		     branch[numbranch].setDestination(findstate(id[startnodeindex[i+1]+nn],n-i-1,G,n,k));
		     branch[numbranch].setBit(findbit(i,branch[numbranch].getOriginal(),id[startnodeindex[i+1]+nn],G,n,k));
		     node[startnodeindex[i+1]+nn].setStage(branch[numbranch].getDestination());
		     nn++;
		   }
	       }
	     numbranch++;
	     nodenumber++;
	   }  
       }
   }
 
 
 lastbranch = numbranch;
 lastnode = nodenumber + 1;
 double snr, Rdb, SNR[10] = {0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5}, sigma, Lc, R, EbN0, Lxy[65][65], Lu[65][65], Le_H[65][65], Le_V[65][65], LL_U[65][65];
 double *Y_H_Vptr;
 int *Uptr, U[58][58], U_estimate[58][58];
 double Y_H_V[65][65];
 int counterbranch = 0, nodenum, counter;
 double temp1,temp2, t, tt;
 int d ;
 int times;
 int error;
 double Pe, loop;
 int iteration = 5, iter;
 
 R = ( double) k*k /(n*n - (n - k)*(n - k)) ;
 tout << "Code Rate =  " << R << endl;
 Rdb = 10*log10(R);
 
 for (int w = 0; w < 8; w++){
   
   times = 0; error = 0; 
   
   EbN0 = pow (10, SNR[w]/10);
   sigma = sqrt(1.0/(2.0*R*EbN0));
   tout << "EBN0 = " << SNR[w] << endl;
   Lc = 4*R*EbN0;
   Encoder_double64 encoder(G, k, n, sigma);
   
   for ( i = 0; i < n; i++) // initialize the Lu[i][j] 
     for (j = 0; j < n; j++)
       Lu[i][j] = Le_V[i][j] = 0;
 
   if ( SNR[w] == 0.0) 
   loop = 2e+0;
   else if ( SNR[w] == 0.5) 
   loop = 4e+0;