// 01/24/2007 fix the problem of overflow of inTab arrays (previous define to size of 500), 
//            modified to dynamical number of residue number + 1
#include "SPARTA.h"

SPARTA::SPARTA()
{	
	string libvar;
	if( getenv( "SPARTA_DIR" ) == NULL )	{
		SPARTA_DIR = ".";
		string temp; 
		if( getenv( "PATH" ) != NULL) {
			temp = getenv( "PATH" );
			if(temp.find("/") != string::npos ) slash_char = "/"; // unix
			else if(temp.find("\\") != string::npos ) slash_char = "\\"; // Windows 
		}
		else slash_char = "\\"; //default Windows
	}
	else 
	{
		SPARTA_DIR = getenv( "SPARTA_DIR" );
		if(SPARTA_DIR.find("/") != string::npos ) slash_char = "/"; // unix
		else if(SPARTA_DIR.find("\\") != string::npos ) slash_char = "\\"; // Windows 
		else SPARTA_DIR = ".";
	}

	argList.push_back("spartaDir");
	argList.push_back("shiftDir");
	argList.push_back("pdbDir");
	argList.push_back("tabDir");
	argList.push_back("predDir");
	argList.push_back("in");
	argList.push_back("sum");
	argList.push_back("out"); //same as "-sum"
	argList.push_back("homo");
	argList.push_back("ref");
	argList.push_back("db");
	argList.push_back("w");
	argList.push_back("fitting");
	argList.push_back("pdbList");
	argList.push_back("rand");
	argList.push_back("adj");
	argList.push_back("prev");
	argList.push_back("next");
	argList.push_back("atom");
	argList.push_back("excl");
	argList.push_back("first");
	argList.push_back("last");
	argList.push_back("count");
	argList.push_back("compile");
	argList.push_back("help");
	argList.push_back("time");
	argList.push_back("notime");

	aN[1]="N"; aN[2]="HA"; aN[3]="C"; aN[4]="CA"; aN[5]="CB"; aN[6]="H";
	aN_ALL[1]="N"; aN_ALL[2]="HA"; aN_ALL[3]="C"; aN_ALL[4]="CA"; aN_ALL[5]="CB"; aN_ALL[6]="H";
}


SPARTA::SPARTA(const string &fileName)
{
	SPARTA();
	inName = fileName;
}


void SPARTA::printSyntax()
{
	cerr << "SPARTA: Shifts Prediction from Analogue of Residue type and Torsion Angle." << endl;
	cerr << "Version 2008.01.22 Shen and Bax" << endl << endl;

	cerr << "Predicts chemical shifts from angle and sequence homology " << endl;
	cerr << "Directory Settings: " << endl;
	cerr << " -spartaDir [.]               SPARTA main directory: SPARTA_DIR." << endl;
	cerr << " -tabDir    [tab]             Subdirectory of master tables." << endl;
	cerr << " -predDir   [pred]            Directory for prediction results." << endl;
	cerr << " -shiftDir  [shifts]          Subdirectory of master shifts (For compiling database only)." << endl;
	cerr << " -pdbDir    [pdb]             Subdirectory of master structures (For compiling database only)." << endl;
	cerr << "File Name Settings: " << endl;
	cerr << " -in        [None]            Input PDB coordinate file." << endl;
	cerr << " -sum       [pred.tab]        Prediction summary output." << endl;
	cerr << " -db        [sparta.tab]      SPARTA triplet database file." << endl;
	cerr << " -w         [weight.tab]      SPARTA triplet weighting file." << endl;
	cerr << " -homo      [homology.tab]    Residue homology table." << endl;
	cerr << " -fitting   [fitting.tab]     Fitting parameters to calculate predict error." << endl;
	cerr << " -ref       [None]            Refernce(Observed) chemical shift table." << endl;
	cerr << " -pdbList   [None]            Candidate protein name table (For compiling database only)" << endl;
	cerr << "Random Coil and Secondary Shift Tables:" << endl;
	cerr << " -rand      [randcoil.tab]    Random coil shift table." << endl;
	cerr << " -adj       [rcadj.tab]       Adjustment table." << endl;
	cerr << " -prev      [rcprev.tab]      Previous residue adjustment." << endl;
	cerr << " -next      [rcnext.tab]      Next residue adjustment." << endl << endl;
	cerr << "Other Options:" << endl;
	cerr << " -atom      [N HA C CA CB H]  List of shift types to predict. (Default: All)" << endl;
	cerr << " -excl      [None]            List of proteins to exclude." << endl;
	cerr << " -first     [First]           First RESID to use." << endl;
	cerr << " -last      [Last]            Last RESID to use." << endl;
	cerr << " -count     [20]              Max match count per triplet." << endl;
	cerr << "Examples:" << endl;
	cerr << " SPARTA -in bpti.pdb          Form shifts predictions from a given PDB coordinate file." << endl;
	cerr << " SPARTA -compile -pdbList list.tab " << endl ;
	cerr << "                              Compile new SPARTA database using proteins named in list.tab." << endl;

	exit(0);
}


//preset the args form command line SHIFT_DIR
void SPARTA::setArgs( map<string, string> args )
{
	
	if( args["help"].length() > 0 ) printSyntax();
	
	//Directory setting
	if( args["spartaDir"].length() > 0 ) SPARTA_DIR = args["spartaDir"];

	
	TAB_DIR = SPARTA_DIR + slash_char+ "tab";
	if( args["tabDir"].length() > 0 ) TAB_DIR = args["tabDir"];
	
	SHIFT_DIR = SPARTA_DIR + slash_char+ "shifts";
	if( args["shiftDir"].length() > 0 ) SHIFT_DIR = args["shiftDir"];

	PDB_DIR = SPARTA_DIR + slash_char+ "pdb";
	if( args["pdbDir"].length() > 0 ) PDB_DIR = args["pdbDir"];

	PRED_DIR = "pred";
	if( args["predDir"].length() > 0 ) PRED_DIR = args["predDir"];	


	//File names setting
	if( args["in"].length() > 0 ) inName = args["in"];

	tripFileName = TAB_DIR + slash_char+ "sparta.tab";
	if( args["db"].length() > 0 ) tripFileName = args["db"];
	
	weightFileName = TAB_DIR + slash_char + "weight.tab";
	if( args["w"].length() > 0 ) weightFileName = args["w"];
	
	homoFileName = TAB_DIR + slash_char + "homology.tab";
	if( args["homo"].length() > 0 ) homoFileName = args["homo"];
	
	fitFileName = TAB_DIR + slash_char + "fitting.tab";
	if( args["fitting"].length() > 0 ) fitFileName = args["fitting"];

	sumName = PRED_DIR + slash_char +"pred.tab";
	if( args["sum"].length() > 0 ) sumName = args["sum"];
	if( args["out"].length() > 0 ) sumName = args["out"];

	if( args["ref"].length() > 0 ) refCSFileName = args["ref"];
	if( args["pdbList"].length() > 0 ) pdbListName = args["pdbList"];

	// Random Coil and Secondary Shift Tables
	rcFileName = TAB_DIR + slash_char + "randcoil.tab";
	if( args["rc"].length() > 0 ) rcFileName = args["rc"];
	
	adjFileName = TAB_DIR + slash_char + "rcadj.tab";
	if( args["adj"].length() > 0 ) adjFileName = args["adj"];
	
	prevFileName = TAB_DIR + slash_char + "rcprev.tab";
	if( args["prev"].length() > 0 ) prevFileName = args["prev"];
	
	nextFileName = TAB_DIR + slash_char + "rcnext.tab";
	if( args["next"].length() > 0 ) nextFileName = args["next"];
	


	//Other Options
	EXCLUDED="";
	if(args["excl"].length() > 0 ) EXCLUDED = args["excl"];
	
	if(args["atom"].length() > 0 ) {
		aN.clear();
		vector<string> temp = GDB::split(" ", args["atom"]);
		int cnt = 1;
		for(int i = 0; i < temp.size(); i++) 
		{
			if( temp[i]!="N" && temp[i]!="HA"&& temp[i]!="C"&& temp[i]!="CA"&& temp[i]!="CB"&& temp[i]!="H")
			{	
				cerr << "\tInvalid atom -" << temp[i] << endl;
				exit(0);
			}
			aN[cnt++] = temp[i];
		}
	}

	matchCount = 20;
	if( args["count"].length() > 0 ) matchCount = atoi(args["count"].c_str());
	if( matchCount == 0 ) matchCount = 20;
	
	tVal = 500.0; // not used
	if( args["thresh"].length() > 0 ) tVal = atof(args["thresh"].c_str());
	if( tVal == 0.0 ) tVal = 500.0;
	
	firstRes = -9999;
	if( args["first"].length() > 0 ) firstRes = atoi(args["first"].c_str());
	//if( firstRes == 0 ) firstRes = 1;

	lastRes = 9999;
	if( args["last"].length() > 0 ) lastRes = atoi(args["last"].c_str());
	//if( lastRes == 0.0 ) lastRes = 9999;

}



