/*************************************************************** 
	The GPDD Sensitivity Project (GPDD_Sens) 
 
	GPDD stands for "Graph-Pair Decision Diagram" 
 
Contributors: 
 
   CHEN, Weiwei		(MS 2007 graduate) 
   MENG, Xiaoxuan		(MS 2009 graduate) 
   ZENG, Mei			(MS 2010 graduate) 
   LI, Ji				(MS 2010 graduate) 
   MA, Diming			(MS 2011 graduate) 
 
   SHI, Guoyong 		(Prof)  shiguoyong@ic.sjtu.edu.cn 
 
   School of Microelectronics 
   Shanghai Jiao Tong University 
   (c) 2005 - 2010 
******************************************************************/ 

/**
  Analysis_AC.cpp
  08-29-2010
*/

#include  <stdlib.h>
#include	<iostream>
#include	<fstream>
#include	<iomanip>
#include	<string>
#include	<cstdlib>
#include	<iomanip>
#include  <stack>

#include "option.h"
#include	"Analysis_AC.h"

#define		USE_HASHER	1

const char *AC_StepTypeName[] = {"NO_TYPE", "DEC", "OCT", "LIN"};

Analysis_AC::Analysis_AC(GPDD *gpdd)
{
	pGPDD = gpdd;
	
	ACstartFreq = 0.0;
	ACstopFreq = 0.0;
	ACfreqDelta = 0.0; 

	ACstepType = NO_TYPE;    
	// to be set to [DECADE | OCTAVE | LINEAR].
	ACnumberSteps = 0;	
	// to be set to [ND|NO|NP] in the ".AC" command
	
	RelTol = 1e-3;		// Relative Tolerance

	// The TF curve is obtained from the following two vectors:
	NumFreqPoints 	= 0;		// Num of freq points.
	FreqVec 		= NULL;		// A list of freq points.
	RespVec 		= NULL;		// A list of response points.
	MagVec 		= NULL;
	PhaseVec 	= NULL;
	SensVec 		= NULL;

	CurrTraversalMark = false;

	cout << LINE_INFO << "AC analysis constructed." << endl;
	
}


Analysis_AC::~Analysis_AC()
{
	if (FreqVec) {
		delete [] FreqVec;
	}
		
	if (RespVec) {
		delete [] RespVec;
	}
		
	if (MagVec) {
		delete [] MagVec;
	}
		
	if (PhaseVec) {
		delete [] PhaseVec;
	}
		
	if (SensVec) {
		delete [] SensVec;
	}

	cout << LINE_INFO << "Deconstructed." << endl;
	
}


const char*	Analysis_AC::GetACstepTypeName()   
{ 
	return  AC_StepTypeName[ACstepType]; 
}


// InitAnalysis()
// Initializes analysis; check the frequency range specified and 
// prepare for the vector of frequency points.
// Returns OK if no problem.
int  Analysis_AC::InitAnalysis()
{
	int 	i, k;
	double freq;
	double freqTol; /* tolerence parameter for finding the final frequency */

#ifdef	TRACE_
	cout << LINE_INFO << endl;
	cout << "\n\tACstepType is : " << AC_StepTypeName[ACstepType] << endl;
#endif

	/** At least one freq point for AC analysis. */
	if (ACnumberSteps < 1)	   
		ACnumberSteps = 1;

	// Next we deal with the cases of DECADE and OCTAVE.
	// ---- Standard SPICE code -----------------

	switch (ACstepType) 
	{
		case DECADE:
			// 10^(1/ACnumberSteps) : the multiplication factor
			// .AC DEC ND FSTART FSTOP
			// ND = # freq points / decade
			ACfreqDelta = exp(log(10.0)/ACnumberSteps);
			break;
		case OCTAVE:
			// 2^(1/ACnumberSteps)  : the multiplication factor
			ACfreqDelta = exp(log(2.0)/ACnumberSteps);
			break;
		case LINEAR:
			///cout << LINE_INFO << "[ACstepType == LINEAR]" << endl;
			if (ACnumberSteps == 1) {
				ACfreqDelta = 0;
			} else  {  // NumFreqPoints >= 2
				ACfreqDelta = (ACstopFreq - ACstartFreq)/(ACnumberSteps - 1);
				// Do this way to avoid roundoff error for the last freq point.
			}
			break;
		default:
			cout << LINE_INFO << "[ACstepType] not recognized.";
			return(STOP);
	}

	// Fix <NumFreqPoints>
	switch(ACstepType) 
	{
		case DECADE:
		case OCTAVE:
			freqTol = ACfreqDelta * ACstopFreq * RelTol;

			// First count the number of frequency points requested
			// so that we can allocate memory for FreqVec[] and RespVec[].
			// It looks like brute force, but ... (why not) ...
			
			freq = ACstartFreq;
			k = 0;	// counting the number of freq points

			///cout << LINE_INFO << "freqTol = " << freqTol << endl;
			
			while (freq <= ACstopFreq+freqTol) 
			{
				if (ACfreqDelta == 1) 	
					break;
				freq *= ACfreqDelta;
				++k;
			}

			NumFreqPoints = k;		// Number of freq points obtained.
			break;
			
		case LINEAR:
			///cout << LINE_INFO << "[ACstepType == LINEAR]" << endl;
			NumFreqPoints = ACnumberSteps;
			break;
	}

	
	// Allocate storage memory ...
	FreqVec 	= new double[NumFreqPoints];
	RespVec 	= new complex_t[NumFreqPoints];
	MagVec 	= new double[NumFreqPoints];
	PhaseVec = new double[NumFreqPoints];
	SensVec 	= new complex_t[NumFreqPoints];

	// Now fill in the frequency points in FreqVec[].
	FreqVec[0] = ACstartFreq;
	freq = ACstartFreq;

	switch (ACstepType) 
	{
		case DECADE:
		case OCTAVE:

			for (k = 1; k < NumFreqPoints; k++) {
				freq *= ACfreqDelta;
				FreqVec[k] = freq;
			}
			break;

		case LINEAR:
			///cout << LINE_INFO << "[ACstepType == LINEAR]" << endl;
			NumFreqPoints = ACnumberSteps;

			//FreqVec[0] = ACstartFreq;
			for (i = 1; i < NumFreqPoints-1; i++) {
				FreqVec[i] = FreqVec[i-1] + ACfreqDelta;
			}
			FreqVec[NumFreqPoints-1] = ACstopFreq;
			// Do this way to avoid roundoff error for the last freq point.
			break;
	}

	cout << "\n\tAnalysis_AC::InitAnalysis() -- done";

	return (OK);
} 


