#include "CSCE_UA.H"

#include <math.h>
#include <time.h>
#include <limits.h>
#include <algorithm>

#define MAX_NUM 255
#define TWENTY 20			// なぜ20なのか不明

using namespace std;

CSCE_UA::CSCE_UA(double dErrorValue)
{
	m_dErrorValue = dErrorValue;
	m_iOptimizedPrm = 0;
}

bool CSCE_UA::sceua(model_input& cModel, const char strOutFile[], int iCurrentNumber,
		void (*model_func) (std::vector<double>& , const std::vector< std::vector<double> >& ,  std::vector< std::vector<double> >&) )
{
	// SHUFFLED COMPLEX EVOLUTION METHOD FOR GLOBAL OPTIMIZATION
	// -- Version 2.1  by QINGYUN DUAN

	// STATEMENT BY AUTHOR:
	// --------------------
	//
	// This general purpose global optimization program is developed at the Department of Hydrology & Water Resources of the University of Arizona.  
	// urther information regarding the SCE-UA method can be obtained from Dr. Q. Duan, Dr. S. Sorooshian or Dr. V.K. Gupta at the address and phone number listed above.  
	// We request all users of this program make proper reference to the paper entitled 
	// 'Effective and Efficient Global Optimization for Conceptual Rainfall-runoff Models' by Duan, Q., S. Sorooshian, and V.K. Gupta, Water Resources Research, Vol 28(4), pp.1015-1031, 1992.
	//
	//
	// LIST OF SCE ALGORITHMIC CONTROL PARAMETERS:
	//
	// ngs = number of complexes in the initial population
	// npg = number of points in each complex
	// npt = total number of points in initial population (npt=ngs*npg)
	// nps = number of points in a sub-complex
	// nspl = number of evolution steps allowed for each complex before complex shuffling
	// mings = minimum number of complexes required, if the number of complexes is allowed to reduce as the optimization proceeds
	// iseed = initial random seed
	// iniflg = flag on whether to include the initial point in population
	//    = 0, not included
	//    = 1, included
	// iprint = flag for controlling print-out after each shuffling loop
	//    = 0, print information on the best point of the population
	//    = 1, print information on every point of the population
	//
	// CONVERGENCE CHECK PARAMETERS
	//
	// maxn = max no. of trials allowed before optimization is terminated
	// kstop = number of shuffling loops in which the criterion value must chang by the given percentage before optimization is terminated
	// pcento = percentage by which the criterion value must change in given number of shuffling loops
	// ipcnvg = flag indicating whether parameter convergence is reached  (i.e., check if gnrng is less than 0.001)
	//    = 0, parameter convergence not satisfied
	//    = 1, parameter convergence satisfied
	//
	// LIST OF LOCAL VARIABLES
	// x(.,.) = coordinates of points in the population
	// xf(.) = function values of x(.,.)
	// xx(.) = coordinates of a single point in x
	// cx(.,.) = coordinates of points in a complex
	// cf(.) = function values of cx(.,.)
	// s(.,.) = coordinates of points in the current simplex
	// sf(.) = function values of s(.,.)
	// bestx(.) = best point at current shuffling loop
	// dBestFunc = function value of bestx(.)
	// worstx(.) = worst point at current shuffling loop
	// dWorstFunc = function value of worstx(.)
	// xnstd(.) = standard deviation of parameters in the population
	// gnrng = normalized geometric mean of parameter ranges
	// lcs(.) = indices locating position of s(.,.) in x(.,.)
	// bound(.) = bound on ith variable being optimized
	// ngs1 = number of complexes in current population
	// ngs2 = number of complexes in last population
	// iseed1 = current random seed
	// criter(.) = vector containing the best criterion values of the last 10 shuffling loops

	//  LOCAL ARRAYS
	//      dimension x(2000,16),xx(16),bestx(16),worstx(16),xf(2000)
	//      dimension s(50,16),sf(50),lcs(50),cx(2000,16),cf(2000)
	//      dimension xnstd(16),bound(16),criter(20),unit(16)

	vector< vector<double> > x;		// coordinates of points in the population;
	vector<double> xf;				// function values of x(.,.)
	vector<double> xx;				// coordinates of a single point in x
	vector<double> bound;			// bound on ith variable being optimized
	vector<double>  unit;			// 
	vector<double>  bestx;			// best point at current shuffling loop
	vector<double>  worstx;			// worst point at current shuffling loop
	double dBestFunc;				// function value of bestx(.)
	double dWorstFunc;				// function value of worstx(.)
	vector<double>  xnstd;			// standard deviation of parameters in the population
	double  gnrng;					// normalized geometric mean of parameter ranges
	int ipcnvg;						// flag indicating whether parameter convergence is reached
	vector< vector<double> > cx;	// coordinates of points in a complex
	vector<double> cf;				// function values of cx(.,.)
	vector<int> lcs;				// indices locating position of s(.,.) in x(.,.)
	vector< vector<double> > s;		// coordinates of points in the current simplex
	vector<double> sf;				// function values of s(.,.)
	vector<double> criter;			// vector containing the best criterion values of the last 10 shuffling loops
	int ngs2;						// number of complexes in last population

//	char strLine[MAX_NUM];	
	int i, j, k, k1, k2, l, lpos;
	double rand;
	int iTerminatedError = 0;
	bool bIsOK = true;

	printf(" ENTER THE SCEUA SUBROUTINE --- \n"); 

	x  = vector< vector<double> >(m_iRecordNum, vector<double>(m_nopt, 0) );
	cx = vector< vector<double> >(m_iRecordNum, vector<double>(m_nopt, 0) );
	xf.assign(m_iRecordNum);
	cf.assign(m_iRecordNum);
	xx.assign(m_nopt);
	unit.assign(m_nopt);
	bestx.assign(m_nopt);
	worstx.assign(m_nopt);
	xnstd.assign(m_nopt);
	lcs.assign(m_nps);
	s.assign(m_nopt);
	s  = vector< vector<double> >(m_nps, vector<double>(m_nopt, 0) );
	sf.assign(m_nps);
	criter.assign(TWENTY);
	
	// *******************************************************
	// INITIALIZE VARIABLES
	// *******************************************************
	int nloop = 0;
	int loop = 0;
	int igs = 0;
	int nopt1 = 8;
	if (m_nopt < 8) nopt1 = m_nopt;
	int nopt2 = 12;
	if (m_nopt < 12) nopt2 = m_nopt;
	
	// *******************************************************
	// INITIALIZE RANDOM SEED TO A NEGATIVE INTEGER
	// *******************************************************
	srand(m_jseed[iCurrentNumber]);
	int iseed1 = m_jseed[iCurrentNumber];

	// *******************************************************
	// COMPUTE THE TOTAL NUMBER OF POINTS IN INITIAL POPUALTION
	// *******************************************************
	int npt = m_ngs * m_npg;
	int ngs1 = m_ngs;
	int npt1 = npt;

	// *******************************************************
	// open output file
	// *******************************************************
//	FILE *ipr;
//	if((ipr  = fopen(strOutFile, "at")) == NULL){
//		printf("Output file could not opened.\n");
//		bIsOK = false;
//		return bIsOK;
//	}

//	fprintf(ipr, " \n\n");
//	fprintf(ipr, "  ==================================================\n");
//	fprintf(ipr, "  ENTER THE SHUFFLED COMPLEX EVOLUTION GLOBAL SEARCH\n");
//	fprintf(ipr, "  ==================================================\n");
//	fprintf(ipr, " \n\n");
//	printf(" ***  Evolution Loop Number %d\n", nloop);

	// *******************************************************
	// COMPUTE THE BOUND FOR PARAMETERS BEING OPTIMIZED
	// *******************************************************
	for(j = 0; j < m_nopt; j++){
		bound.push_back(m_aryHighPrm[j] - m_aryLowPrm[j]);
		unit[j] = 1.0;
	}

	// *******************************************************
	// COMPUTE THE FUNCTION VALUE OF THE INITIAL POINT
	// *******************************************************
	double fa = functn(m_nopt, m_npar, m_iobj, cModel.m_aryObsData, 
		m_aryIniPrm, cModel.m_aryInputData, m_aryIniPrm, cModel.m_aryLocationOfOptParameter, model_func);

	// PRINT THE INITIAL POINT AND ITS CRITERION VALUE
//	fprintf(ipr, "*** PRINT THE INITIAL POINT AND ITS CRITERION VALUE ***\n");
//	fprintf(ipr, "\n");

//	fprintf(ipr, " CRITERION");
//	for(j = 0; j < m_nopt; j++){
//		sprintf(strLine, "X%d", j + 1);
//		fprintf(ipr, "%10s", strLine);
//	}
//	fprintf(ipr, "\n");

//	fprintf(ipr, "%9.2lf", fa);			// write initial value of objective function
//	for(j = 0; j < m_nopt; j++){
//		fprintf(ipr, "%10.2lf", m_aryIniPrm[j]);
//	}
//	fprintf(ipr, " \n");
//	fprintf(ipr, " \n\n");	
 
	// *******************************************************
	// GENERATE AN INITIAL SET OF npt1 POINTS IN THE PARAMETER SPACE
	// *******************************************************
	// IF iniflg IS EQUAL TO 1, SET x(1,.) TO INITIAL POINT a(.) 
	// パラメタに初期値を利用する場合
	if (m_iniflg == 1){
		for(j = 0; j < m_nopt; j++){
			x[0][j] = m_aryIniPrm[j];
		}
		xf[0] = fa;
	}
	// ELSE, GENERATE A POINT RANDOMLY AND SET IT EQUAL TO x(1,.)
	else{
		get_point(m_nopt,1,iseed1,xx, m_aryLowPrm, m_aryHighPrm, unit, m_aryLowPrm);
		for(j = 0; j < m_nopt; j++){
			x[0][j] = xx[j];
		}
		xf[0] = functn(m_nopt, m_npar, m_iobj, cModel.m_aryObsData, xx, 
			cModel.m_aryInputData, m_aryIniPrm, cModel.m_aryLocationOfOptParameter, model_func);
	}
	int icall = 1;
	if (icall >= m_maxn){
		iTerminatedError = 9000;
	}

	// GENERATE npt1-1 RANDOM POINTS DISTRIBUTED UNIFORMLY IN THE PARAMETER
	// SPACE, AND COMPUTE THE CORRESPONDING FUNCTION VALUES
    for(i = 1; i < npt1; i++){
		get_point(m_nopt, 1, iseed1, xx, m_aryLowPrm, m_aryHighPrm, unit, m_aryLowPrm);
		for(int j = 0; j < m_nopt; j++){
			x[i][j] = xx[j];
		}
        xf[i] = functn(m_nopt, m_npar, m_iobj, cModel.m_aryObsData, 
			xx, cModel.m_aryInputData, m_aryIniPrm, cModel.m_aryLocationOfOptParameter, model_func);
		// printf("xf[%d] = %lf\n", i, xf[i] );
        icall = icall + 1;
        if (icall >= m_maxn){
          npt1 = i;
		  break;
        }
	}

	// ARRANGE THE POINTS IN ORDER OF INCREASING FUNCTION VALUE
	sort_sce(npt1, m_nopt, x, xf);

	// RECORD THE BEST AND WORST POINTS
	for( j = 0; j < m_nopt; j++){
		bestx[j] = x[0][j];
		worstx[j] = x[npt1 - 1][j];
	}
	dBestFunc  = xf[0];
	dWorstFunc = xf[npt1 - 1];

	// COMPUTE THE PARAMETER RANGE FOR THE INITIAL POPULATION
//	parstt(npt1, m_nopt, x, xnstd, bound, gnrng, ipcnvg);

	// PRINT THE RESULTS FOR THE INITIAL POPULATION
//	fprintf(ipr, " *** PRINT THE RESULTS OF THE SCE SEARCH *** \n");
//	fprintf(ipr, "\n");

//	fprintf(ipr, " LOOP TRIALS COMPLXS  BEST F   WORST F   PAR RNG ");
//	for(j = 0; j < m_nopt; j++){
//		sprintf(strLine, "X%d", j + 1);
//		fprintf(ipr, "%10s", strLine);
//	}
//	fprintf(ipr, "\n");

//	fprintf(ipr, "%5d %5d %5d %10.3lf %10.3lf %10.3lf", nloop, icall, ngs1, dBestFunc, dWorstFunc, gnrng);
//	for(j = 0; j < m_nopt; j++){
//		fprintf(ipr, "%10.3lf", bestx[j]);
//	}
//	fprintf(ipr, "\n");

//	if (m_iprint == 1){
//		fprintf(ipr, " POPULATION AT LOOP %3d\n", nloop);
//		fprintf(ipr, " ----------------------\n");
//		for(i = 0; i < npt; i++){
//			fprintf(ipr, "               %10.3lf                    ", xf[i]);
//			for(j = 0; j < m_nopt; i++){
//				fprintf(ipr, "%10.3", x[i][j]);
//			}
//			fprintf(ipr, "\n");
//		}
//	}

//	if (icall >= m_maxn){
//		iTerminatedError = 9000;
//	}
//	if (ipcnvg >= 1){
//		iTerminatedError = 9200;
//	}

	// *******************************************************
	// BEGIN THE MAIN LOOP ----------------
	// *******************************************************
	while(iTerminatedError == 0){
		nloop = nloop + 1;
		printf(" ***  Evolution Loop Number %d\n", nloop);


		// *******************************************************
		// BEGIN LOOP ON COMPLEXES
		// *******************************************************
		for(igs = 0; igs < ngs1; igs++){

			// ASSIGN POINTS INTO COMPLEXES
			for( k1 = 0; k1 < m_npg; k1++){
				k2 = k1 * ngs1 + igs;
				for( j = 0; j < m_nopt; j++){
					cx[k1][j] = x[k2][j];
				}
				cf[k1] = xf[k2];
			}

			// *******************************************************
			// BEGIN INNER LOOP - RANDOM SELECTION OF SUB-COMPLEXES ---------------
			// *******************************************************
			for(loop = 0; loop < m_nspl; loop++){

				// CHOOSE A SUB-COMPLEX (nps points) ACCORDING TO A LINEAR
				// PROBABILITY DISTRIBUTION
				if (m_nps == m_npg){
					for(k = 0; k < m_nps; k++){
						lcs[k] = k;
					}
				}
				else{
					rand = rand_normal();
					lcs[0] = 1 + static_cast<int>(m_npg + 0.5 - sqrt( pow((m_npg+.5), 2) - m_npg * (m_npg+1) * rand ));
					for( k = 1; k < m_nps; k++){
						bool bLoop = true; 
						while(bLoop){
							rand = rand_normal();
							lpos = 1 + static_cast<int>(m_npg + 0.5 - sqrt(pow((m_npg+.5), 2) - m_npg * (m_npg+1) * rand ));
							for( k1 = 0; k1 < k-1; k1++){
								if (lpos == lcs[k1]){
									break;
								}
							}
							bLoop = false;
						}
						lcs[k] = lpos;
					}

					// ARRANGE THE SUB-COMPLEX IN ORDER OF INCEASING FUNCTION VALUE
					sort_sce(m_nps, lcs);
				}

				// CREATE THE SUB-COMPLEX ARRAYS
				for( k = 0; k < m_nps; k++){
					for( j = 0; j < m_nopt; j++){
						s[k][j] = cx[ lcs[k] ][j];
					}
					sf[k] = cf[ lcs[k] ];
				}

				// USE THE SUB-COMPLEX TO GENERATE NEW POINT(S)
				// cce(m_nopt, m_nps, s, sf, m_aryLowPrm, m_aryHighPrm, xnstd, icall, m_maxn, iseed1);
				cce( s, sf, xnstd, icall, iseed1, cModel, model_func);

				// IF THE SUB-COMPLEX IS ACCEPTED, REPLACE THE NEW SUB-COMPLEX
				// INTO THE COMPLEX
				for( k = 0; k < m_nps; k++){
					for( j = 0; j < m_nopt; j++){
						cx[ lcs[k] ][j] = s[k][j];
					}
					cf[lcs[k] ] = sf[k];
				}

				// SORT THE POINTS
				sort_sce(m_npg, m_nopt, cx, cf);

				// IF MAXIMUM NUMBER OF RUNS EXCEEDED, BREAK OUT OF THE LOOP
				if (icall >= m_maxn) break;

				// END OF INNER LOOP ------------
			}

			// REPLACE THE NEW COMPLEX INTO ORIGINAL ARRAY x(.,.)
			for( k1 = 0; k1 < m_npg; k1++){
				k2 = k1 * ngs1 + igs;
				for( j = 0; j < m_nopt; j++){
					x[k2][j] = cx[k1][j];
				}
				xf[k2] = cf[k1];
			}
			if (icall >= m_maxn) break;

			// END LOOP ON COMPLEXES
		}

		// RE-SORT THE POINTS
		sort_sce(npt1, m_nopt, x, xf);

		// RECORD THE BEST AND WORST POINTS
		for( j = 0; j < m_nopt; j++){
			bestx[j] = x[0][j];
			worstx[j] = x[npt1][j];
		}
		dBestFunc  = xf[0];
		dWorstFunc = xf[npt1];

		// TEST THE POPULATION FOR PARAMETER CONVERGENCE
		parstt(npt1, m_nopt, x, xnstd, bound, gnrng, ipcnvg);

		// PRINT THE RESULTS FOR CURRENT POPULATION
//		if ( (nloop % 5) == 0){
//			fprintf(ipr, " LOOP TRIALS COMPLXS  BEST F   WORST F   PAR RNG ");
//			for(j = 0; j < m_nopt; j++){
//				sprintf(strLine, "X%d", j + 1);
//				fprintf(ipr, "%10s", strLine);
//			}
//			fprintf(ipr, "\n");
//		}

//		fprintf(ipr, "%5d %5d %5d %10.3lf %10.3lf %10.3lf", nloop, icall, ngs1, dBestFunc, dWorstFunc, gnrng);
//		for(j = 0; j < m_nopt; j++){
//			fprintf(ipr, "%10.3lf", bestx[j]);
//		}
//		fprintf(ipr, "\n");

		// TEST IF MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED
		if (icall >= m_maxn){
			iTerminatedError = 9000;
			break;
		}

		// COMPUTE THE COUNT ON SUCCESSIVE LOOPS W/O FUNCTION IMPROVEMENT
		criter[TWENTY - 1] = dBestFunc;
		if (nloop >= (m_kstop + 1)){
			double denomi = fabs(criter[(TWENTY - 1) - m_kstop] + criter[19]) / 2.0;
			double timeou = fabs(criter[(TWENTY - 1) - m_kstop] - criter[19]) / denomi;

			if (timeou < m_dEpsPercent){
				iTerminatedError = 9100;
				break;
			}
		}
		for( l = 0; l < (TWENTY - 1); l++){
			criter[l] = criter[l+1];
		}

		// IF POPULATION IS CONVERGED INTO A SUFFICIENTLY SMALL SPACE
		if (ipcnvg == 1){
			iTerminatedError = 9200;
			break;
		}

		// NONE OF THE STOPPING CRITERIA IS SATISFIED, CONTINUE SEARCH

		// CHECK FOR COMPLEX NUMBER REDUCTION
		if (ngs1 > m_mings){
			ngs2 = ngs1;
			ngs1 = ngs1 - 1;
			npt1 = ngs1 * m_npg;
			comp(m_nopt, npt1, ngs1, ngs2, m_npg, x, xf, cx, cf);
		}

		// END OF MAIN LOOP -----------
	}

	// *******************************************************
	// SEARCH TERMINATED
	// *******************************************************
//	if(iTerminatedError != 0){
//		if(iTerminatedError == 9000){
//			fprintf(ipr, "\n\n");
//			fprintf(ipr, " *** OPTIMIZATION SEARCH TERMINATED BECAUSE THE LIMIT ON THE MAXIMUM\n");
//			fprintf(ipr, "     NUMBER OF TRIALS %5d EXCEEDED.  SEARCH WAS STOPPED AT\n", m_maxn);
//			fprintf(ipr, "     SUB-COMPLEX %3d OF COMPLEX %3d IN SHUFFLING LOOP %3d ***\n", loop, igs, nloop);
//		}
//		if(iTerminatedError == 9100){
//			fprintf(ipr, "\n\n");
//			fprintf(ipr, " *** OPTIMIZATION TERMINATED BECAUSE THE CRITERION  VALUE HAS NOT CHANGED\n");
//			fprintf(ipr, "     %5.2lf PERCENT IN %3d SHUFFLING LOOPS ***\n", m_dEpsPercent * 100.0, m_kstop);
//		}
//		if(iTerminatedError == 9200){
//			fprintf(ipr, "\n\n");
//			fprintf(ipr, " *** OPTIMIZATION TERMINATED BECAUSE THE POPULATION HAS CONVERGED INTO '\n");
//			fprintf(ipr, "     %5.2lf PERCENT OF THE FEASIBLE SPACE ***'\n", gnrng * 100.0);
//		}
//	}

	// *******************************************************
	// PRINT THE FINAL PARAMETER ESTIMATE AND ITS FUNCTION VALUE
	// *******************************************************
//	fprintf(ipr, "\n\n");
//	fprintf(ipr, " *** PRINT THE FINAL PARAMETER ESTIMATE AND ITS CRITERION VALUE ***'\n");
//	fprintf(ipr, "\n\n");
//
//	fprintf(ipr, " CRITERION");
//	for(j = 0; j < m_nopt; j++){
//		sprintf(strLine, "X%d", j + 1);
//		fprintf(ipr, "%10s", strLine);
//	}
//	fprintf(ipr, "\n");
//	fprintf(ipr, "------------------------------------------------------------\n");

//	fprintf(ipr, "%10.3lf", dBestFunc);
//	for(j = 0; j < m_nopt; j++){
//		sprintf(strLine, "%10.3lf", bestx[j]);
//		fprintf(ipr, "%s", strLine);
//	}
//	fprintf(ipr, "\n");

	// close file pointer
//	fclose(ipr);

	// set best parameter
	m_aryBestPrm.clear();
	for(i = 0; i < bestx.size(); i++){
		m_aryBestPrm.push_back(bestx[i]);
	}
	m_dBestFunc = dBestFunc;

//	cModel.m_aryBestSimulationData.clear();
//	cModel.m_aryBestSimulationData.assign( m_iRecordNum );
//	model_func(bestx, cModel.m_aryInputData, cModel.m_aryBestSimulationData);

	// *******************************************************
	return bIsOK;
}