void SPARTA::init()
{
	cout << "Reading PDB Coordinates from " << inName << endl;
	inPDB.loadPDB(inName);
	
	residList = inPDB.residListOne;

	r1 = inPDB.r1;
	rN = inPDB.rN;

	if (firstRes < r1) firstRes = r1;
    	if (lastRes  < r1) lastRes = r1;

	if (firstRes > rN) firstRes = rN;
    	if (lastRes  > rN) lastRes = rN;

    	if (firstRes > lastRes) 
	{
        	int itemp = firstRes;
        	firstRes  = lastRes;
        	lastRes   = itemp;
	}

	cout << "Reading Random Coil Shifts from " << rcFileName << endl;
	RC_Tab.loadGDB( rcFileName );

	cout << "Reading RC Adjustments from " << adjFileName << endl;
	ADJ_Tab.loadGDB( adjFileName );

	cout << "Reading Previous Residue RC Adjustments from " << prevFileName << endl;
	PREV_Tab.loadGDB( prevFileName );

	cout << "Reading Next Residue RC Adjustments from " << nextFileName << endl;
	NEXT_Tab.loadGDB( nextFileName );

	// init Weight
	cout << "Reading Weighting Factors from " << weightFileName << endl;
	WEIGHT_Tab.loadGDB(weightFileName );
	
	// init homology talble
	cout << "Reading Residue Homology Table from " << homoFileName << endl;
	HOMO_Tab.loadGDB( homoFileName );
	// init Triplet table

	// init fitting parameter talble
	cout << "Reading Fitting Parameter Table from " << fitFileName << endl;
	FIT_Tab.loadGDB( fitFileName );
	// init Triplet table

	cout << "Reading " << tripFileName << ", " ;
	TRIPLET_Tab.loadGDB( tripFileName );
	//TRIPLET_Tab.loadGDB( tripFileName, "SOURCE", EXCLUDED );
	cout << TRIPLET_Tab.getEntryCount() << " Triplets" << endl;

}