/**
// Numerical frequency response are computed and managed
// (saving/plotting/outputing) here.
*/
int Analysis_AC::RunACAnalysis()
{
	int 			error;
	double 		freq;
	GPDDnode 	*pRoot;

//#ifdef	TRACE
	cout << LINE_INFO << endl;
//#endif

	if (!pGPDD)
	{
		cout << LINE_INFO << "No GPDD for analysis." << endl;
		return (NOPE);
	}

	if (InitAnalysis() ==  STOP) {
		cout << LINE_INFO << "InitAnalysis() failed." << endl;
		return (STOP);
	}
	

#ifdef	S_EXPANDED_ANALYSIS
	cout << LINE_INFO << "Evaluation by s-expansion." << endl;
	EvaluateCoefficients();	// Evaluate the s-coefficients a_n and b_n only once.
#else
	cout << LINE_INFO << "Evaluation without s-expansion." << endl;
#endif

	// Run through all the requested freq points
	for (int i = 0; i < NumFreqPoints; i++) 
	{
		freq = FreqVec[i];

#ifdef	S_EXPANDED_ANALYSIS
		RespVec[i] = ExpandedEvaluationAtFreq(freq);
#else
		RespVec[i] = EvaluateAtFreq(freq);
#endif
		MagVec[i] = magn(RespVec[i]);
		PhaseVec[i] = arg(RespVec[i]);
	}

	PrintFreqResponse(NULL);
	
	 // Dump or save data for plotting etc. HERE ...
	
//#ifdef	TRACE
	cout << LINE_INFO << "-- done" << endl;
//#endif
	
	return(OK);
}


/**
	Caution: Before calling this function, the s-expanded coefficients must have been
	evaluated once. (G. Shi)

	It is suggeted to use this function for difference sensitivity analysis.

	The memory needed for <resp> must be allocated a priori.
*/
void  Analysis_AC::EvaluateResponseVector(complex_t *resp)
{
	cout << LINE_INFO << endl;
	double 		freq;

	// Run through all the requested freq points
	for (int i = 0; i < NumFreqPoints; i++) 
	{
		freq = FreqVec[i];

		resp[i] = ExpandedEvaluationAtFreq(freq);
		
///		RespVec[i] = ExpandedEvaluationAtFreq(freq);
///		RespVec[i] = EvaluateAtFreq(freq);

//		MagVec[i] = magn(RespVec[i]);
//		PhaseVec[i] = arg(RespVec[i]);
	}
}


