/*********************************************************
	MAND
	source code and documentation copyright 1995 and 1996
	Stephan Pelikan, Steven H. Rogstad  and
	the University of Cincinnati.

	We're giving this program away.
	You alone must decide if it is correct
	and suitable for your purposes. If you do use
	it, you do so at your own risk.

	You may redistribute the source code or executable
	versions of this program provided that you
		  1) include this copyright notice,
		  2) do not charge money for it, and
		  3) distribute the program without modifications.

**********************************************************/
/*
CHANGES and HISTORY

September 1997
Added a new routine that checks all possible pairs of
parents for each lane of offspring and lists potential pairs. This routine
is used if there's two groups (0=offspring,1=potential parents)

We should include some sort of probabilistic calculutions.
*/

//#include <iostream.h>  //not needed with include of fstream.h
#include <fstream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include "arrays.h"
//#include "maths.h"
#include "llist.h"
#define MIN(A,B) ((A <= B)? A:B)
#define PROGRAM_NAME "MAND"
#define VERSION "0.2"
#define BUILD_DATE "25 September 97"
/***********************************

***********************************/




/***********************************
some global variables where data is stored
**************************************/
int N,G,P;
imatrix data;

/*
iVector number_of_bands;
matrix frequency_of_bands;
Vector avg_number_bands;
iVector number_in_group;
*/



/***********************************
FORMAT OF DATA FILE
There are N "population" bands and P lanes.
The data is grouped into G groups.
The data is stored as an ascii file with N+1
rows and P columns that follow a line on which
N, P, and G appear (in that order). The first row is used
to specify which group the lane belongs to.
You must include the first row. If you don't
want to define subgroups, make G=1 and the first row
all 0's. The groups are identified by consecutive
integers starting with 0 and ending with G-1.

For this program there are three groups, offspring (group 0)
maternal parent (group 1) and potential paternal parent (group 2).

All the program does is this:

For each offspring, it lists all possible paternal parents by 
finding all those lanes in group 2 with the property that all bands in the
offspring are represented in the maternal + paternal lanes.
***********************************/

/*
read_data() is called by process_command_line()
It reads in data from fstream f,
tests it for validity, stores the dimensions of the data
(N,P,G) in global variables, and returns an imatrix
of the data.
*/
imatrix read_data(fstream& f,int& N,int& P,int& G)
{
	int i,j;
	f >> N >> P >> G;
	imatrix data(N+1,P);
	for(i=0;i<N+1;i++)
	{
		for(j=0;j<P;j++)
		{
			f >> data[i][j];
			if(f.eof())
			{
				cerr << "Problem in read_data(): can't seem to read enough data" << endl;
				cerr << "EOF occurred before entire data array was read." << endl;
				exit(37);
			}
		}
	}
	cout << "Here's the data I've just read\n";
	cout << "\n\n-----------------------------------------------\n";
	for(i=0;i<N+1;i++)
	{
		for(j=0;j<P;j++)
		{
			cout << setw(2) << data[i][j] ;
		}
		cout << endl;
	}
	cout << "\n\n-----------------------------------------------\n";
	cout << "Number of bands = " << N << "\nNumber of lanes = " << P << "\nNumber of groups = " << G << endl;
	cout << "\n\n";

//Check to make sure the main array has only 0's and 1's
	for(i=1;i<N+1;i++)
	{
		for(j=0;j<P;j++)
		{
			if( data[i][j] <0 || data[i][j] > 1)
			{
				cerr << "Problem: your data file contains values other that 0 and 1\n";
				cerr << "For example, data["<< i << "][" << j << "] = " << data[i][j] << endl;
				cerr << "Check this.\n";
				exit(7);
			}
		}
	}
/// Check to make sure labels for groups in first row of data array are okay.
	for(j=0;j<(P);j++)
	{
		if(data[0][j] <0 || data[0][j] >(G)-1)
		{
			cerr << "Problem: your data file mentions more groups than expected\n";
			cerr << "With G = " << G << " only the values 0,...," << (G-1) << " should appear in\n";
			cerr << "the first row of data. But I found data[0][" << j << "]=" << data[0][j] << endl;
			cerr << "You'd better fix it before running the program on this data\n";
			exit(7);
		}
	}
	return data;
}