void SPARTA::runPredict()
{
	init();

	map<int, string>::iterator itN;

	// mkdir for prediction
	if (PRED_DIR.find_last_of(slash_char) == PRED_DIR.length()-1 ) 
		PRED_DIR = PRED_DIR.substr(0,PRED_DIR.length()-1);
	for( itN = aN.begin(); itN != aN.end(); itN++ )
		mkdir_pred(PRED_DIR+slash_char+itN->second);

	//get residue and structure information from input coordinates
	cout << "Analyzing " << inName << " " ;
	getResInfo();
	cout << residList.size() << " residues read " << endl;

	//init the prediction summary table and prediction tables
	GDB PRED_SUM;
	PRED_SUM.VARS_str_parser("  RESID RESNAME ATOMNAME SS_SHIFT SHIFT RC_SHIFT HM_SHIFT SIGMA SOURCE");
	PRED_SUM.FORMAT_str_parser(" %4d %4s %4s %9.3f %9.3f %9.3f %9.3f %9.3f %s");
	string str = "PDB "+inName+", residues "+(string)itoa(r1, buf)+" to "+ 
		(string)itoa(rN, buf)+", "+(string)itoa(rN-r1+1, buf)+" total."; 
	PRED_SUM.addRemark(str);
	str = itoa(r1, buf);
	PRED_SUM.setData("FIRST_RESID", str+"\n");
	PRED_SUM.setData("SEQUENCE", sequence);

	GDB PRED;
	PRED.VARS_str_parser("  INDEX SS_SHIFT SCORE W R1 R2 R3 SOURCE");
	PRED.FORMAT_str_parser(" %4d %9.3f %9.3f %9.3f %-4s %-4s %-4s %s");

	// allocation
	int MAX_DB_SIZE = TRIPLET_Tab.Entries.size()+1;
	int m = 10;
	float** T_ANGLES = new float* [MAX_DB_SIZE];
	//used to store the phi/psi/chi1 angles for triplets in database, for the purpose of speed
	T_ANGLES[0] = new float [MAX_DB_SIZE*m];
	for(int i = 1; i < MAX_DB_SIZE; ++i)
		T_ANGLES[i] = T_ANGLES[i-1] + m;

	float** T_SHIFTS = new float* [MAX_DB_SIZE];
	//used to store the shifts for centre residue of the triplets in database
	T_SHIFTS[0] = new float [MAX_DB_SIZE*(aN.size()+1)];
	for(int i = 1; i < MAX_DB_SIZE; ++i)
		T_SHIFTS[i] = T_SHIFTS[i-1] + aN.size()+1;

	m =4;
	string** T_NAME = new string* [MAX_DB_SIZE];
	//used to store the resNames for triplets in database
	T_NAME[0] = new string [MAX_DB_SIZE*m];
	for(int i = 1; i < MAX_DB_SIZE; ++i)
		T_NAME[i] = T_NAME[i-1] + m;
	
	// read in database tables
	// 're-format' the database entries uisng simple 2D arrays, which could save the running time
	// compared to directly iterate over the database using its origianl format (STD::map)
	map< int, map<string, string> >::iterator it, itA;
	for ( it = TRIPLET_Tab.Entries.begin(); it != TRIPLET_Tab.Entries.end(); it++ )
	{
		string source = it->second["SOURCE"];
		if( EXCLUDED.find( source ) != string::npos ) continue;

		for(int j = 1; j <=3; j++)
		{
			str = itoa(j,buf,10);
			T_ANGLES[it->first][(j-1)*3+1] = atof( (it->second)["PHI_R"+str].c_str() );
   			T_ANGLES[it->first][(j-1)*3+2] = atof( (it->second)["PSI_R"+str].c_str() );
			T_ANGLES[it->first][(j-1)*3+3] = atof( (it->second)["CHI1_R"+str].c_str() );
			for(itN = aN.begin(); itN != aN.end(); itN++)
				T_SHIFTS[it->first][itN->first] = atof( (it->second)[itN->second].c_str() );

      			T_NAME[it->first][j] = (it->second)["RESNAME_R"+str];
		}

	}//


	//Read in sequence homology table
	map<string, map<string, float> > homology;
	map<string, string>::iterator itQ;
	for ( it = HOMO_Tab.Entries.begin(); it != HOMO_Tab.Entries.end(); it++ )
	{
		string resN = it->second["RESNAME"];

		for ( itQ = it->second.begin(); itQ != it->second.end(); itQ++ )
		{
			if( itQ->first != "RESNAME" )
				homology[resN][itQ->first] = atof( itQ->second.c_str() ) ;
		}
	} //

	
	//read in weights table for atom type itN->second
	float Wts[7][15];
	for(itN = aN.begin(); itN != aN.end(); itN++) //prediction for each defined atom type
	{ 
		int ccc = 0;
		for ( it = WEIGHT_Tab.Entries.begin(); it != WEIGHT_Tab.Entries.end(); it++ )
		{
			if( it->second["RESNAME"] == itN->second ) 
				for ( itQ = it->second.begin(); itQ != it->second.end(); itQ++ ) 
					Wts[itN->first][ccc++] = atof( itQ->second.c_str() );
		}
	}//


	map<string, map<string, float> > Fitting;
	for ( it = FIT_Tab.Entries.begin(); it != FIT_Tab.Entries.end(); it++ )
	{
		Fitting[ it->second["NAME"] ]["SLOPE"] = atof( it->second["SLOPE"].c_str() );
		Fitting[ it->second["NAME"] ]["INTERCEPT"] = atof( it->second["INTERCEPT"].c_str() );
	}

	cout << "Predicting ..." << endl;
	for( itN = aN.begin(); itN != aN.end(); itN++ )
	{
		sprintf(buf, "%8s ", itN->second.c_str());
		cout << buf ;
	}
	cout << endl; 

	//Start of search for each residue from query protein
	for( it = inTab.Entries.begin(); it != inTab.Entries.end(); it++ )
	{
		int U = it->first; // index U, indentical to (residue numner of center residue - r1)

		if( U+r1 < firstRes || U+r1 > lastRes ) continue;
		
		//classify the chi1 angles for query triplet
		float chi1_s = U_ANGLES[U][6];
		if( chi1_s > -90.0 && chi1_s < -30.0) chi1_s = -60.0;
		if( chi1_s > 30.0 && chi1_s < 90.0) chi1_s = 60.0;
		if( chi1_s < 999.0 && (chi1_s > 150.0 || chi1_s < -150.0) ) chi1_s = 180.0;
		
		float chi1_p = U_ANGLES[U][3];
		if( chi1_p > -90.0 && chi1_p < -30.0) chi1_p = -60.0;
		if( chi1_p > 30.0 && chi1_p < 90.0) chi1_p = 60.0;
		if( chi1_p < 999.0 && (chi1_p > 150.0 || chi1_p < -150.0) ) chi1_p = 180.0;	
		
		float chi1_n = U_ANGLES[U][9];
		if( chi1_n > -90.0 && chi1_n < -30.0) chi1_n = -60.0;
		if( chi1_n > 30.0 && chi1_n < 90.0) chi1_n = 60.0;
		if( chi1_n < 999.0 && (chi1_n > 150.0 || chi1_n < -150.0) ) chi1_n = 180.0;		
		
		map< float, int > bList[7];

		//search the database
		for(itA = TRIPLET_Tab.Entries.begin(); itA != TRIPLET_Tab.Entries.end(); itA++)
		{
			int T = itA->first;

			if( T_ANGLES[T][6] < 999.0 && getDiff( chi1_s, T_ANGLES[T][6] ) > 60 ) continue;
			//chi1 specific for all center residue

			//if( T_SHIFT[T]> 999 || U_SOURCE[U] == T_SOURCE[T] ) continue;
			
			float Dphi_P = getDiff( U_ANGLES[U][1], T_ANGLES[T][1] )/15;
			float Dpsi_P = getDiff( U_ANGLES[U][2], T_ANGLES[T][2] )/15;
			float Dphi_S = getDiff( U_ANGLES[U][4], T_ANGLES[T][4] )/15;
			float Dpsi_S = getDiff( U_ANGLES[U][5], T_ANGLES[T][5] )/15;
			float Dphi_N = getDiff( U_ANGLES[U][7], T_ANGLES[T][7] )/15;
			float Dpsi_N = getDiff( U_ANGLES[U][8], T_ANGLES[T][8] )/15;

			float H_P = homology[ U_NAME[U][1] ][ T_NAME[T][1] ];
			float H_C = homology[ U_NAME[U][2] ][ T_NAME[T][2] ];
			float H_N = homology[ U_NAME[U][3] ][ T_NAME[T][3] ];
			
			for(itN = aN.begin(); itN != aN.end(); itN++){
				int N = itN->first;

				if( T_SHIFTS[T][N] > 999.0 ) continue;
				
				if( Wts[N][0] == 1 && T_ANGLES[T][3] < 999.0 && getDiff( chi1_p, T_ANGLES[T][3] ) > 60 ) continue;
				//chi1 specific for previous residue if required
				if( Wts[N][2] == 1 && T_ANGLES[T][9] < 999.0 && getDiff( chi1_n, T_ANGLES[T][9] ) > 60 ) continue;
				//chi1 specific for following residue if required

				float score = 0.0;

				score += Wts[N][6]*Dphi_P*Dphi_P; 
				score += Wts[N][9]*Dpsi_P*Dpsi_P; 
				score += Wts[N][7]*Dphi_S*Dphi_S; 
				score += Wts[N][10]*Dpsi_S*Dpsi_S; 
				score += Wts[N][8]*Dphi_N*Dphi_N; 
				score += Wts[N][11]*Dpsi_N*Dpsi_N;
			
				score += Wts[N][3] * H_P * H_P;
				score += Wts[N][4] * H_C * H_C;
				score += Wts[N][5] * H_N * H_N;
			
				(bList[N])[ score ] = T;
			}
		}//end of search database


		//selection of the best 20 matched triplets
		float pred[7], sum[7];
		string resName = inTab.Entries[U]["RESNAME_R2"];
		int resID = U; //residue index number, not real residue ID. -> real_resID - r1
		map< float, int >::iterator itB;
		for( itN = aN.begin(); itN != aN.end(); itN++ )
		{
			int N = itN->first;
			string aName = itN->second;
			pred[N] = 0.0; 
			sum[N] = 0.0; 
			int cnt = 0;

			// exlude the prediction for HN of Proline, CB for Gly
			if( (resName == "P" && aName == "H") || (resName == "G" && aName == "CB") )  
			{
				pred[N] = 9999.0; 
				sprintf(buf, "%8.3f ", pred[N]);
				cout << buf;
				continue;
			}

			PRED.REMARKS.clear();
			PRED.DATA.clear();
			str = aName + " Shifts Prediction Matches for Residue "+ resName +inTab.Entries[U]["RESID_R2"];
			PRED.addRemark(str);
			str = " "+(string)itoa(resID+r1-1,buf,10)+" "+(string)itoa(resID+r1,buf,10)+" "+(string)itoa(resID+r1+1,buf,10);
			PRED.setData("RESIDS  ", str);
			str = " "+inTab.Entries[U]["RESNAME_R1"]+" "+resName+" "+inTab.Entries[U]["RESNAME_R3"];
			PRED.setData("RESNAMES", str);

			vector<float> pList;

			for( itB = bList[N].begin(); itB != bList[N].end(); itB++ )
			{
				int B = itB->second;
				string source = TRIPLET_Tab.Entries[B]["SOURCE"];
				if( inTab.Entries[U]["SOURCE"] == source ) continue;

				int rID = atoi( TRIPLET_Tab.Entries[B]["RESID_R2"].c_str() );

				cnt++;
				if(cnt>matchCount) break;

	
				float S = itB->first; if( S==0 ) S = 0.00010001; // mininum value: 1/9999.0
				pred[N] += ( T_SHIFTS[B][N]/S );
				sum[N] += 1.0/S;

				pList.push_back(T_SHIFTS[B][N]);

				//add match into PRED table
				PRED.setEntry(cnt,"INDEX",itoa(cnt,buf));
				PRED.setEntry(cnt,"SS_SHIFT",ftoa(T_SHIFTS[B][N],buf) );
				PRED.setEntry(cnt,"SCORE",ftoa(itB->first,buf) );
				PRED.setEntry(cnt,"W",ftoa(1.0/itB->first,buf) );
				sprintf(buf, "%1s%-3s",TRIPLET_Tab.Entries[B]["RESNAME_R1"].c_str(), \
				TRIPLET_Tab.Entries[B]["RESID_R1"].c_str() );
				PRED.setEntry(cnt,"R1",buf);
				sprintf(buf, "%1s%-3s",TRIPLET_Tab.Entries[B]["RESNAME_R2"].c_str(), \
				TRIPLET_Tab.Entries[B]["RESID_R2"].c_str() );
				PRED.setEntry(cnt,"R2",buf);
				sprintf(buf, "%1s%-3s",TRIPLET_Tab.Entries[B]["RESNAME_R3"].c_str(), \
				TRIPLET_Tab.Entries[B]["RESID_R3"].c_str() );
				PRED.setEntry(cnt,"R3",buf);
				PRED.setEntry(cnt,"SOURCE", TRIPLET_Tab.Entries[B]["SOURCE"] );
			}
			pred[N] /= sum[N];

			// add scaled ring current shift to predicted shifts
			pred[N] += (0.6*U_RING_SHIFTS[ resID ][N]);

			// add H-Bond corrections to HA and HN predicted shifts
			if( aName == "H" && U_HN_HB[resID] > 0 ) pred[N] -= Fitting["H_HBOND"]["SLOPE"]/(U_HN_HB[resID]*U_HN_HB[resID]*U_HN_HB[resID]) + Fitting["H_HBOND"]["INTERCEPT"];
			if( aName == "HA" && U_HA_HB[resID] > 0 ) pred[N] -= Fitting["HA_HBOND"]["SLOPE"]/(U_HA_HB[resID]*U_HA_HB[resID]*U_HA_HB[resID]) + Fitting["HA_HBOND"]["INTERCEPT"];

			PRED.setData("_AVERAGE_SHIFTS_WEIGHTED",ftoa(pred[N],buf)); //the weighted predcited shifts after addition of ring shifts and H-bond corrections
			PRED.setData("_AVERAGE_SHIFTS",ftoa(getAVG(pList),buf)); // the average predicted shifts before addition of ring shifts and H-bond corrections
			PRED.setData("_STANDARD_DEVIATION",ftoa(getSTD(pList),buf)); 
			// the std.dev.of of 20 "best matched shifts" (before addition of ring shifts and H-bond shifts corrections)
			str = inTab.Entries[U]["RESID_R2"];
			if(str.length() == 1)str = "00" + str; 
			if(str.length() == 2)str = "0" + str; 
			string savName = PRED_DIR+slash_char+aName+slash_char+"res"+str+".tab";
			PRED.saveGDB(savName.c_str());


			//add items to the predcition summary table
			int index = PRED_SUM.Entries.size()+1;
			float RC = getRC(inTab.Entries[U]["RESNAME_R2"],aName);
			float RCadj = getRCadj(inTab.Entries[U]["RESNAME_R2"],aName);
			float PrevRCadj = getPrevRCadj(inTab.Entries[U]["RESNAME_R1"],aName);
			float NextRCadj = getNextRCadj(inTab.Entries[U]["RESNAME_R3"],aName);
			float p = pred[N] + RC + RCadj + PrevRCadj + NextRCadj;
			if( p > 999.0) p = 9999.0;

			PRED_SUM.setEntry(index, "RESID", inTab.Entries[U]["RESID_R2"] );
			PRED_SUM.setEntry(index, "RESNAME", inTab.Entries[U]["RESNAME_R2"] );
			PRED_SUM.setEntry(index, "SS_SHIFT", ftoa(pred[N],buf) );
			PRED_SUM.setEntry(index, "SHIFT", ftoa(p,buf) );
			PRED_SUM.setEntry(index, "SIGMA", ftoa(getSTD(pList)*Fitting[aName+"_ERROR"]["SLOPE"]+Fitting[aName+"_ERROR"]["INTERCEPT"],buf) );
			// the estimated prediction error from the std.dev.of of 20 "best matched shifts"
			PRED_SUM.setEntry(index, "RC_SHIFT", ftoa(RC+RCadj+PrevRCadj+NextRCadj ,buf) );
			PRED_SUM.setEntry(index, "HM_SHIFT", inTab.Entries[U][aName+"_HM"] );
			PRED_SUM.setEntry(index, "SOURCE", inTab.Entries[U]["SOURCE"] );
			PRED_SUM.setEntry(index, "ATOMNAME", aName );
			if( aName == "HA" && resName == "G" ) PRED_SUM.setEntry(index, "ATOMNAME", "HA2" );

			if( aName == "HA" && resName == "G" ) {
				float hm = U_RING_SHIFTS[resID][7]; //inPDB.getOrbitalShift(1, resID, "HA1"); // **************** HA1 or HA3 *************
				index++;
				PRED_SUM.Entries[index] = PRED_SUM.Entries[index-1];
				PRED_SUM.setEntry(index-1, "ATOMNAME", aName+"1" ); // **************** HA1 or HA3 *************
				PRED_SUM.setEntry(index-1, "HM_SHIFT", ftoa(hm, buf) );
				PRED_SUM.setEntry(index-1, "SHIFT", ftoa(p+0.6*(hm-U_RING_SHIFTS[resID][N]),buf) );
			}

			sprintf(buf, "%8.3f ", p);
			cout << buf;

		} // end of selection of 20 best and output
		
		sprintf(buf, "%4s %4s %8s", inTab.Entries[U]["RESID_R2"].c_str(), 
			U_NAME[U][2].c_str(), inTab.Entries[U]["SOURCE"].c_str()); 
		cout << buf << endl ;
	}
	
	PRED_SUM.saveGDB(sumName);

	cout << endl;

	if( refCSFileName.length() > 0 ) compareRef( PRED_SUM );

	// deallocation
	delete [] T_ANGLES[0];
	delete [] T_ANGLES;
	delete [] T_NAME[0];
	delete [] T_NAME;
}