/**
  	Calculate sensitivity by difference approximation.
  	DiffSens(H, p) = (p/H) (H(p+dp) - H(p)) / (dp).

  	Hence, it requires a parameter variation (dp) and two evaluations of
  	H(s, p1) and H(s, p2). That is, to evaluate the s-polynomial coefficients twice.
  	(G. Shi, 09-05-2010)

  	<percent> is the percentage of variation.
*/
void  Analysis_AC::DifferenceSensitivity(const string &paraname, double percent)
{
	cout << LINE_INFO << endl;

	int sens_index = pGPDD->LookupSymbolIndex(paraname.c_str());
	
	if (sens_index == NOPE) {
		cout << LINE_INFO << "\n\tCould not find the symbol \"" 
			<< paraname << "\"" << endl;
		exit(0);
	}

	if (NumFreqPoints == 0) {
		cout << LINE_INFO << "Nothing has to be done." << endl;
		return;
	}
	
	complex_t	*respvec1 = new complex_t[NumFreqPoints];
	complex_t	*respvec2 = new complex_t[NumFreqPoints];
//	complex_t	*sensvec = new complex_t[NumFreqPoints];
	
	double time_1 = clock();
	
	cout<< LINE_INFO << "\n\tFirst evaluation of coefficients." << endl;
	EvaluateCoefficients();   // evaluate all coefficients of s-expanded H(s)
	EvaluateResponseVector(respvec1);

//	PrintFreqResponse2(respvec1);
	
	double time_2 = clock();
	
	cout << LINE_INFO << "\n\tTime for the 1st evaluation is: " << (time_2-time_1)/1000000 << endl;

	double 	old_value = pGPDD->GetSymbolValue(sens_index);
	// Symbol value to be restored
		
	cout<< LINE_INFO << "\n\tChange one parameter value." << endl;
	if (pGPDD->VariationByPercentage(paraname, percent) == STOP) {
		cout<< LINE_INFO << "Program terminated." << endl;
		exit(0);
	}
	
	pGPDD->PrintSymbols(); 
	
	time_1 = clock();
	
	cout<< LINE_INFO << "\n\tSecond evaluation of coefficients." << endl;
	EvaluateCoefficients(); // evaluate all coefficients of s-expanded H(s)
	EvaluateResponseVector(respvec2);

//	PrintFreqResponse2(respvec2);
	
	time_2 = clock();
	
	cout<< LINE_INFO << "\n\tTime for the 2nd evaluation is: " << (time_2-time_1)/1000000 << endl;

	cout<< LINE_INFO << "\n\tCalculate the difference sensitivity." << endl;

	for (int i = 0; i < NumFreqPoints; i++) {
		SensVec[i] = (respvec2[i] - respvec1[i]) / (complex_t(percent) * respvec1[i]);
	}

	cout<< LINE_INFO << "\n\tPrint the \"difference\" sensitivity." << endl;

	PrintSensResponse(NULL);

	pGPDD->SetSymbolValue(sens_index, old_value);

	delete [] respvec1;
	delete [] respvec2;
	
}

/**
  Direct evaluate the GPDD at a given frequency without s-expansion.
*/
complex_t Analysis_AC::EvaluateAtFreq(double freq)
{
#ifdef	TRACE
	cout << LINE_INFO << endl;
#endif

	GPDDnode *pRoot = pGPDD->GetGPDDroot();

	GPDDnode *num = pRoot->GetVertex_IN();
	GPDDnode *den = pRoot->GetVertex_OUT();

	complex_t  val_den, val_num;

	// Set the composite symbol values.
	// This only can be done after a frequency is selected.
	
	pGPDD->EvaluateAllSymbolsAtFreq(freq);

#if 0
	HasherTree  *ht = new HasherTree();

	val_num = EvaluateH(num, freq, ht);
	val_den = EvaluateH(den, freq, ht);
	
	delete  (HasherTree *) ht;
	//delete sometimes causes problem.
#else

///	pGPDD->UnMarkAllNodes(pGPDD->GetGPDDroot(), 1);

	CurrTraversalMark = ! pRoot->GetMark();  // Flip the existing mark
	
	val_num = Evaluate_R(num, freq);
	val_den = Evaluate_R(den, freq);
	
	pRoot->SetMark(CurrTraversalMark);
	
#endif

	if (val_den == complex_t(0))
	{
		cout << LINE_INFO << "Divide by zero.";
		return  -val_num * DBL_MAX;
	}
	else
		return  -val_num * pRoot->GetSign_IN() / val_den * pRoot->GetSign_OUT();
}