void CSCE_UA::cce(std::vector< std::vector<double> >& s, std::vector<double>& sf, 
		 std::vector<double>& xnstd, int icall, int iseed, model_input& cModel,
		 void (*model_func) (std::vector<double>& , const std::vector< std::vector<double> >& ,  std::vector< std::vector<double> >&))
{
	// ALGORITHM GENERATE A NEW POINT(S) FROM A SUB-COMPLEX using competitive complex evolution
	// 
	// SUB-COMPLEX VARIABLES
	// parameter (c1=0.8,c2=0.4)

	// LIST OF LOCAL VARIABLES
	// sb(.) = the best point of the simplex
	// sw(.) = the worst point of the simplex
	// w2(.) = the second worst point of the simplex
	// fw = function value of the worst point
	// ce(.) = the centroid of the simplex excluding wo
	// snew(.) = new point generated from the simplex
	// iviol = flag indicating if constraints are violated
	//      = 1 , yes
	//     = 0 , no

	// EQUIVALENCE OF VARIABLES FOR READABILTY OF CODE
	const int n = m_nps;
	const int m = m_nopt;
	const double alpha = 1.0;
	const double beta = 0.5;

	vector<double> sw, sb, ce, snew;
	double fw;		// function value of the worst point
	double fnew;
	int i, j, ibound;

	sb.assign(m_nopt);
	sw.assign(m_nopt);
	ce.assign(m_nopt);
	snew.assign(m_nopt);

	// IDENTIFY THE WORST POINT wo OF THE SUB-COMPLEX s
	// COMPUTE THE CENTROID ce OF THE REMAINING POINTS
	// COMPUTE step, THE VECTOR BETWEEN wo AND ce
	// IDENTIFY THE WORST FUNCTION VALUE fw
	for (j = 0; j < m; j++){
		sb[j] = s[0][j];
		sw[j] = s[n-1][j];
		ce[j] = 0.0;
		for (i = 0; i < n - 1; i++){
			ce[j] = ce[j] + s[i][j];
		}
		ce[j] = ce[j] / static_cast<double>(n - 1);
	}
	fw = sf[n - 1];

	// COMPUTE THE NEW POINT snew

	// FIRST TRY A REFLECTION STEP
	for(j = 0; j < m; j++){
		snew[j] = ce[j] + alpha * (ce[j] - sw[j]);
	}

	// CHECK IF snew SATISFIES ALL CONSTRAINTS
	chkcst(m_nopt, snew, m_aryLowPrm, m_aryHighPrm, ibound);

		  
	// snew IS OUTSIDE THE BOUND,
	// CHOOSE A POINT AT RANDOM WITHIN FEASIBLE REGION ACCORDING TO
	// A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEX
	//  AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
	if (ibound >= 1) get_point(m_nopt, 2, iseed, snew, m_aryLowPrm, m_aryHighPrm, xnstd, sb);


	// COMPUTE THE FUNCTION VALUE AT snew
	fnew = functn(m_nopt, m_npar, m_iobj, cModel.m_aryObsData, 
		snew, cModel.m_aryInputData, m_aryIniPrm, cModel.m_aryLocationOfOptParameter, model_func);
	icall = icall + 1;
 
	// COMPARE fnew WITH THE WORST FUNCTION VALUE fw

	// fnew IS LESS THAN fw, ACCEPT THE NEW POINT snew AND RETURN
	if (fnew > fw){
		if (icall >= m_maxn){
			return;
		}

		// fnew IS GREATER THAN fw, SO TRY A CONTRACTION STEP
		for( j = 0; j < m; j++){
			snew[j] = ce[j] - beta * (ce[j] - sw[j]);
		}

		// COMPUTE THE FUNCTION VALUE OF THE CONTRACTED POINT
		fnew = functn(m_nopt, m_npar, m_iobj, cModel.m_aryObsData, 
			snew, cModel.m_aryInputData, m_aryIniPrm, cModel.m_aryLocationOfOptParameter, model_func);
		icall = icall + 1;

		// COMPARE fnew TO THE WORST VALUE fw
		// IF fnew IS LESS THAN OR EQUAL TO fw, THEN ACCEPT THE POINT AND RETURN
		if (fnew > fw){
			if (icall >= m_maxn) return;

			// IF BOTH REFLECTION AND CONTRACTION FAIL, CHOOSE ANOTHER POINT
			// ACCORDING TO A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEX
			// AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
			get_point(m_nopt, 2 , iseed, snew, m_aryLowPrm ,m_aryHighPrm, xnstd, sb);

			// COMPUTE THE FUNCTION VALUE AT THE RANDOM POINT
			fnew = functn(m_nopt, m_npar, m_iobj, cModel.m_aryObsData, 
				snew, cModel.m_aryInputData, m_aryIniPrm, cModel.m_aryLocationOfOptParameter, model_func);
			icall = icall + 1;
		}
	}

	// REPLACE THE WORST POINT BY THE NEW POINT
	for( j = 0; j < m; j++){
		s[n - 1][j] = snew[j];
	}
	sf[n - 1] = fnew;

	// END OF SUBROUTINE CCE
}