//compare the OBS and PRED chemical shifts if the "-ref" option is active
//and calculate RMS deviation between OBS and PRED chemical shifts for each shifts type
void SPARTA::compareRef(GDB & Pred_Sum)
{
	int len = sumName.length();
	int pos = sumName.find_last_of(".");
	Pred_Sum.saveGDB( sumName.substr(0, pos) + "_full.tab" ); //save the original prediction summary file to a new name

	GDB REF_CS_Tab(refCSFileName);

	GDB COMP_Tab;
	
	COMP_Tab.REMARKS = Pred_Sum.REMARKS;
	COMP_Tab.DATA = Pred_Sum.DATA;
	
	COMP_Tab.VARS_str_parser("  RESID RESNAME ATOMNAME CS_OBS SHIFT RC_SHIFT CS_DIFF SIGMA W");
	COMP_Tab.FORMAT_str_parser("  %4d %4s %4s %9.3f %9.3f %9.3f %9.3f %9.3f %.2f");
	COMP_Tab.addRemark( "Observed chemical shift from: " + refCSFileName ); 

	map<int, string>::iterator itN;
	map< int, map<string, string> >::iterator it;

	vector<float> OBS_V, PRED_V, DIFF_V, OBS_V_CORRECTED;
	float obs_shift, pred_shift;

	for( itN = aN.begin(); itN != aN.end(); itN++ )
	{
		string aName = itN->second;
		
		for ( it = REF_CS_Tab.Entries.begin(); it != REF_CS_Tab.Entries.end(); it++ )
		{
			if( aName == "HA" && it->second["RESNAME"] == "G") 
			{
				if( it->second["ATOMNAME"].find("HA") != 0 ) continue;
			}
			else if( it->second["ATOMNAME"] != aName ) continue;
				
			GDB_Entry temp = Pred_Sum.getEntry("RESID",it->second["RESID"],"ATOMNAME",it->second["ATOMNAME"],1);
			
			if(temp["SHIFT"].length() <= 0) continue;

			obs_shift = atof( (it->second["SHIFT"]).c_str() );
			pred_shift = atof( (temp["SHIFT"]).c_str() );

			if( aName == "HA" && it->second["RESNAME"] == "G") //for HA of Gly
			{
						
				if( it->second["ATOMNAME"] == "HA2" ) // assign HA2 to the one with smaller shift
				{
					float shift_HA3 = 9999.000;
					GDB_Entry HA3 = REF_CS_Tab.getEntry("RESID",it->second["RESID"],"ATOMNAME","HA3",1);
					if( HA3["ATOMNAME"] == "HA3") shift_HA3 = atof( (HA3["SHIFT"]).c_str() );
					
					if( shift_HA3 < obs_shift) obs_shift = shift_HA3;
					
					shift_HA3 = 9999.000;
					HA3 = Pred_Sum.getEntry("RESID",it->second["RESID"],"ATOMNAME","HA3",1);
					if( HA3["ATOMNAME"] == "HA3") shift_HA3 = atof( (HA3["SHIFT"]).c_str() );
					
					if( shift_HA3 < pred_shift) pred_shift = shift_HA3;			
				}
				else if( it->second["ATOMNAME"] == "HA3" ) // assign HA3 to the one with larger shift
				{
					float shift_HA2 = -9999.000;
					GDB_Entry HA2 = REF_CS_Tab.getEntry("RESID",it->second["RESID"],"ATOMNAME","HA2",1);
					if( HA2["ATOMNAME"] == "HA2") shift_HA2 = atof( (HA2["SHIFT"]).c_str() );
					
					if( shift_HA2 > obs_shift) obs_shift = shift_HA2;
					
					shift_HA2 = -9999.000;
					HA2 = Pred_Sum.getEntry("RESID",it->second["RESID"],"ATOMNAME","HA2",1);
					if( HA2["ATOMNAME"] == "HA2") shift_HA2 = atof( (HA2["SHIFT"]).c_str() );
					
					if( shift_HA2 > pred_shift) pred_shift = shift_HA2;
				}
			}


			if( pred_shift > 999.0 ) continue;



			OBS_V.push_back( obs_shift );
			PRED_V.push_back( pred_shift );
			DIFF_V.push_back( pred_shift-obs_shift );

			COMP_Tab.setEntry(it->first,"RESID",it->second["RESID"]);
			COMP_Tab.setEntry(it->first,"RESNAME",it->second["RESNAME"]);
			COMP_Tab.setEntry(it->first,"ATOMNAME",it->second["ATOMNAME"]);
			COMP_Tab.setEntry(it->first,"CS_OBS",ftoa(obs_shift, buf));
			COMP_Tab.setEntry(it->first,"SHIFT",ftoa(pred_shift, buf));
			COMP_Tab.setEntry(it->first,"RC_SHIFT",temp["RC_SHIFT"]);
			// rc+adj+rxprev+rcnext, (+ HM?)
			
			COMP_Tab.setEntry(it->first,"CS_DIFF",ftoa(obs_shift-pred_shift, buf) );
			COMP_Tab.setEntry(it->first,"SIGMA",temp["SIGMA"] );
			COMP_Tab.setEntry(it->first,"W","1.00" );
		}

		string rms_str = ftoa( getRMS(OBS_V, PRED_V), buf);
		string size_str = itoa( OBS_V.size(), buf);

		string str = "RMS(OBS, PRED) for " + aName + " (n=" + size_str + "): " + rms_str +" ppm";

		COMP_Tab.setData( "NUM_"+aName, size_str);
		COMP_Tab.setData( "RMS_"+aName, rms_str);

		float avg_diff = getAVG(DIFF_V);
		float error_diff = getSTD(DIFF_V); 
		int size = DIFF_V.size();
		
		if( avg_diff >=  3*error_diff/sqrt((float)size) ) // check if the referencing correction is required
		{
		   cout << "Std.dev. of prediction error for " << aName << ": " << error_diff << " ppm" << endl;
		   cout << "Average  of prediction error for " << aName << ": " << avg_diff 
		   	<< " (> 3*" << error_diff << "/sqrt(n) )" << endl;
		   for ( it = REF_CS_Tab.Entries.begin(); it != REF_CS_Tab.Entries.end(); it++ )
		   {
			if( aName == "HA" && it->second["RESNAME"] == "G") 
			{
				if( it->second["ATOMNAME"].find("HA") != 0 ) continue;
			}
			else if( it->second["ATOMNAME"] != aName ) continue;

			obs_shift = atof( (it->second["SHIFT"]).c_str() );
			
			REF_CS_Tab.Entries[ it->first ]["SHIFT"] = ftoa( obs_shift+avg_diff, buf);
		
		   }
		}
		
		if( DIFF_V.size() > 0 )
			COMP_Tab.setData( "AVG_DIFF_"+aName, ftoa(avg_diff, buf) );
		else 
			COMP_Tab.setData( "AVG_DIFF_"+aName, "0" );
		cout << str << endl;

		for(int i = 0; i< OBS_V.size(); i++)
			OBS_V_CORRECTED.push_back( OBS_V[i] + avg_diff);
		COMP_Tab.setData( "RMS_"+aName+"_CORRECTED", ftoa( getRMS(OBS_V_CORRECTED, PRED_V), buf) );
		// the predicted shifts data after correction using <Pred-Obs>

		OBS_V.clear(); PRED_V.clear(); DIFF_V.clear(); OBS_V_CORRECTED.clear();
	}

	COMP_Tab.saveGDB( sumName );
	//COMP_Tab.saveGDB( PRED_DIR + slash_char +"comp.tab" );
	REF_CS_Tab.saveGDB( PRED_DIR + slash_char +"ref.tab" );
}