Complex_H  Analysis_AC::EvaluateFractionAtFreq(double freq)
{
	GPDDnode *pRoot = pGPDD->GetGPDDroot();
	
	pRoot->Mark();

	complex_t val_den, val_num;

	pGPDD->EvaluateAllSymbolsAtFreq(freq);

//	pGPDD->UnReferenceAll(pGPDD->GetGPDDroot(), 1);
//	pGPDD->UnMarkAllNodes(pGPDD->GetGPDDroot(), 1);
	
	GPDDnode *num = pRoot->GetVertex_IN();
	GPDDnode *den = pRoot->GetVertex_OUT();

	CurrTraversalMark = ! pRoot->GetMark();
///	pGPDD->SetTraversalMark(CurrTraversalMark);
	
	val_num = Evaluate_R(num, freq);
	val_den = Evaluate_R(den, freq);
	
	pRoot->SetMark(CurrTraversalMark);

	Complex_H 	fraction;
	
	fraction.num = val_num;
	fraction.den = val_den;
	//cout << "val_num is " << val_num << endl;	
	
	return fraction;	
}


	
complex_t	Analysis_AC::Evaluate_R(GPDDnode *P, double freq)
{
	if ( P->GetSymbolIndex() == -2 )		// 'zero' node
	{
		return complex_t(0.0);
	}
	else if (P->GetSymbolIndex() == -1)	// 'one' node
	{
		return complex_t(1.0);
	}

	if (P->GetMark() == CurrTraversalMark)
		return P->GetValue();
	else
		P->SetMark(CurrTraversalMark);

	complex_t	V, V1, V0;


	double	value = pGPDD->GetSymbolValue(P->GetSymbolIndex());
	double	imag;

	switch ( pGPDD->GetSymbolName(P->GetSymbolIndex())[0] )
	{
	case 'R':
		if (value == 0)
			V = complex_t(DBL_MAX);
		else
			V = complex_t(1/value);
		//	V = 1.0;
		break;
	case 'L':
		if (freq == 0)
			V = complex_t(0, -DBL_MAX);
		else
			V = complex_t(0, -1/(2*PI*freq*value));
		//	V = 1.0;
		break;
	case 'C':
		//V = complex_t(0, 1.0);
		V = complex_t(0, 2*PI*freq*value);
		break;
	case '^':	// Composite symbol
		imag = pGPDD->GetSymbolValue_i(P->GetSymbolIndex());
		V = complex_t(value, imag);
		//cerr << "\n  V = " << V;
		break;
	default:		// Gain, or +/-: "+" = 1, "-" = -1.
		V = complex_t(value);
		break;
	}

	V1 = Evaluate_R(P->GetVertex_IN(), freq);
	V0 = Evaluate_R(P->GetVertex_OUT(), freq);
	
	V  = V * V1 * P->GetSign_IN() + V0 * P->GetSign_OUT();

	if ( P->GetValueP() ) 	// to avoid repeated memory creation
	{
		P->SetValue(V);
	} 
	else
	{	// Memory is allocated for only once.
		complex_t *val = new complex_t(V);
		P->SetValueP(val); // allocate memory for the node value
	}

	return V;
}

/**
	Evaluate with the help of a hash tree.
*/
complex_t  Analysis_AC::EvaluateH(GPDDnode *P, double freq, HashTree *ht)
{
	if ( P->GetSymbolIndex() == -2 )		// 'zero' node
		return complex_t(0.0);
	else if (P->GetSymbolIndex() == -1)	// 'one' node
		return complex_t(1.0);

	if ( !ht->Put(P, *P) )	// Node value already computed.
	{
		/*
		cout << "\nHit hasher.";
		cout << "\nValue = " << P->GetValue();
		*/
		return P->GetValue();
	}

	complex_t	V, V1, V0;
	// Get the symbol value.
	double	value = pGPDD->GetSymbolValue(P->GetSymbolIndex());
	double 	imag;
	
//cout << "\nvalue = " << value;

//cout << "\nFile: " << __FILE__ << ", Line: " << __LINE__;
	// Get the element value: R/L/C/Gain.
	switch ( pGPDD->GetSymbolName(P->GetSymbolIndex())[0] )
	{
	case 'R':
		if (value == 0)
			V = complex_t(DBL_MAX);
		else
			V = complex_t(1/value);
		break;
	case 'L':
		if (freq == 0)
			V = complex_t(0, -DBL_MAX);
		else
			V = complex_t(0, -1/(2*PI*freq*value));
		break;
	case 'C':
		V = complex_t(0, 2*PI*freq*value);
		break;
	case '^':	// Composite symbol
		imag = pGPDD->GetSymbolValue_i(P->GetSymbolIndex());
		V = complex_t(value, imag);
		//cerr << "\n  V = " << V;
		break;
	default:		// Gain or "+" or "-".
		V = complex_t(value);
		break;
	}

	V1 = EvaluateH(P->GetVertex_IN(), freq, ht);
	V0 = EvaluateH(P->GetVertex_OUT(), freq, ht);
	
	V  = V * V1 * P->GetSign_IN() + V0 * P->GetSign_OUT();
	
#if DEBUG
	//cerr << "\nP->index = " << P->Index() << ", V = " << V;
#endif

	if ( P->GetValueP() )
		P->SetValue(V);
	else
	{	// Must allocate memory for the first time use.
		complex_t *val = new complex_t(V);
		P->SetValueP(val); // allocate memory for the node value
	}

	return V;
}