double CSCE_UA::functn(const int& nopt, const int& npar, const int& iobj,
			  std::vector< std::vector<double> >& obsq, std::vector<double>& x, std::vector< std::vector<double> >& p, std::vector<double>& pt, std::vector<int>& loc,
			  void (*model_func) (std::vector<double>& , const std::vector< std::vector<double> >& ,  std::vector< std::vector<double> >&) )
{
	double functn = -9999.0;
	int i, j, k, ndata;
	int iCount;
	vector<double> par;
	vector< vector<double> > q;

	if(obsq[0].size() != 0){
		ndata = obsq[0].size();
	}
	else{
		return functn;
	}

	// q.assign(obsq.size());
	q  = vector< vector<double> >(obsq.size(), vector<double>(obsq[0].size(), 0) );

	// *******************************************************
	// ASSIGN x TO APPROPRIATE PARAMETERS IN A HYDROLOGIC MODEL
	// *******************************************************
	for( j = 0; j < npar; j++){
		par.push_back(pt[j]);
	}
	for( k = 0; k < nopt; k++){
		par[ loc[k]-1 ] = x[k];
	}

	// *******************************************************
	// COMPUTE THE STREAMFLOW FOR PARAMETER SET par
	// *******************************************************
	model_func(par, p, q);

	// *******************************************************
	// COMPUTE THE SIMPLE LEAST SQUARE OBJECTIVE FUNCTION VALUE
	// *******************************************************
	if (iobj == 0){
		functn = 0.0;
		for( i = 0; i < ndata; i++){
			if(obsq[0][i] <= int(m_dErrorValue) || q[m_iOptimizedPrm][i] <= int(m_dErrorValue)){
				continue;
			}
			functn = functn + pow((obsq[m_iOptimizedPrm][i] - q[m_iOptimizedPrm][i]), 2);
		}
	}

	// *******************************************************
	// COMPUTE THE RMSE VALUE using single constraint
	// *******************************************************
	if (iobj == 1){
		functn = 0.0;
		iCount = 0;
		for(i = 0; i < ndata; i++){
			if(obsq[m_iOptimizedPrm][i] <= int(m_dErrorValue) || q[m_iOptimizedPrm][i] <= int(m_dErrorValue)){
				continue;
			}
			//if( fabs(obsq[m_iOptimizedPrm][i] - q[m_iOptimizedPrm][i]) > fabs(obsq[m_iOptimizedPrm][i]) * 0.15 ){
			//	functn = functn + pow( (obsq[m_iOptimizedPrm][i] - q[m_iOptimizedPrm][i]), 2);
			//}
			functn = functn + pow( (obsq[m_iOptimizedPrm][i] - q[m_iOptimizedPrm][i]), 2);
			iCount++;
		}
		if(iCount <= 1){
			functn = - m_dErrorValue;
		}
		else{
			functn = sqrt(functn / iCount);
		}
	}

	// *******************************************************
	// COMPUTE THE RMSE VALUE using two or more constraints
	// *******************************************************
	if(iobj > 1){
		functn = 0.0;
		iCount = 0;
		vector<double> aryAverage_obs, aryTxx, aryStdev_obs;
		aryAverage_obs.assign(obsq.size(), 0.0);
		aryStdev_obs.assign(obsq.size(), 0.0);
		aryTxx.assign(obsq.size(), 0.0);

		bool bError = false;
		for(i = 0; i < ndata; i++){
			bError = false;
			for(int j = 0; j < obsq.size(); j++){
				if(obsq[j][i] == m_dErrorValue ||  q[j][i] == m_dErrorValue){
					bError = true;
				}
			}
			if(bError == false){
				for(j = 0; j < obsq.size(); j++){
					aryAverage_obs[j] += obsq[j][i];
					aryTxx[j] += obsq[j][i] * obsq[j][i];
					iCount++;
				}
			}
		}
		if(iCount <= 1){
			functn = - m_dErrorValue;
			return functn;
		}
		for(j = 0; j < obsq.size(); j++){
			aryStdev_obs[j] = sqrt((aryTxx[j] - aryAverage_obs[j] * aryAverage_obs[j] / iCount) / iCount);
			aryAverage_obs[j] /= iCount;
		}

		vector<double> aryFunc;
		aryFunc.assign(obsq.size(), 0.0);
		for(i = 0; i < ndata; i++){
			bError = false;
			for(int j = 0; j < obsq.size(); j++){
				if(obsq[j][i] == m_dErrorValue ||  q[j][i] == m_dErrorValue){
					bError = true;
				}
			}
			if(bError) continue;

			for(j = 0; j < obsq.size(); j++){
				aryFunc[j] = functn + pow( (obsq[j][i] - q[j][i]), 2);
			}
		}
		for(j = 0; j < obsq.size(); j++){
			functn += sqrt(aryFunc[j] / iCount) / aryAverage_obs[j];
		}
	}

	return functn;
}

