/**
* \mainpage
* \section Usage
* @code
* PowFit2D [-v] [-p maxdegree] inputfilename > outputfile
* @endcode
* The @c inputfilename is a list of function values on
* a finite discrete set of 2D Cartesian coordinates as described in PowFit2D::readf().
* The standard output contains the fitting coefficients as described in
* PowFit2D::display(). If one wants to delete the comment lines that are dispersed
* in the standard output, one can use a pipe like
* @code
* PowFit2D [-v] [-p maxdegree] inputfilename | fgrep -v '#' | awk '{print $1 $2 $3}' > outputfile
* @endcode
* One can extract the original input and the fit by using the -v option
* and looking at the lines
* between FITTED START and FITTED END:
* @code
* PowFit2D -v inputfilename | sed -n '/FITTED START/,/FITTED END/p' | sed '/FITTED/d' | tr "#" " " > outputfile
* @endcode
* and this output is immediately useful for gnuplot(1)'s splot command.
* @code
* gnuplot
* gnuplot> splot "outputfile" title "original", "outputfile" using 1:2:4 title "fitted"
* gnuplot> quit
* @endcode
*
* \subsection Command line options
* See PowFit2D::main() .
*
* \section Description
* See PowFit2D::fit() for the fitting formul.
* 
* \section Compilation
* A prerequisite is that the LAPACK library of <a href="http://www.netlib.org/lapack">netlib.org</a>
* has been installed.
*  @code
*    g++ -o PowFit2D -I<dir_of_clapack.h> PowFit2D.cpp -L<dir_of_liblapack> -llapack
*  @endcode
*  The makefile entry is
*  @code
*    LIBFLAGS=-L<dir_of_liblapack>
*    INCFLAGS=-I<dir_of_clapack.h>
*    PowFit2D: PowFit2D.cpp
*         $(CXX) -o $@ $(INCFLAGS) $< $(LIBFLAGS) -llapack
*  @endcode
* IF the <a href="http://www.gnu.org/software/gsl/">GSL</a> library is available, this ought be indicated by
* setting preprocessor variable HAVE_GSL, then the LAPACK library is not needed
*  @code
*    g++ -o PowFit2D -I<dir_of_gsl> -DHAVE_GSL PowFit2D.cpp -L... -lgslcblas -lgsl
*  @endcode
*
* @author Richard J. Mathar
* @todo compute error bars on the fitting coefficients
* @since 2007-07-12
* <a href="http://www.strw.leidenuniv.nl/~mathar">Richard J. Mathar homepage</a>.
* @see
* \cite{GuoOptik117}
*/
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <unistd.h>
#include <string.h>
#include <regex.h>
#include <time.h>

#ifdef HAVE_GSL
#include <gsl/gsl_sf.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

#else /* HAVE_GSL */

extern "C" {
/** The f2c header file f2c.h.
* A prerequisite to include clapack.h below.
* This is obtained from <a href="http://www.netlib.org/f2c">www.netlib.org/f2c</a>.
*/
#include <f2c.h>
}

extern "C" {
/** The clapack header file \c clapack.h.
* This is obtained from <a href="http://www.netlib.org/clapack">www.netlib.org/clapack</a>.
*/
#include <clapack.h>
}

#endif /* HAVE_GSL */

using namespace std ;

/** The class contains a list of 2D points \f$(x_k,y_k)\f$, \f$k=1,..,N\f$, and scalar
* functions values \f$f_k\f$ of these, and performs a least squares fit on a polynomial
* basis in the sense of \f$f_k\approx \sum_{i,j} \alpha_{i,j} x_k^i y_k^j\f$.
* @since 2007-07-12
*/
class PowFit2D {
public:

	/** Ctor.
	* This initializes the PowFit2D::x, PowFit2D::y and PowFit2D::f vectors with
	* the triples of floating point numbers that are to be fitted.
	* @param[in] fname the file name of the input file with the tabulated (x,y,f) triples.
	* @since 2007-07-12
	*/
	PowFit2D(char *fname, const bool verbos=false) : verb(verbos)
	{
		readf(fname,verbos) ;
	}