/**
	s-expanded transfer function takes the following form
    	<BR><UL>					a0 + a1 * s + a2 * s^2 + ... + am * s^m </UL> 
			TF =    --------------------------------------------------------
	<BR><UL>				    b0 + b1 * s + b2 * s^2 + ... + bn * s^n</UL>

	In this form, the symbol values are evaluated only once for a set of frequency points.

	sCoeffList[ ] keeps the record of all s-coefficients corresponding to the powers of s.
	sPowerList[ ] keeps the record of the powers of s.

	The two vectors are with the two children nodes at the s-expanded GPDD root. 
*/

complex_t Analysis_AC::
ExpandedEvaluationAtFreq(double freq)
{	
	int 	i, j;
	double 		coeff_val;
	complex_t	termval;
	complex_t  	val_den, val_num;
	complex_t	power_of_s;
	complex_t	H;	// Frequency response

#ifdef	TRACE
	cout << LINE_INFO << endl;
#endif

	ExpGPDDnode *pExpRoot = pGPDD->GetExpGPDDroot();

	// Compute the numerator of H(s) = N(s) / D(s)
	vector<sCoeffNode*> 	*psCoeff = pExpRoot->Vertex_IN->sCoeffList;
	vector<int>         		*psPower = pExpRoot->Vertex_IN->sPowerList;

	val_num = complex_t(0.0);

	for(i = 0; i < psCoeff->size(); i ++)
	{
		coeff_val = (*psCoeff)[i]->GetValue(); 	// s-coefficient is supposed to be real.
		
		power_of_s = complex_t(1.0);
		
		/*printf("\n[%s : line %d] \n\tNum[%d] value = (%g, %g)", 
					__FILE__, __LINE__, i, real(coeff_val), imag(coeff_val)); */

		// Calculate (j*\omega)^power
		for (j = 0; j < abs((*psPower)[i]); j ++)	// power could be negative, so take abs( ).
		{
			power_of_s = power_of_s * complex_t(0, 2 * PI * freq);  // (j*2*pi*f)
		}
		
		// Power could be negative for circuits with L (i.e., 1/(L*s)).
		if ((*psPower)[i] < 0)	// b_n * s^n
			termval = complex_t(coeff_val) / power_of_s;
		else		// b_n * s^(-n)
			termval = complex_t(coeff_val) * power_of_s;
		
		val_num += termval;
	}

	// Compute the denominator of H(s) = N(s) / D(s)
	psCoeff	= pExpRoot->Vertex_OUT->sCoeffList;
	psPower	= pExpRoot->Vertex_OUT->sPowerList;

	val_den = complex_t(0.0);

	for (i = 0; i < psCoeff->size(); i ++)
	{
		coeff_val = (*psCoeff)[i]->GetValue();
		power_of_s = complex_t(1.0);
		
		/*printf("\n[%s : line %d] \n\tDen[%d] value = (%g, %g)", 
					__FILE__, __LINE__, i, real(coeff_val), imag(coeff_val)); */

		for (j = 0; j < abs((*psPower)[i]); j ++)
		{
			power_of_s = power_of_s * complex_t(0, 2 * PI * freq);
		}
		
		if ((*psPower)[i] < 0) 
			termval = complex_t(coeff_val) / power_of_s;
		else
			termval = complex_t(coeff_val) * power_of_s;
		
		val_den += termval;
	}

	if (val_den == complex_t(0)) 
	{
		cout << LINE_INFO << "Divide by zero.";
		exit(0);
//		return  -val_num * DBL_MAX;
	} 
	else 
	{
		H =  -val_num * pExpRoot->Sign_IN / val_den * pExpRoot->Sign_OUT;
		/*
		printf("\n[%s : line %d]   val_num = (%g, %g)", __FILE__, __LINE__, real(val_num), imag(val_num));
		printf("\n[%s : line %d]   val_den = (%g, %g)", __FILE__, __LINE__, real(val_den), imag(val_den));
		printf("\n[%s : line %d]   signs = (%d, %d)", __FILE__, __LINE__, 
							pExpRoot->Sign_IN, pExpRoot->Sign_IN);
		printf("\n[%s : line %d]   H = (%g, %g)", __FILE__, __LINE__, real(H), imag(H));
		*/
		return (H);
	}
}


/**	EvaluateCoefficients()
	evaluates the coefficients of the s-expanded transfer function (TF) 
	with the symbols substituted by their numerical values.

	The purpose of s-expanded evaluation is for speedup in frequency response (FR)
	evaluation where a range for s is to be scaned for the same set of symbol values.
	
	The s-expanded TF takes the following form:
	
      				b0 + b1 * s + b2 * s^2 + ... + bm * s^m  
	     H(s) =  ---------------------------------------------
		    		a0 + a1 * s + a2 * s^2 + ... + an * s^n
*/