void CSCE_UA::create_random_seed_defualt(int nrun, bool bUseTime)
{
	// initilizaton jseed	

	int iTime = 0;
	if(bUseTime){
		iTime = ((unsigned) time(NULL) ) % INT_MAX;
	}

	for(int i = 0; i < nrun; i++){
		m_jseed.push_back(i + 1 + iTime);
	}
}

void CSCE_UA::set_default_sce_parameters(int nopt, int ngs)
{	
	m_npg = 2 * nopt + 1;
	m_nps = nopt + 1;
	m_nspl = m_npg;
	m_mings = ngs;
	m_iniflg = 0;
	m_iprint = 0;
}

double CSCE_UA::max(double iA, double iB)
{
	if(iA >= iB){
		return iA;
	}
	else{
		return iB;
	}
}

double CSCE_UA::min(double iA, double iB)
{
	if(iA >= iB){
		return iB;
	}
	else{
		return iA;
	}
}

double CSCE_UA::rand_normal()
{
	double dRand1;

	dRand1 = ((double)rand()+1.0)/((double)RAND_MAX+2.0);
	return dRand1;
}

double CSCE_UA::rand_Gaussian_dev(int& idum)
{
	// THIS SUBROUTINE IS FROM "NUMERICAL RECIPES" BY PRESS ET AL.
	double dGasdev;
//	double fac;
	double gset = -1.0;

	// http://www.cs.miyazaki-u.ac.jp/~date/lectures/random/
	static int iset = 0;
	static double v1, v2, r;
	if (iset == 0){
		iset = 1;
		 do{
			v1 = (2. * rand_normal()) - 1.0;
			v2 = (2. * rand_normal()) - 1.0;
			r = v1 * v1 + v2 * v2;
		}while(r >=1. || r == 0.0);
		r = sqrt(- 2.0 * log(r) / r);
		dGasdev = v1 * r;
	}
	else{
		iset = 0;
		dGasdev = v2 * r;
	}

	return dGasdev;
}