	/** Do the fitting.
	* Minimize
	* \f$ \sum_k [f(x_k,y_k)-f_k]^2\rightarrow \min \f$.
	* with the ansatz
	* \f$ f(x,y)= \sum_{p=0}^n \sum_{i=0}^p \alpha_{i,p-i} x^i y^{p-i} \f$.
	* which leads to the linear system of equations for the unknown \f$ \alpha_{i,p-i}\f$
	* \f$ \sum_{p,i}\alpha_{i,p-i} \sum_k x_k^iy_k^{p-i} x_k^l y_k^m = \sum_k f_k x_k^l y_k^m \f$.
	* @param[in] maxP the maximum power in the fitting procedure.
	*  The default of -1 means that it is chosen automatically to be the maximum
	*  allowed by the number of 2D points that have been read in.
	* @since 2007-07-12
	*/
	void fit(int maxP = -1)
	{
		/* the degrees of freedom as determined by the number of points.
		*/
		const int deg = x.size() ;

		if( deg <=0)
		{
			cerr << "Number of points " << deg << " insufficient\n"; 
			exit(EXIT_FAILURE) ;
		}

		/* determine maxP automatically. deg=0 does not allow fitting.
		* deg=1,2 means maxP=0, deg=3..5 means maxP=1 etc.
		*/
		if ( maxP < 0 )
		{
			maxP = 0 ;
			while ( triang(maxP+2) <= deg)
				maxP++ ;
		}
		else
		{
			while ( triang(maxP+1) > deg)
				maxP-- ;
		}

		maxHybrP = maxP ;

		/* the number of fitting degrees of freedom
		* We keep this an overdetermined system, n<=deg. The case of n=deg
		* means the fit is actually an interpolation (ie, exact).
		*/
		int n=triang(maxP+1) ;

		if (verb)
		{
			cout << "# maximum total power " << maxP << " for " << deg << " input points; " 
				<< n << " fitting degrees\n" ;
		}

#ifdef HAVE_GSL
		/* the n by n matrix of the linear system of equations */
		gsl_matrix *A=gsl_matrix_alloc(n,n) ;
		gsl_matrix_set_all(A,0.) ;

		/* the RHS of the linear system of equations */
		gsl_vector *rhs=gsl_vector_alloc(n) ;
		gsl_vector_set_all(rhs,0.) ;

		/* fill the matrix and the right hand side
		*/
		for(int row=0; row<n ;row++)
		{
			int rowi, rowj ;
			strIndxInv(row,rowi,rowj) ;
			for(int col=0; col< n ;col++)
			{
				int coli, colj ;
				strIndxInv(col,coli,colj) ;
				/* debugging
				* cout << "# mapped " << row << " " << rowi << " " << rowj << endl ;
				*/
				for(int k=0 ; k < deg ;k++)
					* gsl_matrix_ptr(A,row,col) += pow(x[k],rowi+coli)*pow(y[k],rowj+colj) ;
			}
			for(int k=0 ; k < deg ;k++)
				*gsl_vector_ptr(rhs,row) += f[k]*pow(x[k],rowi)*pow(y[k],rowj) ;
		}

		/* solve the linear system of equations
		*/
		gsl_linalg_HH_svx(A,rhs) ;

		/* pack the expansion coefficients into the general space 
		*/
		alpha.clear() ;
		for(int row=0; row < n ; row++)
			alpha.push_back(gsl_vector_get(rhs,row)) ;

		gsl_vector_free(rhs) ;
		gsl_matrix_free(A) ;
#else /* HAVE_GSL */

/* the macro maps a pair of Fortran77 0-based coordinates to a 
* 0-based 1D C index.
*/
#define POWFIT2D_F77_TO_C(row,col,N) ((col)*(N)+(row))
		/* the n by n matrix of the linear system of equations */
		double *A=new double[n*n] ;

		/* the RHS of the linear system of equations */
		double *rhs=new double[n] ;

		/* fill the matrix and the right hand side
		*/
		for(int row=0; row<n ;row++)
		{
			int rowi, rowj ;
			strIndxInv(row,rowi,rowj) ;
			for(int col=0; col< n ;col++)
			{
				int coli, colj ;
				strIndxInv(col,coli,colj) ;
				/* debugging
				* cout << "# mapped " << row << " " << rowi << " " << rowj << endl ;
				*/
				A[POWFIT2D_F77_TO_C(row,col,n)] =0. ;
				for(int k=0 ; k < deg ;k++)
					A[POWFIT2D_F77_TO_C(row,col,n)] += pow(x[k],rowi+coli)*pow(y[k],rowj+colj) ;
			}
			rhs[row]=0. ;
			for(int k=0 ; k < deg ;k++)
				rhs[row] += f[k]*pow(x[k],rowi)*pow(y[k],rowj) ;
		}

		/* solve the linear system of equations
		*/

		integer nrhs=1 ;
		integer lda= n ;
		integer *ipiv = new integer[n] ;
		integer info = 0 ;
		dgesv_(&lda,&nrhs,A,&lda,ipiv,rhs,&lda,&info) ;
		if ( info )
			cerr << __FILE__ << " " << __LINE__ << " dgesv(1) info= " << info << endl ;
		delete [] ipiv ;

		/* pack the expansion coefficients into the general space 
		*/
		alpha.clear() ;
		for(int row=0; row < n ; row++)
			alpha.push_back(rhs[row]) ;

		delete [] rhs ;
		delete [] A ;
#undef POWFIT2D_F77_TO_C
#endif /* HAVE_GSL */
	}