//Get the list of angles\ring shifts\h-bond information from coordinates for all possbile residues
//**************** NOT ABLE TO HANDLE PROTEIN WITH MULTIPLE CHAINS ****************
void SPARTA::getResInfo()
{
	inTab.Entries.clear();

	// allocation
	int n = rN-r1+1, m = 10;
	U_ANGLES = new float* [n];
	U_ANGLES[0] = new float [n*m];
	for(int i = 1; i < n; ++i)
		U_ANGLES[i] = U_ANGLES[i-1] + m;

	U_RING_SHIFTS = new float* [n];
	U_RING_SHIFTS[0] = new float [n*(aN.size()+1)];
	for(int i = 1; i < n; ++i)
		U_RING_SHIFTS[i] = U_RING_SHIFTS[i-1] + aN.size()+1;

	n = rN-r1+1; m = 4;
	U_NAME = new string* [n];
	U_NAME[0] = new string [n*m];
	for(int i = 1; i < n; ++i)
		U_NAME[i] = U_NAME[i-1] + m;

	U_HN_HB = new float [n];
	U_HA_HB = new float [n];
	U_CO_HB = new float [n];

	int pos0 = inName.find_last_of(slash_char)+1;
	int pos1 = inName.find_last_of(".");

	map<int, string>::iterator itN;
	int cnt = 0;
	// format the sequence read from PDB coordinates
	for(itN = residList.begin(); itN != residList.end(); itN++){
		sequence += itN->second;
		cnt++;
		if( cnt%10 == 0 ) sequence += " "; //separator for each 10 residues

		itN++; 
		if(itN != residList.end()) //add "?" if sequence numbers are not consecutive
		{
			int j = itN->first;
			itN--;
			for(int i = 1; i< j - itN->first; i++)
			{
				sequence += "?"; cnt++;
				if( cnt%10 == 0 ) sequence += " ";
			}
		}
		else itN--;
	}
	
	inPDB.initOrbitalShift();
	inPDB.initHBond();

	inTab.setData("SEQUENCE", sequence);
	inTab.VARS_str_parser("  RESID_R1 RESNAME_R1 PHI_R1 PSI_R1 CHI1_R1 RESID_R2 RESNAME_R2 PHI_R2 PSI_R2 CHI1_R2 RESID_R3 RESNAME_R3 PHI_R3 PSI_R3 CHI1_R3 N_HM HA_HM C_HM CA_HM CB_HM H_HM H_HB HA_HB CO_HB SOURCE");
	inTab.FORMAT_str_parser("%4d %s %8.3f %8.3f %8.3f %4d %s %8.3f %8.3f %8.3f %4d %s %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s");

	//calculate H-bond for the first residue
	float dist = inPDB.HBDistList[r1]["HN"];
	U_HN_HB[0] = dist;
	dist = inPDB.HBDistList[r1]["HA"];
	U_HA_HB[0] = dist;
	dist = inPDB.HBDistList[r1]["O"];
	U_CO_HB[0] = dist;	

	float shift;

	//loop for the polypeptide chain
	for ( int i = r1+1; i < rN; i++ )
	{
		int index = i-r1;
		
		if( residList.find(i-1) == residList.end() ||
			residList.find(i) == residList.end() || 
			residList.find(i+1) == residList.end()) continue;


		shift = inPDB.getPhi(1,i+1); 
		U_ANGLES[index][7] = shift;
		shift = inPDB.getPsi(1,i+1); 
		U_ANGLES[index][8] = shift;
		shift = inPDB.getChi1(1,i+1); 
		U_ANGLES[index][9] = shift;

		U_NAME[index][3] = residList[i+1];

		inTab.Entries[index]["PHI_R3"] = ftoa(U_ANGLES[index][7], buf);
		inTab.Entries[index]["PSI_R3"] = ftoa(U_ANGLES[index][8], buf);
		inTab.Entries[index]["CHI1_R3"] = ftoa(U_ANGLES[index][9], buf);
		inTab.Entries[index]["RESID_R3"] = itoa(i+1, buf);
		inTab.Entries[index]["RESNAME_R3"] = residList[i+1];
		inTab.Entries[index]["SOURCE"] = inName.substr(pos0,pos1-pos0);

		//Ring current shifts
		for(itN = aN.begin(); itN != aN.end(); itN++)
		{
			string name = itN->second;
			if( name == "H" ) 
			{ 
				name = "HN";
				if( residList[i] == "P" ) continue;
			}
			else if( name == "HA" && residList[i] == "G" ) 
			{
				U_RING_SHIFTS[index][7] = inPDB.getOrbitalShift(1,i,"HA1") ;
				name = "HA2";
			}
			else if( name == "CB" && residList[i] == "G" ) continue;
			
			U_RING_SHIFTS[index][itN->first] = inPDB.getOrbitalShift(1,i,name) ;
			inTab.Entries[index][itN->second+"_HM"] = ftoa(U_RING_SHIFTS[index][itN->first], buf);
		}

		//H-Honds
		dist = inPDB.HBDistList[i]["HN"];
		inTab.Entries[index]["H_HB"] = ftoa(dist, buf);
		U_HN_HB[index] = dist;

		dist = inPDB.HBDistList[i]["HA"];
		inTab.Entries[index]["HA_HB"] = ftoa(dist, buf);
		U_HA_HB[index] = dist;

		dist = inPDB.HBDistList[i]["O"];
		inTab.Entries[index]["CO_HB"] = ftoa(dist, buf);
		U_CO_HB[index] = dist;	

		if( inTab.Entries.find(index-1) != inTab.Entries.end() ) { 
			//if tripet i-1 exist
			//assign the values of positions 1 and 2 of triplet i using the postions 2 and 3 of triplet i-1
			U_ANGLES[index][1] = U_ANGLES[index-1][4];
			U_ANGLES[index][2] = U_ANGLES[index-1][5];
			U_ANGLES[index][3] = U_ANGLES[index-1][6];
			U_ANGLES[index][4] = U_ANGLES[index-1][7];
			U_ANGLES[index][5] = U_ANGLES[index-1][8];
			U_ANGLES[index][6] = U_ANGLES[index-1][9];
			U_NAME[index][1] = U_NAME[index-1][2];
			U_NAME[index][2] = U_NAME[index-1][3];

			inTab.Entries[index]["PHI_R1"] = inTab.Entries[index-1]["PHI_R2"];
			inTab.Entries[index]["PSI_R1"] = inTab.Entries[index-1]["PSI_R2"];
			inTab.Entries[index]["CHI1_R1"] = inTab.Entries[index-1]["CHI1_R2"];
			inTab.Entries[index]["RESID_R1"] = inTab.Entries[index-1]["RESID_R2"];
			inTab.Entries[index]["RESNAME_R1"] = inTab.Entries[index-1]["RESNAME_R2"];

			inTab.Entries[index]["PHI_R2"] = inTab.Entries[index-1]["PHI_R3"];
			inTab.Entries[index]["PSI_R2"] = inTab.Entries[index-1]["PSI_R3"];
			inTab.Entries[index]["CHI1_R2"] = inTab.Entries[index-1]["CHI1_R3"];
			inTab.Entries[index]["RESID_R2"] = inTab.Entries[index-1]["RESID_R3"];
			inTab.Entries[index]["RESNAME_R2"] = inTab.Entries[index-1]["RESNAME_R3"];
		}
		else {	//else, calculate the values from coordinates	
			
			shift = inPDB.getPhi(1,i-1); 
			U_ANGLES[index][1] = shift;
			shift = inPDB.getPsi(1,i-1); 
			U_ANGLES[index][2] = shift;
			shift = inPDB.getChi1(1,i-1); 
			U_ANGLES[index][3] = shift;
			U_NAME[index][1] = residList[i-1];

			inTab.setEntry(index, "PHI_R1", ftoa(U_ANGLES[index][1], buf) );
			inTab.setEntry(index, "PSI_R1", ftoa(U_ANGLES[index][2], buf) );
			inTab.setEntry(index, "CHI1_R1", ftoa(U_ANGLES[index][3], buf) );
			inTab.setEntry(index, "RESID_R1", itoa(i-1, buf) );
			inTab.setEntry(index, "RESNAME_R1", residList[i-1] );

			shift = inPDB.getPhi(1,i); 
			U_ANGLES[index][4] = shift;
			shift = inPDB.getPsi(1,i); 
			U_ANGLES[index][5] = shift;
			shift = inPDB.getChi1(1,i); 
			U_ANGLES[index][6] = shift;
			U_NAME[index][2] = residList[i];

			inTab.Entries[index]["PHI_R2"] = ftoa(U_ANGLES[index][4], buf);
			inTab.Entries[index]["PSI_R2"] = ftoa(U_ANGLES[index][5], buf);
			inTab.Entries[index]["CHI1_R2"] = ftoa(U_ANGLES[index][6], buf);
			inTab.Entries[index]["RESID_R2"] = itoa(i, buf);
			inTab.Entries[index]["RESNAME_R2"] = residList[i];
		}
	}



	//calculate H-bond for the last residue
	dist = inPDB.HBDistList[rN]["HN"];
	U_HN_HB[rN-r1] = dist;
	dist = inPDB.HBDistList[rN]["HA"];
	U_HA_HB[rN-r1] = dist;
	dist = inPDB.HBDistList[rN]["O"];
	U_CO_HB[rN-r1] = dist;	
	
	inTab.saveGDB(PRED_DIR+slash_char+inName.substr(pos0,pos1-pos0) + "_in.tab");
}