void process_command_line(int argc,char *argv[])
{
	cerr << PROGRAM_NAME << " version " << VERSION << " built on " << BUILD_DATE << endl;
	if(argc < 2)
	{
		cout << "\n----------------------------------------------------------\n";
		cout << "Use this program by giving the command "<< endl;
		cout << argv[0] << " <datafilename>" << endl;
		cout << "\n----------------------------------------------------------\n";
		exit( 0);
	}
	fstream fin(argv[1],ios::in);

	if( fin.fail())
	{
		cout << "Sorry: can't seem to open file " << argv[1] << endl;;
		exit(11);
	}
// 	if(argc==3)
// 	{
// 		if(isdigit(argv[2][0])) Iterations = atoi(argv[2]);
// 		else
// 		{
// 			cerr << "You seem to be trying to set the number of iterations, but\n";
// 			cerr << "I can't treat the third argument, " << argv[2] << ", on command line as an integer\n";
// 		}
// 	}
	cout << "Output from program " << PROGRAM_NAME << " version " << VERSION << " built " << BUILD_DATE << endl;;

	cout << "Reading data from file: " << argv[1] << endl;

	data=read_data(fin,N,P,G);
	fin.close();
	return;
}
int main(int argc, char *argv[])
{
  int l1,l2,l3,misses,band,okay;
  LList<int> L;
	cout << setiosflags(ios::showpoint| ios::fixed| ios::left) ;
	cout << setprecision(4);

	
	process_command_line(argc,argv) ;

if(G==3)
{

	/*	
	  for each lane l1  in group 0
	  for each lane l2 in group 1
	  for each lane l3  in gropup 2
	  count the bands in l1 that don't  appear in either l2 or l3 and output
	  "l1 l2 l3 count " followed by a list of band numbers that don't appear in l2 or l3
	  */
	for(l1 =0; l1 < P; l1++)
	  {
		 if( data[0][l1] == 0)
	      {
		for(l2=0;l2<P;l2++)
		  {
		    if( data[0][l2]==1)
		      {
		
			for(l3=0;l3<P;l3++)
			  {
			    L.Clear();
			    if(data[0][l3]==2)
			      {
				misses=0;
				band=1;
				while(band <=N )
				  {
				    if(data[band][l1] ==1 && ( data[band][l2]==0 && data[band][l3]==0))
				     {
				       misses++;
				       L.Insert(band);
					  }
				    band++;
				  }
			    
				cout << l1 <<" " << l2 << " " << l3 << " " << misses << " ";
				if(misses <4) cout << L << endl; else cout <<endl;
				if(misses==0) cout << "<<<<<<<<-----------" << endl;
			      }
			      
			  }	
		      }
		  }
	      }
	  }
}
if(G==2)
{
	//for each lane in group 0, consider all possible pairs
	//from group one (including self crosses) and list those that
	//could be parents
  cout << "For each offspring band (group 0) we list all possible pairs" <<endl;
  cout <<"of lanes from group 1 that could be parents" << endl;
	for(l1=0;l1<P;l1++)
	{
		if(data[0][l1]==0) //if l1 is group 0
		{
		cout <<"\nOffspring lane " << l1 <<": ";
			for(l2=0;l2<P;l2++)
			{
				if(data[0][l2]==1)
				{
					for(l3=l2; l3<P;l3++)
					{
						if(data[0][l3]==1)
						{
							//test
							okay=1;
							for(band=1;band<=N;band++) 
							  {
								if(data[band][l1]==1 && 
									data[band][l2]==0 &&data[band][l3]==0)
								{ 	
									okay=0;
									break;
								}
							  }
							if(okay==1) cout <<"("<<l2<<","<<l3<<") ";
						}
					}
				}
			}
		}
	}
}      
	cout << "DONE" << endl;
	return 1;
}