void CSCE_UA::chkcst(int nopt, std::vector<double>& x, std::vector<double>& bl, std::vector<double>& bu, int& ibound)
{
	// This subroutine check if the trial point satisfies all constraints.
	//
	// ibound - violation indicator
	//       = -1 initial value
	//       = 0  no violation
	//        = 1  violation
	// nopt = number of optimizing variables
	// ii = the ii'th variable of the arrays x, bl, and bu
	
	ibound = -1;

	// Check if explicit constraints are violated
	for(int ii = 0; ii < nopt; ii++){
		if(x[ii] == bl[ii]){
			continue;
		}
		if (x[ii] < bl[ii] || x[ii] > bu[ii]){
			ibound = 1;
			return;
		}
	}
	if (nopt == 1){
		ibound = 0;
		return;
	}
}

void CSCE_UA::chkcst(int nopt, double& x, double& bl, double& bu, int& ibound)
{
	// This subroutine check if the trial point satisfies all constraints.
	//
	// ibound - violation indicator
	//       = -1 initial value
	//       = 0  no violation
	//        = 1  violation
	// nopt = number of optimizing variables
	// ii = the ii'th variable of the arrays x, bl, and bu
	
	ibound = -1;

	// Check if explicit constraints are violated
	if(x == bl){
		return;
	}
	if (x < bl || x > bu){
		ibound = 1;
		return;
	}
	if (nopt == 1){
		ibound = 0;
		return;
	}
}