void Analysis_AC::EvaluateCoefficients()
{
//#ifdef	TRACE
	cout << LINE_INFO << endl;
//#endif

//	if (!pGPDD->ExpGPDDroot) 
	if (!pGPDD->LumpExpanded()) 
	{
		if (pGPDD->ExpandLumpedGPDD() == STOP) {
			cout << LINE_INFO << "Failed to expand lumped GPDD." << endl;
			exit(0);
		}
	}


	ExpGPDDnode *pExpRoot = pGPDD->ExpGPDDroot;

	CurrTraversalMark = ! pExpRoot->GetMark();  // Flip the existing mark
	// pGPDD->FlipTraversalMark();
	
///	pGPDD->UnMarkAllNodes(pExpRoot, 0);
	
//	cout << LINE_INFO << "\n\tEvaluating the Numerator coefficients: " << endl;
	EvaluateCoefficients_R(pExpRoot->Vertex_IN);
	
//	cout << LINE_INFO << "\n\tEvaluating the Denominator coefficients: " << endl;
//	pGPDD->UnMarkAllNodes(pGPDD->ExpGPDDroot, 0);
	EvaluateCoefficients_R(pExpRoot->Vertex_OUT);

	pExpRoot->SetMark(CurrTraversalMark);

//	PrintCoefficients();

}

/**
	Evaluate the expanded GPDD using the values set to the ExpGPDDnode's.
*/
void Analysis_AC::EvaluateCoefficients_R(ExpGPDDnode *P)
{
	if ( P->SymbolIndex == -2 )		// GPDD vertex 'zero' 
	{
		return;
	}
	else if (P->SymbolIndex == -1)	// GPDD vertex 'one' 
	{
		return;
	}

	if (P->GetMark() == CurrTraversalMark) {
		///cout << LINE_INFO << "Return at [#" << P->SymbolIndex << "]" << endl;
		return;
	}
	else {
		P->SetMark(CurrTraversalMark);
		///cout << LINE_INFO << "Marked vertex [#" << P->SymbolIndex << "]" << endl;
	}

	EvaluateCoefficients_R(P->Vertex_IN);
	EvaluateCoefficients_R(P->Vertex_OUT);

	int 	i;
	///complex_t   V, V1, V0;
	double   V, V1, V0;		// suffices to used double evaluation.
	SymbNode  *pSymbN;

	char eType = pGPDD->GetLumpedSymbList()[P->SymbolIndex]->symbname[0];
	
	vector<sCoeffNode*> *psCoeff = P->sCoeffList;

	/*cout << LINE_INFO << "\n\tWorking with vertex [#" << P->SymbolIndex 
		<< "] with #coeff = " 
		<< psCoeff->size() << endl;*/
	
	if (eType == '^')
	{
		eType = pGPDD->GetLumpedSymbList()[P->SymbolIndex]->symbname[1];
		pSymbN = pGPDD->GetLumpedSymbList()[P->SymbolIndex]->nextP;
	}
	else
	{
		pSymbN = pGPDD->GetLumpedSymbList()[P->SymbolIndex];
	}

	double value = pSymbN->value;

	switch ( eType)
	{
	case 'R':
//		cout << LINE_INFO << "\n\tEvaluate [R] : value = " << value << endl;
		if (value == 0)
			V = DBL_MAX;
		else
			V = 1/value;
		break;
	case 'L':
//		cout << LINE_INFO << "Evaluate [L]" << endl;
		if (value == 0)
			V = DBL_MAX;
		else
			V = 1/value;
		break;
	case 'C':
//		cout << LINE_INFO << "Evaluate [C]" << endl;
		V = value;
		break;
	default:		// Gain, or +/-: "+" = 1, "-" = -1.
//		cout << LINE_INFO << "Evaluate others" << endl;
		V = value;
		break;
	}

	double  Vtop;

	for (i = 0; i < psCoeff->size(); i ++)
	{
		V1 = (*psCoeff)[i]->Vertex_IN->GetValue();
		V0 = (*psCoeff)[i]->Vertex_OUT->GetValue();
		
		/*cout << LINE_INFO << "\n\t[i = " << i << "]  The value triple is (" 
			<< V << ", " << V1 << ", " << V0 << ")" << endl; */
		
		Vtop   = (V) * (V1) * (*psCoeff)[i]->pExpNode->Sign_IN + 
					(V0) * (*psCoeff)[i]->pExpNode->Sign_OUT;

		(*psCoeff)[i]->Value = Vtop;

		/*if ( !(*psCoeff)[i]->val_ptr ) {
			(*psCoeff)[i]->val_ptr = new complex_t(V);
			cout << LINE_INFO << "\n\tCreated coeff value = " << V << endl;
		}
		else
			*((*psCoeff)[i]->val_ptr) = V; */
	}

}