	/** Show the fitted function on stdout.
	* The output format is:
	* - lines that start with hashes (#) are comment lines. Their number
	*   depends on whether the -v command line option was used or  not.
	*   - If the -v option was used, they contain a list of the original (x,y,f)
	*     triples, each time followed by the value of the fitted bivariate and
	*     the deviation between the fit and the original f. This optional section
	*     is started with a line POWFIT2D FITTED START and continues up to a line
	*     that contains POWFIT2D FITTED END.
	*   - The comment lines contain a section started with POWFIT2D SPHERICAL START
	*     and concluded with POWFIT2D SPHERICAL END which show the fitting function
	*     transformed to spherical coordinates, \f$x=r\cos(\varphi)\f$,
	*     \f$y=r\sin(\varphi)\f$. Each of these lines is a term in the expansion
	*     \f$f=\sum \beta_k r^p (\sin,\cos)(m\varphi)\f$ encoded as follows: the
	*     first number is the integer \f$p\f$. The second number is the value of \f$m\f$.
	*     If this value is positive or zero, the term is \f$\propto \cos(m\varphi)\f$,
	*     if this value is negative, the term is \f$\propto\sin(|m|\varphi)\f$. The
	*     third number is the value of \f$\beta\f$ for this term, which may be zero.
	*     This functionality is equivalent to a lookup via the table in
	*     <a href="http://www.strw.leidenuniv.nl/~mathar/progs/Zernike.txt">Zernike.txt</a>.
	* - For lines that are not comment lines, the first non-negative integer is the power \f$i\f$
	*   in the x coordinate, the next
	*   non-negative integer is the power \f$j\f$ in the y coordinate, and the following number
	*   is the expansion coefficient \f$\alpha_{i,j}\f$ in \f$f\approx \sum_{i,j}\alpha_{i,j}x^iy^j\f$.
	*   The rest of the line is an expression summarizing this in
	*   more human-readable fashion.
	* @warning this only makes sense if PowFit2D.fit() had been called before.
	* @parameter[in] verb verbosity level
	* @since 2007-07-12
	*/
	void display(const bool verb=false)
	{
		/* Display the main result, which is the expansion coefficients
		* and their associated powers in the two cartesian coordinates.
		*/
		for(int i=0; i < alpha.size() ; i++)
		{
			int powx, powy ;
			strIndxInv(i,powx,powy) ;
			cout << powx << " " << powy << " " << alpha[i] << " " 
				<< alpha[i] << "*x^" << powx << "*y^" << powy << endl ;
		}

		/* If verbose output is demanded, also show the fit of the quality
		* by providing the original triples together with their fitted
		* approximation by the bi-variate polynomial and the error between
		* the fit and the input value.
		*/
		if ( verb)
		{
			cout << "# POWFIT2D FITTED START\n" ;
			for(int i=0; i < x.size() ; i++)
			{
				const double fi = fitted(x[i],y[i]) ;
				cout << "# " << x[i] << " " << y[i] << " " << f[i] << " " << fi << " " << fi-f[i] << endl ;
			}
			cout << "# POWFIT2D FITTED END\n" ;
		}

		/* Finally convert the expansion coefficients in the Cartesian
		* coordinates to expansion coefficients in the spherical system.
		*/
		cout << "# POWFIT2D SPHERICAL START x=r*cos phi, y=r*sin phi\n" ;

		/* loop over the powers r^p
		*/
		for(int p=0; p <= maxHybrP ; p++)
		{
			/* loop over the sin(m*phi) with m<0 and the cos(m*phi) with m>=0.
			*/
			for(int m= p; m >= -p ; m--)
			{
				double co = 0. ;
				/* gather the contribution of all mixed powers individually.
				* The powx+powy=p that contribute with x^powx*y^powy are all in a consecutiave
				* range of straightened indices.
				*/
				for(int s= triang(p) ; s < triang(p+1) ; s++)
				{
					int powx, powy ;
					strIndxInv(s,powx,powy) ;
					const int degfact = ( m==0) ? 1 : 2 ;
					co += degfact*alpha[s]*cosFlatn(powx,powy,m) ;
				}
				cout << "# " << p << " " << m << " " << co << " " << co << "*r^" << p ;
				if ( m > 0)
					cout << "*cos(" << m << "*phi)" << endl ;
				else if ( m < 0 )
					cout << "*sin(" << -m << "*phi)" << endl ;
				else
					cout << endl ;
			}
		}
		cout << "# POWFIT2D SPHERICAL END\n" ;
	}

protected:
private:

	/** list of x coordinates
	*/
	vector<double> x ;

	/** list of y coordinates, one per x coordinate
	*/
	vector<double> y ;

	/** list of function values, one per (x,y) pair.
	*/
	vector<double> f ;

	/** verbosity rised if true.
	*/
	bool verb ;

	/** list of expansion coefficients, 1D version (straightened)
	*/
	vector<double> alpha ;

	/** the maximum combined power of the x and y coordinates for the fit
	*/
	int maxHybrP ;

	/** classify the input line as comment or non-comment.
	*/
	static bool isComment(const string inpli, const regex_t & preg)
	{
		/* debugging
		* cout << __FILE__ << " " << __LINE__ << " " << inpli << endl ;
		*/
		if ( regexec(&preg,inpli.c_str(),0,0,0) == 0)
			return true ;
		else
			return false ;
	}

	/** The binomial C(a,b).
	* \latexonly
	* \cite[3.1.2]{AS}
	* \endlatexonly
	* @param[in] a the non-negative numerator of the binomial
	* @param[in] b the non-negative denomiator, less or equal to a
	* @return  the binomial coefficient \f$C(a,b)=a!/b!/(a-b)!\f$
	* @since 2007-07-12
	* @see <a href="http://research.att.com/~njas/sequences/A007318">OEIS A007318</a>
	*/
	static double binomial(int a, int b)
	{
#ifdef HAVE_GSL
		return gsl_sf_choose(unsigned(a),unsigned(b)) ;
#else

		if ( a<0 || b<0 || b>a)
			return 0. ;
		if (a-b < b)
			b = a-b ;
		if ( b == 0 )
			return 1. ;
		double resul = a ;
		while ( b-- > 1)
			resul *= double(a-b+1)/double(b) ;
		return a ;
#endif /* HAVE_GSL */
	}

#if 0
	/** sign function.
	* @param[in] x the value of the argument
	* @return 1 if x>=0, -1 if x<0.
	* @since 2007-07-12
	*/
	inline static int sign(const double x)
	{
		return (x> 0.) ? 1 (x==0. ? 0 : -1) ;
	}
#endif