bool CSCE_UA::compare(const sort_set &s1, const sort_set &s2) 
 {
	return s1.m_value < s2.m_value;
 }


void CSCE_UA::sort_sce(int n, int m, std::vector< std::vector<double> >& rb, std::vector<double>& ra)
{
	int i, j;
//	vector<int> aryIndex;
//	vector<double> arySorted;
//	vector< vector<double> > arySortedMatrix;

	vector<sort_set> arySort;
	arySort.assign(n);

	for(i = 0; i < n; i++){
		arySort[i].m_index = i;
		arySort[i].m_value = ra[i];

		for(j = 0; j < m; j++){
			arySort[i].m_aryParamters.push_back(rb[i][j]);
		}
	}

	sort(arySort.begin(), arySort.end(), compare );

	for(i = 0; i < n; i++){
		ra[i] = arySort[i].m_value;

		for(j = 0; j < m; j++){
			rb[i][j] = arySort[i].m_aryParamters[j];
		}
	}

	return;
}

void CSCE_UA::sort_sce(int n, std::vector<int>& ra)
{
   int i, j, temp;

    for (i = 0; i < n - 1; i++) {
        for (j = n - 1; j > i; j--) {
            if (ra[j - 1] > ra[j]) {  /* 前の要素の方が大きかったら */
                temp = ra[j];        /* 交換する */
                ra[j] = ra[j - 1];
                ra[j - 1]= temp;
            }
        }
    }
}