void SPARTA::Compile()
{
	
	GDB pdbList;
	
	if( pdbListName.length() > 0) {
		cout << "Reading Candidate Proteins Names from " << pdbListName ;
		pdbList.loadGDB( pdbListName );
		cout << ", " << pdbList.Entries.size() << " proteins read." << endl;
	}
	else {
		cout << "PDB name list is required " << endl ;
		cout << "\t sparta -compile -pdbList list.tab " << endl ;
		
		exit(0);
	}
	
	
	cout << "Reading Random Coil Shifts from " << rcFileName << endl;
	RC_Tab.loadGDB( rcFileName );

	cout << "Reading RC Adjustments from " << adjFileName << endl;
	ADJ_Tab.loadGDB( adjFileName );

	cout << "Reading Previous Residue RC Adjustments from " << prevFileName << endl;
	PREV_Tab.loadGDB( prevFileName );

	cout << "Reading Next Residue RC Adjustments from " << nextFileName << endl;
	NEXT_Tab.loadGDB( nextFileName );


	inTab.VARS_str_parser(" RESID_R1 RESNAME_R1 PHI_R1 PSI_R1 CHI1_R1 RESID_R2 RESNAME_R2 PHI_R2 PSI_R2 CHI1_R2 RESID_R3 RESNAME_R3 PHI_R3 PSI_R3 CHI1_R3 N HA C CA CB H SOURCE ");
	inTab.FORMAT_str_parser(" %4d %s %8.3f %8.3f %8.3f %4d %s %8.3f %8.3f %8.3f %4d %s %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s ");
	inTab.Entries.clear();


	map< int, map<string, string> >::iterator it;

	cout << "Checking files ...." << endl; // checking if files are exsit
	bool checked = true;
	for(it = pdbList.Entries.begin(); it != pdbList.Entries.end(); it++)
	{
		string pname = it->second["PDB_NAME"];
		string name = pname.substr(0,pname.find_last_of("."));
		string pdbName = PDB_DIR + slash_char + name + ".pdb";
		string shiftName = SHIFT_DIR + slash_char + name + ".tab";

		ifstream file(shiftName.c_str());
		if(! file.is_open() ){
			cerr << "\tshifts file " << shiftName << " not found." << endl;
			checked = false;
		}
		file.close();
		
		ifstream file2(pdbName.c_str());
		if(! file2.is_open() ){
			cerr << "\tpdb file " << pdbName << " not found." << endl;
			checked = false;		
		}
		file2.close();
	}
	if(!checked) exit(0);
	
	//compile each protein
	for(it = pdbList.Entries.begin(); it != pdbList.Entries.end(); it++)
	{
		string name = it->second["PDB_NAME"];
		//cout << "\tCompiling protein " << name << "..." << endl;
		Compile( name.substr(0,name.find_last_of(".")) );
	}
	

	string str = tripFileName.substr(0, tripFileName.find_last_of('.')) + "_new.tab";
	cout << "Writing new triplet database " << str << endl;
	inTab.saveGDB(str);

}