int  Analysis_AC::SetParameter(int index, G_Value *value)
{
#ifdef	TRACE_
			cout << LINE_INFO << endl;
#endif

	switch(index) {
	    case AC_DEC:
			ACstepType = DECADE;  
#ifdef	TRACE_
			cout << LINE_INFO << "[ACstepType = DECADE]" << endl;
			printf("\nACstepType = %s\n", AC_StepTypeName[ACstepType]);
#endif
			break;

	    case AC_OCT:
			ACstepType = OCTAVE;  
#ifdef	TRACE_
			cout << LINE_INFO << "[ACstepType = OCTAVE]" << endl;
			printf("\nACstepType = %s\n", AC_StepTypeName[ACstepType]);
#endif
			break;

	    case AC_LIN:
			ACstepType = LINEAR;  
#ifdef	TRACE_
			cout << LINE_INFO << "[ACstepType = LINEAR]" << endl;
			printf("\nACstepType = %s\n", AC_StepTypeName[ACstepType]);
#endif
			break;

	    case AC_STEPS:
			ACnumberSteps = value->iValue;
#ifdef	TRACE_
			cout << LINE_INFO << "ACnumberSteps = " << ACnumberSteps << endl;
#endif
			break;

	    case AC_START:
			// frequency must be > 0
			if (value->rValue <= 0.0) {
#ifdef	TRACE_
			printf("\n[%s : line %d] Set_Parameter(): AC_START <= 0 is invalid.", 
					__FILE__, __LINE__);
#endif
				ACstartFreq = 1.0;
				return  (-1);
			}

			ACstartFreq = value->rValue;
#ifdef	TRACE_
			printf("\n[%s : line %d] Set_Parameter(): ACstartFreq = %g", 
					__FILE__, __LINE__, ACstartFreq);
#endif
			break;

	    case AC_STOP:
		if (value->rValue <= 0.0) {
#ifdef	TRACE_
			printf("\n[%s : line %d] Set_Parameter(): AC_STOP <= 0 is invalid.", 
					__FILE__, __LINE__);
#endif
			ACstartFreq = 1.0;
			return  (-1);
		}

	        ACstopFreq = value->rValue;
#ifdef	TRACE_
			printf("\n[%s : line %d] Set_Parameter(): ACstopFreq = %g", 
					__FILE__, __LINE__, ACstopFreq);
#endif
	        break;

	    default:
#ifdef	TRACE_
			printf("\n[%s : line %d] Set_Parameter(): Parameter not defined.", 
					__FILE__, __LINE__);
#endif
		return  (-1);
    }
    return  (OK);
}

/**
   PrintFreqResponse()
   Prints freq response to the console and an output file simultaneously.
*/
void
Analysis_AC::PrintFreqResponse(char *tofilename)
{
	cout << LINE_INFO << endl;

	if (tofilename) {
		ListFile.open(tofilename);
	}
	
	cout << "\n  Frequency response:" << endl;
	
	cout << setw(4) << "#" << setw(14) << "Hz" 
		<< setw(30) << "Complex Value" << endl;
//		<< setw(20) << "Mag" << setw(30) << "Complex Value" << endl;

	if (tofilename) {
		ListFile << setw(4) << "#" << setw(14) << "Hz" 
			<< setw(20) << "Mag" << setw(30) << "Complex Value" << endl;
	}

#if 00
	// Reference code !!

	cout << "\nFrequency response:" << endl;
	cout << setw(10) << "Hz" << setw(20) << "Mag" 
		 << setw(20) << "Phase"
		 << setw(20) << "db20"
		 << setw(30) << "Complex Value" << endl;*/

	dataFile << setw(10) << "Hz" << setw(20) << "Mag" 
		 << setw(30) << "Complex Value" << endl;*/
	
	//******************************************
	// magn(), arg(), db20() etc. are defined in "complex.h".
	//******************************************
	cout << setw(20) << setprecision(4) 
		 << setiosflags(ios::fixed) << magn(H)
		 << setw(20) << setprecision(4) 
		 << setiosflags(ios::fixed) << arg(H)
		 << setw(20) << setprecision(4) 
		 << setiosflags(ios::fixed) << db20(H);
	dataFile << setw(20) << setprecision(4) 
		 << setiosflags(ios::fixed) << magn(H)
		 << setw(20) << setprecision(4) 
		 << setiosflags(ios::fixed) << arg(H)
		 << setw(20) << setprecision(4) 
		 << setiosflags(ios::fixed) << db20(H) << endl;

	cout << setw(20) << setprecision(6) << H.re <<" +j" << H.im << endl;
	dataFile << setw(20) << setprecision(6) << H.re <<" +j" << H.im << endl;
#endif

	for (int i = 0; i < NumFreqPoints; i++)
	{
		cout << setiosflags(ios::scientific) 
			 << setw(4) << setprecision(0) << i+1
			 << setw(14) << setprecision(4) << FreqVec[i]
//		 	 << setw(20) << setprecision(4) << MagVec[i]
		 	 << setw(20) << setprecision(6) << RespVec[i] << endl;
		
		if (tofilename) {
			ListFile << setiosflags(ios::scientific) 
			 << setw(4) << setprecision(0) << i+1
			 << setw(14) << setprecision(4) << FreqVec[i]
		 	 << setw(20) << setprecision(4) << MagVec[i]
		 	 << setw(20) << setprecision(6) << RespVec[i] << endl;
		 }
	}
	
	if (tofilename) {
		cout << "\nFrequency response has been written to the file: \"" 
			 << tofilename << "\"" << endl;
		ListFile.close();
	}

} // PrintFreqResponse()