	/** The contribution of \f$ \cos^l(\varphi) \sin^k(\varphi)\f$ to \f$\cos(m \varphi)\f$
	*  or \f$\sin(m\varphi)\f$. With the Euler formula
	* \latexonly
	* \cite[4.3]{AS}
	* \endlatexonly
	* and the binomial expansion
	* \latexonly
	* \cite[3.1.3]{AS}
	* \endlatexonly
	* we have
	* \f$\cos^l\varphi \sin^k\varphi=(e^{i\varphi}+e^{-i\varphi})^l/(2^l) (e^{i\varphi}-e^{-i\varphi})^k/(2i)^k\f$
	* \f$\sim\frac{1}{2^{k+l}}(-)^{\langle k/2\rangle}\sum_{s=0}^l{l \choose s}\sum_{t=0}^k{k \choose t}(-)^t e^{i(2s-l+k-2t)\varphi}\f$
	* to be divided by \f$i\f$ if \f$k\f$ is odd. The algorithm is simply to match
	* the integer in the exponential with \f$m\f$.
	* @param[in] l the power of the cosine
	* @param[in] k the power of the sine
	* @param[in] m the cycle number of the cosine (m>=0) or the sine (m<0)
	* for which the coefficient is created.
	* @return the parameter beta in the expansion cos^l phi sin^k phi = .... beta*cos(m*phi)
	*   or = .. beta*sin(|m|*phi)
	* @since 2007-07-12
	*/
	static double cosFlatn(const int l, const int k, int m)
	{
		double resul=0. ;
		/* for nonzero result, an even k must meet a cosine expansion coeffieicnt
		* and an odd k must meet a sine expansion coefficient.
		*/
		if ( m >=0  && k % 2 == 0  || m < 0 && k % 2)
		{
			m = abs(m) ;

			/* k/2 if k is even , floor(k/2) if k is odd
			*/
			const int khalf = k /2 ;
			if ( (m+l+k) % 2 ==0)
			{
				for(int s=0 ; s <= l ; s++)
				{
					const int t= s-(m+l-k)/2 ;
					if ( t>=0 && t <= k)
						if ( t % 2 )
							resul -= binomial(l,s)*binomial(k,t) ;
						else
							resul += binomial(l,s)*binomial(k,t) ;
				}
				if ( khalf % 2)
					resul *= -1. ;
				resul /= pow(2.,l+k) ;
			}
		}
		return resul ;
	}

	/** Calculate a value of the fitting function.
	* @param[in] x Cartesian x coordinate of the 2D point
	* @param[in] y Cartesian y coordinate of the 2D point
	* @return the mixed polynomial (the fit) at the 2D point.
	* @since 2007-07-12
	*/
	double fitted(const double x, const double y)
	{
		double f=0. ;
		for(int i=0; i < alpha.size() ; i++)
		{
			int powx, powy ;
			strIndxInv(i,powx,powy) ;
			f += alpha[i]*pow(x,powx)*pow(y,powy) ;
		}
		return f ;
	}

	/** Read the (x,y,f) triples from an ASCII file.
	* The syntax of the lines in the ASCII file is as follows:
	* - Lines that start with optional white space followed by a hash (#) or
	*   containing only white space are ignored (a.k.a. comment lines)
	* - All other lines start with three floating point numbers separated
	*   by white space. Letters and numbers in the rest of these lines are
	*   ignored. The three floating point numbers are the x Cartesian coordinate,
	*   then the y Cartesian coordinate, then the function value f.
	*   The count of lines of this form implies the number of input triples.
	* @param[in] fname the UNIX file name
	* @param[in] verb if true turns on verbose commenting on the inputs
	* @since 2007-07-12
	*/
	void readf(char *fname, bool verb=false)
	{
		ifstream inf(fname) ;
		/* complain if the file is not readable.
		*/
		if ( !inf)
		{
			cerr << "Cannot open " << fname << endl ;
			exit(EXIT_FAILURE) ;
		}
		
	

		/* a regular expression characterizing comment lines
		*/
		regex_t preg;
		if ( int rege = regcomp(&preg,"(^[[:space:]]*#)|(^[[:space:]]*$)",REG_NOSUB|REG_EXTENDED) )
		{
			char errbuf[257] ;
			regerror(rege,&preg,errbuf,257) ;
			cerr << __FILE__ << " " << __LINE__ << " " << errbuf << endl ;
		}

		/* read the file line by line until EOF.
		*/
		while( !inf.eof() )
		{
			/* one line of the input file
			*/
			string inpli ;
			getline(inf,inpli) ;
			if ( !isComment(inpli,preg) )
			{
				double vals[3] ;
				istringstream s(inpli) ;
				s >> vals[0] >> vals[1] >> vals[2] ;
				if ( verb)
					cout << "# " << x.size() << " " << vals[0] << " " << vals[1] << " " << vals[2] << endl ;
				x.push_back(vals[0]) ;
				y.push_back(vals[1]) ;
				f.push_back(vals[2]) ;
			}
		}
		regfree(&preg) ;
	}