void CSCE_UA::get_point(int& nopt, int idist, int& iseed, std::vector<double>& x, std::vector<double>& bl, std::vector<double>& bu, std::vector<double>& std, std::vector<double>& xi)
{
	// This subroutine generates a new point within feasible region
	// 
	// x(.) = new point
	// xi(.) = focal point
	// bl(.) = lower bound
	// bu(.) = upper bound
	// std(.) = standard deviation of probability distribution
	// idist = probability flag
	//      = 1 - uniform distribution
	//      = 2 - Gaussian distribution
/*c
      implicit real*8 (a-h,o-z)
      dimension x(16),bl(16),bu(16),std(16),xi(16)
c*/
	double rand;
	int ibound = 1;

	do{
		for(int j = 0; j < nopt; j++){
			do{
				if (idist == 1) rand = rand_normal();
				if (idist == 2) rand = rand_Gaussian_dev(iseed);
				x[j] = xi[j] + std[j] * rand * (bu[j] - bl[j]);
				
				// Check explicit constraints
				chkcst(1, x[j], bl[j], bu[j], ibound);
			}while(ibound >= 1);
		}

		// Check implicit constraints
		chkcst(nopt, x, bl, bu, ibound);
	}while(ibound >= 1);
}

void CSCE_UA::parstt(int npt, int nopt, std::vector< std::vector<double> >& x, std::vector<double>& xnstd, std::vector<double>& bound, double& gnrng, int& ipcnvg)
{
	// SUBROUTINE CHECKING FOR PARAMETER CONVERGENCE
	const double delta = 1.0 * pow(10, -20);
	const double peps  = 1.0 * pow(10, -3);	
	int i;
	double gsum = 0.0, xsum1, xsum2;
	vector<double> xmax, xmin, xmean;

	// COMPUTE MAXIMUM, MINIMUM AND STANDARD DEVIATION OF PARAMETER VALUES
	xmax.assign(m_nopt);
	xmin.assign(m_nopt);
	xmean.assign(m_nopt);

	for(int k = 0; k < nopt; k++){
		xmax[k] = -1.0 * pow(10, 20);
		xmin[k] =  1.0 * pow(10, -20);
		xsum1 = 0.0;
		xsum2 = 0.0;
		for(i = 0; i < npt; i++){
			xmax[k] = max(x[i][k], xmax[k]);
			xmin[k] = min(x[i][k], xmin[k]);
			xsum1 = xsum1 + x[i][k];
			xsum2 = xsum2 + x[i][k] * x[i][k];
		}
		xmean[k] = xsum1 / static_cast<double>(npt);
		xnstd[k] = (xsum2 / static_cast<double>(npt) - xmean[k] * xmean[k] );
		if (xnstd[k] <= delta) xnstd[k] = delta;
		xnstd[k] = sqrt(xnstd[k]);
		xnstd[k] = xnstd[k] / bound[k];
		gsum = gsum + log( delta + (xmax[k] - xmin[k])/bound[k] );
	}
	gnrng = exp(gsum/static_cast<double>(nopt - 1));

	// CHECK IF NORMALIZED STANDARD DEVIATION OF PARAMETER IS <= eps
	ipcnvg = 0;
	if (gnrng <= peps){
		ipcnvg = 1;
	}
}

void CSCE_UA::comp(int n, int npt, int ngs1, int ngs2, int npg, std::vector< std::vector<double> >& a, std::vector<double>& af, std::vector< std::vector<double> >& b, std::vector<double>& bf)
{
	// THIS SUBROUTINE REDUCE INPUT MATRIX a(n,ngs2*npg) TO MATRIX
	// b(n,ngs1*npg) AND VECTOR af(ngs2*npg) TO VECTOR bf(ngs1*npg)
/*	implicit real*8 (a-h,o-z)
	dimension a(2000,16),af(2000),b(2000,16),bf(2000)*/

	int i, j, igs, ipg, k1, k2;

	for( igs = 0; igs < ngs1; igs++){
		for( ipg = 0; ipg < npg; ipg++){
			k1 = (ipg - 1) * ngs2 + igs;
			k2 = (ipg - 1) * ngs1 + igs;
			for( i = 0; i < n; i++){
				b[k2][i] = a[k1][i];
			}
			bf[k2] = af[k1];
		}
	}

	for( j = 0; j < npt; j++){
		for( i = 0; i < n; i++){
			a[j][i] = b[j][i];
		}
		af[j] = bf[j];
	}
}