void Analysis_AC::PrintCoefficients()
{
	ExpGPDDnode *pExpRoot = pGPDD->ExpGPDDroot;

	vector<sCoeffNode*> 	*psCoeff = pExpRoot->Vertex_IN->sCoeffList;
	vector<int>         		*psPower = pExpRoot->Vertex_IN->sPowerList;

	cout << LINE_INFO << "\n\tPrint the Numerator coefficients: " << endl;
	for (int i = 0; i < psCoeff->size(); i ++)
	{
		cout << "    Coeff. of s^(" << (*psPower)[i] << ") = " <<  (*psCoeff)[i]->GetValue() << endl; 
	}

	psCoeff = pExpRoot->Vertex_OUT->sCoeffList;
	psPower = pExpRoot->Vertex_OUT->sPowerList;

	cout << LINE_INFO << "\n\tPrint the Denominator coefficients: " << endl;
	for (int i = 0; i < psCoeff->size(); i ++)
	{
		cout << "    Coeff. of s^(" << (*psPower)[i] << ") = " <<  (*psCoeff)[i]->GetValue() << endl; 
	}
}

/**
   PrintFreqResponse()
   Prints freq response to the console and an output file simultaneously.
*/
void
Analysis_AC::PrintFreqResponse2(complex_t *respvec)
{
	cout << LINE_INFO << endl;

	if (!respvec) {
		return;
	}
	
	cout << "\n  Frequency response:" << endl;
	
	cout << setw(4) << "#" << setw(14) << "Hz" 
		<< setw(30) << "Complex Value" << endl;
//		<< setw(20) << "Mag" << setw(30) << "Complex Value" << endl;

	/*if (tofilename) {
		ListFile << setw(4) << "#" << setw(14) << "Hz" 
			<< setw(20) << "Mag" << setw(30) << "Complex Value" << endl;
	}*/

	for (int i = 0; i < NumFreqPoints; i++)
	{
		cout << setiosflags(ios::scientific) 
			 << setw(4) << setprecision(0) << i+1
			 << setw(14) << setprecision(4) << FreqVec[i]
//		 	 << setw(20) << setprecision(4) << MagVec[i]
		 	 << setw(20) << setprecision(6) << respvec[i] << endl;
		
		/*if (tofilename) {
			ListFile << setiosflags(ios::scientific) 
			 << setw(4) << setprecision(0) << i+1
			 << setw(14) << setprecision(4) << FreqVec[i]
		 	 << setw(20) << setprecision(4) << MagVec[i]
		 	 << setw(20) << setprecision(6) << RespVec[i] << endl;
		 }*/
	}
	
	/*if (tofilename) {
		cout << "\nFrequency response has been written to the file: \"" 
			 << tofilename << "\"" << endl;
		ListFile.close();
	}*/

} // PrintFreqResponse()



/**
	Print the sensitivity response to the console and to an output file if given.
*/
void
Analysis_AC::PrintSensResponse(char *tofilename)
{
	cout << LINE_INFO << endl;

	if (tofilename) {
		ListFile.open(tofilename);
	}
	
	cout << "\n  Sensitivity response:" << endl;
	
	cout << setw(4) << "#" << setw(14) << "Hz" 
		<< setw(30) << "Complex sensitivity" << endl;

	if (tofilename) {
		ListFile << setw(4) << "#" << setw(14) << "Hz" 
			<< setw(30) << "Complex sensitivity" << endl;
	}

	for (int i = 0; i < NumFreqPoints; i++)
	{
		cout << setiosflags(ios::scientific) 
			 << setw(4) << setprecision(0) << i+1
			 << setw(14) << setprecision(4) << FreqVec[i]
//		 	 << setw(20) << setprecision(4) << MagVec[i]
		 	 << setw(20) << setprecision(6) << SensVec[i] 
		 	 << endl;
		
		if (tofilename) {
			ListFile << setiosflags(ios::scientific) 
			 << setw(4) << setprecision(0) << i+1
			 << setw(14) << setprecision(4) << FreqVec[i]
//		 	 << setw(20) << setprecision(4) << MagVec[i]
		 	 << setw(20) << setprecision(6) << SensVec[i] 
		 	 << endl;
		 }
	}
	
	if (tofilename) {
		cout << "\nSensitivity response has been written to the file: \"" 
			 << tofilename << "\"" << endl;
		ListFile.close();
	}

} // PrintFreqResponse()