	/** Compute  triangular number.
	* <a href="http://research.att.com/~njas/sequences/A000217">OEIS A000217</a>
	* \latexonly
	* \cite{EIS}
	* \endlatexonly
	* @param[in] i
	* @return the i'th triangular number, i(i+1)/2.
	* @since 2007-07-12
	*/
	int triang( const int i)
	{
		return i*(i+1)/2 ;
	}

	/** Derive a straightened single index from a pair index. The result is
	* constructed by reading the infinite 2x2 array of integer indices
	* in the first quadrant along anti-diagonals.
	* @param[in] i first index of the pair, >=0
	* @param[in] j second index of the pair, >=0
	* @return a unique number larger or equal to zero. (0,0) is mapped on 0,
	*  (0,1) on 1, (1,0) on 2, (0,2) on 3, (1,1) on 4, (2,0) on 5 etc. So a group
	*  of common i+j produces a sequential list of integers, and the smallest i
	*  are maped on the lower single indices.
	* @since 2007-07-12
	*/
	int strIndx( const int i, const int j)
	{
		return triang(i+j)+i ;
	}

	/** Decompose a straightened single index into a pair of indices.
	* This is the inverse function of strIndx().
	* @param[in] n the straightened single index.
	* @param[out] i first index of the pair, >=0
	* @param[out] j second index of the pair, >=0
	* @since 2007-07-12
	*/
	void strIndxInv( const int n, int &i, int & j)
	{
		/* the triangular number that starts a sequence of  common sum i+j
		*/
		int p=0 ;
		while ( triang(p) <= n)
			p++ ;
		p-- ;
		i = n-triang(p) ;
		j=p-i ;
	}
} ;

/** Print a usage information.
* @param[in] argv0 the command name
* @since 2007-07-12
*/
void usage(char *argv0)
{
        cout << "usage : " << argv0 << " [-v] [-p maxdegree] inptfile\n" ;
}

/**  The main executable.
* @param[in] -v the option \c -v can be used to produce more
*   verbose output
* @param[in] -p the option \c -p followed by a non-negative integer indicates
*   that the maximum (hybrid) degree of the powers in the fit should be reduced
*   to this integer. If not used, the program will calculate a default, which
*   is to take the maximum fitting degree which still matches a least-squares
*   algorithm. The default means that the degree is 0 for 1 Cartesian input point,
*   the degree is 1 for 2 or 3 input points, the degree is 2 for 4 to 6 input
*   points etc. This can be seen as always including a full set of azimuthal
*   parameters \f$m\f$ for each \f$\sim r^p\cos(m\varphi)\f$ and \f$\sim r^p\sin(m\varphi)\f$,
*   \f$0\le m\le p\f$, such that no directional preference is found in the
*   algorithm.
* @param[in] fname the single command line argument is the ASCII file name
*   with input lines as described in PowFit2D::readf().
* @return 0 if the fitting was apparently successful, -1 if wrong options
*   were used or the input file couldn't be opened.
* @since 2007-07-12
*/
int main(int argc, char *argv[])
{
	/* command line option.
	*/
	int c ;

	/* Maximum radial power in the fit. Initialized with negative value
	* means we automate the calculation.
	*/
	int maxP = -1 ;

	/* verbosity
	*/
	bool verb=false ;

	/* read command line options
	*/
	while( (c=getopt(argc,argv,"vp:")) != -1 )
        {
                switch(c)
                {
                case 'v':
			verb=true ;
			break ;
                case 'p':
			maxP = atoi(optarg) ;
			break ;
                case '?':
                        usage(argv[0]) ;
                        return -1 ;
                        break ;
                }
        }

	if ( optind >= argc)
	{
		usage(argv[0]) ;
		cerr << "Missing file name argument\n" ;
		return -1 ;
	}

	if (verb)
	{
		time_t now=time(NULL) ;
		cout << "# " << getenv("USER") << " " << ctime(&now) << "# " ;
		for(c=0 ; c < argc; c++)
			cout << argv[c] << " " ;
		cout << endl ;
	}

	/* read the 2D points to be fitted from the ASCII File
	*/
	PowFit2D fitProbl(argv[optind],verb) ;

	/* do the fitting.
	*/
	fitProbl.fit(maxP) ;

	/* show the results
	*/
	fitProbl.display(verb) ;

	return 0 ;
}