void SPARTA::Compile(const string &name)
{
	string pdbName = PDB_DIR + slash_char + name + ".pdb";
	string shiftName = SHIFT_DIR + slash_char + name + ".tab";

	GDB shift(shiftName);
	PDB pdb(pdbName);

	map< int, map<string, string> >::iterator it;
	map<int, string>::iterator itN;
	map<string, float>::iterator itS;

	map<int, string> resList;
	map<int, map<string, float> > shifts; //indexed by resID, atomName and shifts
	int r1 = 9999, rN =-9999;

	// pre-read the shift table
	for(it = shift.Entries.begin(); it != shift.Entries.end(); it++)
	{
		string aName = it->second["ATOMNAME"];
		string resName = it->second["RESNAME"];
		if( aName == "HA3" || aName == "HA1") continue;

		for(itN = aN_ALL.begin(); itN != aN_ALL.end(); itN++)
			if( aName == itN->second ) break;

		if( itN == aN_ALL.end() && !(aName == "HA2" && resName == "G") ) continue; //skip entry with invalid atom name

		int resID = atoi( it->second["RESID"].c_str() );
		if( r1 > resID ) r1 = resID;
		if( rN < resID ) rN = resID;

		resList[resID] = resName;
		
		float shf = atof( it->second["SHIFT"].c_str() );

		if( aName == "HA2" && resName == "G" ) //for GLY, using the average HA2, HA3 shifts
		{
			GDB_Entry HA3 = shift.getEntry( "RESID", it->second["RESID"], "ATOMNAME", "HA3", 1 ); //HA1 or HA3
			float shift_HA3= atof( HA3["SHIFT"].c_str() );

			shf = 0.5*( shift_HA3 + shf );
		}
		
		if( fabs(shf) < 0.001 )  shf = 0.001;
		
		shifts[resID][aName] = shf;
	}
	cout << "Compiling protein " << name << ", residue " << r1 << " to " << rN << " ... " << endl;

	pdb.initOrbitalShift();

	map<int, map<string, float> > angles;
	// pre-read the pdb coordinates and calculate the angles for each residue
	for(int i = r1; i <= rN; i++)
	{
		if( !(resList[i].empty()) && !(pdb.residListOne[i].empty()) 
			&& resList[i] != pdb.residListOne[i] )
		{
			if( resList[i] == "C" || resList[i] == "c" ) 
			{
				cerr << "\tWarning: Incompatible CYS oxidation status for residue " << i << endl; 
				resList[i] = pdb.residListOne[i];
			}
			else 
				cerr << "\tWarning: Mutated residues " << i << "(" << resList[i] << "->" << pdb.residListOne[i] << ")" << endl; 
		}
		angles[i]["PHI"] = pdb.getPhi(1,i);
		angles[i]["PSI"] = pdb.getPsi(1,i);
		angles[i]["CHI1"] = pdb.getChi1(1,i);

	}


	for(int i = r1+1; i < rN; i++) //loop over from 2nd residue to the 2nd to the last residue
	{
		if( shifts[i-1].size() <= 2 || shifts[i].size() <= 2 || shifts[i+1].size() <= 2 ) continue;
		//skip resdiue with less than two chemical shifts assignments 
		if( angles[i]["PHI"] > 999.0 && angles[i]["PSI"] > 999.0 ) continue;
		if( angles[i-1]["PHI"] > 999.0 && angles[i-1]["PSI"] > 999.0 ) continue;
		if( angles[i+1]["PHI"] > 999.0 && angles[i+1]["PSI"] > 999.0 ) continue;
		//skip resdiue with no coordinates in the range of triplet
		
		int index = inTab.Entries.size()+1;

		inTab.Entries[index]["PHI_R1"] = ftoa(angles[i-1]["PHI"], buf);
		inTab.Entries[index]["PSI_R1"] = ftoa(angles[i-1]["PSI"], buf);
		inTab.Entries[index]["CHI1_R1"] = ftoa(angles[i-1]["CHI1"], buf);
		inTab.Entries[index]["PHI_R2"] = ftoa(angles[i]["PHI"], buf);
		inTab.Entries[index]["PSI_R2"] = ftoa(angles[i]["PSI"], buf);
		inTab.Entries[index]["CHI1_R2"] = ftoa(angles[i]["CHI1"], buf);
		inTab.Entries[index]["PHI_R3"] = ftoa(angles[i+1]["PHI"], buf);
		inTab.Entries[index]["PSI_R3"] = ftoa(angles[i+1]["PSI"], buf);
		inTab.Entries[index]["CHI1_R3"] = ftoa(angles[i+1]["CHI1"], buf);

		inTab.Entries[index]["RESID_R1"] = itoa(i-1, buf);
		inTab.Entries[index]["RESID_R2"] = itoa(i, buf);
		inTab.Entries[index]["RESID_R3"] = itoa(i+1, buf);

		inTab.Entries[index]["RESNAME_R1"] = resList[i-1];
		inTab.Entries[index]["RESNAME_R2"] = resList[i];
		inTab.Entries[index]["RESNAME_R3"] = resList[i+1];

		//shifts
		for(itN = aN_ALL.begin(); itN != aN_ALL.end(); itN++)
		{
			string aName = itN->second;
			float shft = shifts[i][aName];

			if( aName == "HA" && resList[i] == "G" ) 
				shft = shifts[i]["HA2"];
			
			if( shft == 0.0 ) 
			{
				inTab.Entries[index][aName] = ftoa(MAX, buf);
				continue;
			}

			float RC = getRC(resList[i],aName);
			float RCadj = getRCadj(resList[i],aName);
			float PrevRCadj = getPrevRCadj(resList[i-1],aName);
			float NextRCadj = getNextRCadj(resList[i+1],aName);
			//ring current shifts
			if( aName == "H" ) aName = "HN";
			else if( aName == "HA" && resList[i] == "G" ) aName = "HA2";

			float HM_Ring = pdb.getOrbitalShift(1,i,aName);

			if( aName == "HA" && resList[i] == "G" ) HM_Ring = 0.5*(HM_Ring + pdb.getOrbitalShift(1,i,"HA1"));
			if( fabs(HM_Ring) > 10.0 ) HM_Ring = 0;
			
			shft -= ( RC + RCadj + PrevRCadj + NextRCadj + 0.6*HM_Ring);
			if( shft >= 999.0 || RC + RCadj + PrevRCadj + NextRCadj + 0.6*HM_Ring >= 999.0 ) shft = MAX;
			inTab.Entries[index][itN->second] = ftoa(shft, buf);
		}

		inTab.Entries[index]["SOURCE"] = name;
	}

}



// get random coil chemical shift for atom 'aName' of residue 'resName'
float SPARTA::getRC(const string& resName, const string& aName)
{
	GDB_Entry temp = RC_Tab.getEntry("RESNAME",resName,1);

	if(temp.size() != 0) return atof( temp[aName].c_str() );
	
	return 9999.0;
}



float SPARTA::getRCadj(const string& resName, const string& aName)
{
	GDB_Entry temp = ADJ_Tab.getEntry("RESNAME",resName,1);

	if(temp.size() != 0) return atof( temp[aName].c_str() );
	
	return 0.0;
}



float SPARTA::getPrevRCadj(const string& prev_rName, const string& aName)
{
	GDB_Entry temp = PREV_Tab.getEntry("RESNAME",prev_rName,1);

	if(temp.size() != 0) return atof( temp[aName].c_str() );
	
	return 0.0;
}



float SPARTA::getNextRCadj(const string& next_rName, const string& aName)
{
	GDB_Entry temp = NEXT_Tab.getEntry("RESNAME",next_rName,1);

	if(temp.size() != 0) return atof( temp[aName].c_str() );
	
	return 0.0;
}



float SPARTA::getWeight(const string& Name, const string& aName)
{
	GDB_Entry temp = WEIGHT_Tab.getEntry("RESNAME",Name,1);

	if(temp.size() != 0) return atof( temp[aName].c_str() );
	
	return 9999.0;

}



void SPARTA::mkdir_pred(const string& d)// create a directory for prediction results
{
	if( ! isDirExists( d ) )
	{
		int i = MKDIR(d.c_str())+1;
		if( !i ) { // if not success
			string parentD = d.substr(0,d.find_last_of(slash_char)+1);
			mkdir_pred(parentD.c_str());
			if ( !(MKDIR(d.c_str())+1) ) {
				cerr << "\tCan't create prediction directory " << d.c_str() << endl; 
				exit(0);
			}
		}
	}
}



float SPARTA::getDiff( float ang1, float ang2)
{
	float a = fabs(ang1-ang2);
	
	if( a < 180 ) return a;
	else if( a > 180 && a < 360.1 && ang1*ang2 < 0)
	{	
		return 360.0-a;
	}

	return 0;
}



char * SPARTA::itoa( int n, char *buff, int base )
{
	sprintf(buff, "%d", n);
	return buff;
}



char * SPARTA::ftoa( float n, char *buff, char f, int prec )
{
    if ( !(f=='f' || f=='F' || f=='e' || f=='E' || f=='g' || f=='G') ) {
		f = 'f';
    }
    char format[20];
    char *fs = format;				// generate format string
    *fs++ = '%';					//   "%.<prec>l<f>"
    if ( prec >= 0 ) {
		if ( prec > 99 )			// buf big enough for precision?
			prec = 99;
		*fs++ = '.';
		if ( prec >= 10 ) {
			*fs++ = prec / 10 + '0';
			*fs++ = prec % 10 + '0';
		} else {
		    *fs++ = prec + '0';
		}
    }
    *fs++ = 'l';
    *fs++ = f;
    *fs = '\0';
    sprintf( buff, format, n );

    return buff;
}



int SPARTA::MKDIR(const char *dirName)
{
#ifdef WIN32
  return mkdir(dirName);
#else
  return mkdir(dirName, 0777);
#endif
}



bool SPARTA::isDirExists(const string &Dir) {
#ifdef WIN32
	char oldDir[100];
	if( NULL == getcwd(oldDir, _MAX_PATH) ) return false;

	if( chdir( Dir.c_str() ) == -1 ) return false;

	if( chdir( oldDir ) == -1 ) return false;

	return true;
#else 
       DIR *dp;
       dp = opendir(Dir.c_str());
       if (dp != NULL)
         {
           (void) closedir (dp);
		   return true;
         }
       else
          return false;
#endif
}



float SPARTA::getAVG(vector<float> &v1)
{
	int n = v1.size();

	float sum = 0.0;

	for(int i = 0; i < n; i++)
		sum += v1[i];

	return sum/(float)n;
}



float SPARTA::getSTD(vector<float> &v1)
{
	int n = v1.size();
	float avg = getAVG(v1);

	float sum = 0.0;

	for(int i = 0; i < n; i++)
		sum += (v1[i]-avg)*(v1[i]-avg);

	return sqrt(sum/(float)(n-1));
}


float SPARTA::getRMS(vector<float> &v1, vector<float> &v2)
{
	if(v1.size() != v2.size() ) return -1.0;

	int n = v1.size();

	float sum = 0.0;
	for(int i = 0; i < n; i++)
	{
		sum += ( (v1[i]-v2[i])*(v1[i]-v2[i]) );
	}

	return sqrt(sum/(float)(n-1));
}

