/* Copyright (C) 2000  Dale Schuurmans, Finnegan Southey                     */
/* This work is licensed under the Gnu General Public License (see gpl.txt). */

/* This is a floating-point implementation of SDF. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "satfront.h"
#include "readcnf.h"

#define swap(x,y) temp=x; x=y; y=temp
#define swapint(x,y) tempint=x; x=y; y=tempint

	/* When on, vars occuring in false clauses will be selected over those that
	 * do not, even if they offer a smaller gain in objective value. */
#define PRIORITIZE_FALSE_CLAUSE_VARS 1

	/* Sanity checks that determines if the same variable is flipped twice under
	 * trivial conditions. */
#define CHECK_OSCILLATION 0

	/* Turns on heap-based flip selection. */
	/* No advantage observed at n=3016 (bw_large.c). */
	/* Some advantage observed at n=6325 (bw_large.d). */
#define USE_HEAP 0

#if USE_HEAP
	#define SUMSCORE_X sumscore[varheapind[x]]
	#define SUMSCORE_XFLIP sumscore[varheapind[xflip]]
#else
	#define SUMSCORE_X sumscore[x]
	#define SUMSCORE_XFLIP sumscore[xflip]
#endif

typedef double REAL;

REAL temp;
REAL tempint;
int numsat, flips, r;

int *cumcount;
REAL *satdiff;
REAL *diffdiff;
REAL *satscore;
int *sat;
REAL *weight;
int *yfalse;
int *wherefalse;
int *floodtime;
REAL *falseadd;
REAL *satadd;
REAL *sat1sub;
REAL *satmsub;
REAL *sumscore;

int *xinfalse;
int *wherexinfalse;
int *countxinfalse;
int numxinfalse;

#if USE_HEAP
	int *scorevar;
	int *varheapind;
#endif

int *make;
int *brake;
REAL *rawsatadd;
REAL *rawsatmsub;
int *maxXs;
REAL *betas;
REAL *satmeans;
REAL *rhobeta;
REAL *cummean;

/*************** Incremental SHRUNK NON-BOOSTING FLOOD-SAT *******************/

void preSolve(SATParams *params)
	{
	int T, F, l, litScoreCeil;

	cumcount = malloc((Nmax+1) * sizeof(int));
	satdiff = malloc((Nmax+1) * sizeof(REAL));
	diffdiff = malloc((Nmax+1) * sizeof(REAL));
	satscore = malloc((Nmax+1) * sizeof(REAL));
	sat = malloc(m * sizeof(int));
	weight = malloc(m * sizeof(REAL));
	yfalse = malloc(m * sizeof(int));
	wherefalse = malloc(m * sizeof(int));
	floodtime = malloc(m * sizeof(int));
	falseadd = malloc((n+1) * sizeof(REAL));
	satadd = malloc((n+1) * sizeof(REAL));
	sat1sub = malloc((n+1) * sizeof(REAL));
	satmsub = malloc((n+1) * sizeof(REAL));
	sumscore = malloc((n+1) * sizeof(REAL));
	xinfalse = malloc(n * sizeof(int));
	wherexinfalse = malloc((n+1) * sizeof(int));
	countxinfalse = malloc((n+1) * sizeof(int));

	#if USE_HEAP
		scorevar = malloc((n+1) * sizeof(int));
		varheapind = malloc((n+1) * sizeof(int));
	#endif

	make = malloc((n+1) * sizeof(int));
	brake = malloc((n+1) * sizeof(int));
	rawsatadd = malloc((n+1) * sizeof(REAL));
	rawsatmsub = malloc((n+1) * sizeof(REAL));
	maxXs = malloc((n+1) * sizeof(int));

	T =  params->intParams[0];

	if (T == 0)
		T = 100000;

	F = T/2;

	litScoreCeil = params->intParams[4];
	if (litScoreCeil == 0)
		litScoreCeil = Nmax;
	else
		if (litScoreCeil > Nmax)
			litScoreCeil = Nmax;

	betas = malloc(F * sizeof(REAL));
	satmeans = malloc(F * sizeof(REAL));
	rhobeta = malloc(F * sizeof(REAL));
	cummean = malloc(F * sizeof(REAL));

	for (l = Nmax; l > litScoreCeil; l--) {
		satdiff[l] = 0;
		diffdiff[l] = 0;
	}
	satdiff[litScoreCeil] = 1.0;

	cumcount[Nmax] = lencount[Nmax];
	for (l = Nmax - 1; l >= 0; l--) 
		cumcount[l] = cumcount[l+1] + lencount[l];

	for (l = litScoreCeil - 1; l > 0; l--) {
		diffdiff[l] = satdiff[l+1] * cumcount[l+1];
		satdiff[l] = satdiff[l+1] + diffdiff[l];
	}
	for (l = 2; l <= litScoreCeil; l++) {
		diffdiff[l] /= satdiff[1];
		satdiff[l] /= satdiff[1];
	}

	satdiff[1] = 1.0;
	satscore[0] = 0.0;
	for (l = 1; l <= litScoreCeil; l++) 
		satscore[l] = satscore[l-1] + satdiff[l];
	for (l = litScoreCeil+1; l <= Nmax; l++) 
		satscore[l] = satscore[l-1];
	}

void runSolver(SATParams *params, SATResults *results)
	{
	int R;
	int T;
	int F;
	REAL delta;
	REAL rho;
	REAL oneminusrho;
	int minorCorrPeriod;
	int majorCorrPeriod;

	int t, f, x, y, i, j, xflip =-1 , u;
	int nfalse;
	int nmxs;
	int floods, uphist, tempint;
	REAL maxdiff;
	REAL total;
	REAL sattotal, falsetotal;
	REAL oldsattotal;
	REAL meansat, dsattotal;
	REAL alpha, minalpha, beta;

	#if CHECK_OSCILLATION
		int lastFlip = -1;
	#endif

	#include "mobility1.c"

	T =  params->intParams[0];
	R = params->intParams[1];
	minorCorrPeriod = params->intParams[2];
	majorCorrPeriod = params->intParams[3];
	delta = params->doubleParams[0];
	rho = params->doubleParams[1];

	if (T == 0)
		T = 100000;

	if (R == 0)
		R = 5;

	if (delta == 0)
		delta = 0.05/n;

	if (rho == 0)
		rho = 0.995;

	if (minorCorrPeriod == 0)
		minorCorrPeriod = 10;

	if (majorCorrPeriod == 0)
		majorCorrPeriod = 10;

	majorCorrPeriod *= minorCorrPeriod;

	oneminusrho = 1.0 - rho;

	F = T/2;

	flips = floods = 0;

	for (r = 0; r < R; r++) {  /********* r ***************/

		for (x = 1; x <= n; x++) {
			assign[x] = (rand()%100 < 50) * 2 - 1;
			countxinfalse[x] = 0;
		}

		flips++;

		numsat = nfalse = 0;
		falsetotal = 0.0;
		numxinfalse = 0;

		for (y = 0; y < m; y++) {
			sat[y] = floodtime[y] = 0;
			weight[y] = 1.0 / m;
			for (i = 0; i<nXs[y]; i++) 
				sat[y] += assign[Xs[y][i]] == Xsgns[y][i];
			if (sat[y] > 0) {
				numsat++;
			} else {
				for (i = 0; i<nXs[y]; i++) {
					x = Xs[y][i];
					if (++countxinfalse[x] == 1) {
						xinfalse[numxinfalse] = x;
						wherexinfalse[x] = numxinfalse;
						numxinfalse++;
					}
				}
				yfalse[nfalse] = y;
				wherefalse[y] = nfalse;
				nfalse++;
			}
		}

		f = uphist = 0;
		betas[0] = 1.0;
		rhobeta[0] = 1.0;
		falsetotal = (double)nfalse / m;
		sattotal = 1.0 - falsetotal;
		total = 1.0;
		satmeans[0] = 1.0 / m;
		cummean[0] = oneminusrho / m;

		for (x = 1; x <= n; x++) {
			falseadd[x] = satadd[x] = sat1sub[x] = satmsub[x] = 0.0;
			make[x] = brake[x] = 0;
			rawsatadd[x] = rawsatmsub[x] = 0.0;
			for (j = 0; j < nYs[x]; j++) {
				y = Ys[x][j];
				if (sat[y] == 0) {
					falseadd[x] += weight[y] * satdiff[1];
					make[x]++;
				} else if (assign[x] != Ysgns[x][j]) {
					satadd[x] += weight[y] * satdiff[sat[y]+1];
					rawsatadd[x] += satdiff[sat[y]+1];
				} else if (sat[y] == 1) {
					sat1sub[x] += weight[y] * satdiff[1];
					brake[x]++;
				} else {
					satmsub[x] += weight[y] * satdiff[sat[y]];
					rawsatmsub[x] += satdiff[sat[y]];
				}
			}

			sumscore[x] = (falseadd[x] + satadd[x]) - (sat1sub[x] + satmsub[x]);
			#if USE_HEAP
				scorevar[x] = x;
				varheapind[x] = x;
			#endif
		}

	
		for (t = 1; t < T; t++) { /********** t **************/

			if (numsat == m) break;

			#include "unsat1.c"
			#include "mobility2.c"
			#include "assigndump.c"


			#if USE_HEAP
				buildHeap(sumscore, scorevar, varheapind, n);

				maxdiff = sumscore[1];
				maxXs[0] = scorevar[1];
				nmxs = 1;
			#else
				maxdiff = -1e50;
				nmxs = 0;

				#if PRIORITIZE_FALSE_CLAUSE_VARS
					for (j = 0; j < numxinfalse; j++) {
						x = xinfalse[j];
						if (sumscore[x] > maxdiff) {
							maxdiff = sumscore[x];
							maxXs[0] = x;
							nmxs = 1;
						} else if (sumscore[x] == maxdiff)
							maxXs[nmxs++] = x;
					}
					if (maxdiff <= 0) {
						for (x = 1; x <= n; x++)
							if (sumscore[x] > maxdiff) {
								maxdiff = sumscore[x];
								maxXs[0] = x;
								nmxs = 1;
							} else if (sumscore[x] == maxdiff)
								maxXs[nmxs++] = x;
					}
				#else			
					for (x = 1; x <= n; x++)
						if (sumscore[x] > maxdiff) {
							maxdiff = sumscore[x];
							maxXs[0] = x;
							nmxs = 1;
						} else if (sumscore[x] == maxdiff)
							maxXs[nmxs++] = x;
				#endif
			#endif

			if (maxdiff <= 0.0) {     /*********** FLOOD **********/
				REAL betarho, meansatoneminusrho;

				if ((f+1) % minorCorrPeriod == 0) {
					sattotal = falsetotal = 0.0;
					for (y = 0; y < m; y++) {
						if (floodtime[y] < f) {
							if (floodtime[y] < uphist) {
								for (u = uphist - 1; u >= floodtime[y]; u--) {
									rhobeta[u] = rhobeta[u+1] * betas[u+1] * rho;
									cummean[u] = cummean[u+1] + 
										oneminusrho*rhobeta[u]*satmeans[u];
								}
								uphist = floodtime[y];
							}
							weight[y] = weight[y]*rhobeta[floodtime[y]]
								+ cummean[floodtime[y]+1];
							floodtime[y] = f;
						}
						if (sat[y] == 0)
							falsetotal += weight[y];
						else
							sattotal += weight[y];
					}


					if ((f+1) % majorCorrPeriod == 0) {
					
						total = sattotal + falsetotal;

							/* Normalize weights. */
						for (y = 0; y < m; y++)
							weight[y] /= total;

						falsetotal /= total;
						sattotal = 1.0 - falsetotal;

							/* Recompute scores. */
						for (x = 1; x <= n; x++) {
							falseadd[x] = satadd[x] = sat1sub[x] = satmsub[x] = 0.0;
							rawsatadd[x] = rawsatmsub[x] = 0.0;
							for (j = 0; j < nYs[x]; j++) {
								y = Ys[x][j];
								if (sat[y] == 0) {
									falseadd[x] += weight[y] * satdiff[1];
								} else if (assign[x] != Ysgns[x][j]) {
									satadd[x] += weight[y] * satdiff[sat[y]+1];
									rawsatadd[x] += satdiff[sat[y]+1];
								} else if (sat[y] == 1) {
									sat1sub[x] += weight[y] * satdiff[1];
								} else {
									satmsub[x] += weight[y] * satdiff[sat[y]];
									rawsatmsub[x] += satdiff[sat[y]];
								}
							}
			
							SUMSCORE_X = (falseadd[x] + satadd[x]) - (sat1sub[x] + satmsub[x]);
						}
					}
				}

				floods++;
				f++;
				uphist = f;

				if (f >= F) {
					printf("Error: max number of floods exceeded\n");
					exit(1);
				}

				total = sattotal + falsetotal;
				dsattotal = delta * sattotal;
				minalpha = 1e20;
				nmxs = 0;

					/** MINALPHA SEARCH **/
				for (j = 0; j < numxinfalse; j++) {
					REAL score;

					x = xinfalse[j];

					score = (sat1sub[x]+satmsub[x]-satadd[x]);
					alpha = total * (score+dsattotal) / 
						( falsetotal*score + falseadd[x]*sattotal );

					if (alpha <= 1.0) {
						printf("Error: alpha is messed up, alpha=%g\n", alpha);
						exit(1);
					}

					if (alpha < minalpha) {
						minalpha = alpha;
						maxXs[0] = x;
						nmxs = 1;
					} else if (alpha == minalpha)
						maxXs[nmxs++] = x;
				} 

				falsetotal *= minalpha;

				for (j = 0; j < nfalse; j++) {
					y = yfalse[j];
					weight[y] *= minalpha;
					floodtime[y] = f;
				}

				oldsattotal = sattotal;
				sattotal = total - falsetotal;

				satmeans[f] = meansat = sattotal / numsat;

				betas[f] = beta = sattotal / oldsattotal;

				cummean[f] = oneminusrho * meansat;

				rhobeta[f] = 1.0;

				if (beta >= 1.0 || beta <= 0.0) {
					printf("Error: beta is messed up, beta=%g\n", beta);
					exit(1);
				}

				betarho = beta * rho;
				meansatoneminusrho = meansat * oneminusrho;
					
				for (x = 1; x <= n; x++) {
					if (falseadd[x] > 0.0)
						falseadd[x] *= minalpha;

					if (satadd[x] <= 0.0)
						satadd[x] = 0.0;
					else
						satadd[x] = betarho*satadd[x] + rawsatadd[x]*meansatoneminusrho;

					if (sat1sub[x] > 0.0)
						sat1sub[x] = betarho*sat1sub[x] + brake[x]*satdiff[1]*meansatoneminusrho;

					if (satmsub[x] <= 0.0)
						satmsub[x] = 0.0;
					else
						satmsub[x] = betarho*satmsub[x] + rawsatmsub[x]*meansatoneminusrho;

					SUMSCORE_X = (falseadd[x] + satadd[x]) - (sat1sub[x] + satmsub[x]);
				}

				#if CHECK_OSCILLATION
					/* If we have just flooded, we can allow a duplicate flip. */
					lastFlip = -1;
				#endif
			}                         /****************************/


			if (nmxs == 1)
				xflip = maxXs[0];
			else
				xflip = maxXs[rand()%nmxs];

			#if CHECK_OSCILLATION
				if (lastFlip == xflip) {
					printf("Error: Flipped same variable twice consecutively: var %d flip %d\n",
							lastFlip, flips);
					exit(1);
				}

				lastFlip = xflip;
			#endif

			#ifdef ADD_NOISE
				if(rand() % 100000 < ADDED_NOISE * 100000)
					xflip = rand() % n + 1;
			#endif

			flips++;
			assign[xflip] = -assign[xflip];
			swap(falseadd[xflip], sat1sub[xflip]);
			swap(satadd[xflip], satmsub[xflip]);
			swapint(make[xflip], brake[xflip]);
			swap(rawsatadd[xflip], rawsatmsub[xflip]);

			SUMSCORE_XFLIP = -SUMSCORE_XFLIP;

			for (j = 0; j < nYs[xflip]; j++) {

				y = Ys[xflip][j];

				if (floodtime[y] < f) {
					if (floodtime[y] < uphist) {
						for (u = uphist - 1; u >= floodtime[y]; u--) {
							rhobeta[u] = rhobeta[u+1] * betas[u+1] * rho;
							cummean[u] = cummean[u+1] + 
								oneminusrho * rhobeta[u] * satmeans[u];
						}
						uphist = floodtime[y];
					}

					weight[y] = weight[y] * rhobeta[floodtime[y]]
						+ cummean[floodtime[y]+1];
					floodtime[y] = f;
				}
				if (assign[xflip] != Ysgns[xflip][j]) {
					if (--sat[y] == 0) {
						numsat--;
						yfalse[nfalse] = y;
						wherefalse[y] = nfalse;
						nfalse++;
						falsetotal += weight[y];
						sattotal  = total - falsetotal;
						for (i = 0; i < nXs[y]; i++) {
							REAL t1, t2, t3;
			
							x = Xs[y][i];

							if (++countxinfalse[x] == 1) {
								xinfalse[numxinfalse] = x;
								wherexinfalse[x] = numxinfalse;
								numxinfalse++;
							}

							if (x  == xflip)
								continue;

							t1 = weight[y] * satdiff[1];
							falseadd[x] += t1;
							make[x]++;
							t2 = weight[y] * satdiff[2];
							t3 = satadd[x] - t2;
							if (t3 <= 0.0) {
								SUMSCORE_X += t1 - satadd[x];
								satadd[x] = rawsatadd[x] = 0.0;
							} else {
								SUMSCORE_X += t1 - t2;
								satadd[x] = t3;
								rawsatadd[x] -= satdiff[2];
							}
						}
					} else {
						for (i = 0; i < nXs[y]; i++) {
							if ((x = Xs[y][i]) == xflip) 
								continue;

							if (Xsgns[y][i] != assign[x]) {
								REAL t1;
								t1 = weight[y] * diffdiff[sat[y]+1];
								satadd[x] += t1;
								rawsatadd[x]+= diffdiff[sat[y]+1];
								SUMSCORE_X += t1;
							} else if (sat[y] == 1) {
								REAL t1, t2, t3;
								t1 = weight[y] * satdiff[1];
								sat1sub[x] += t1;
								brake[x]++;
								t2 = weight[y] * satdiff[2];
								t3 = satmsub[x] - t2;
								if (t3 <= 0.0) {
									SUMSCORE_X -= t1 - satmsub[x];
									satmsub[x] = rawsatmsub[x] = 0.0;
								}
								else {
									SUMSCORE_X -= t1 - t2;
									rawsatmsub[x] -= satdiff[2];
									satmsub[x] = t3;
								}
							} else {
								REAL t1;
								t1 = weight[y] * diffdiff[sat[y]];
								satmsub[x] += t1;
								rawsatmsub[x] += diffdiff[sat[y]];
								SUMSCORE_X -= t1;
							}
						}
					}
				} else {
					if (++sat[y] == 1) {
						numsat++;
						nfalse--;
						yfalse[wherefalse[y]] = yfalse[nfalse];
						wherefalse[yfalse[nfalse]] = wherefalse[y];
						falsetotal -= weight[y];
						sattotal = total - falsetotal;
						for (i = 0; i < nXs[y]; i++) {
							REAL t1;

							x = Xs[y][i];

							if (--countxinfalse[x] == 0) {
								int wheredrop, endx;
								numxinfalse--;
								endx = xinfalse[numxinfalse];
								wheredrop = wherexinfalse[x];
								xinfalse[wheredrop] = endx;
								wherexinfalse[endx] = wheredrop;
							}

							if (x == xflip)
								continue;

							t1 = weight[y] * satdiff[2];
							satadd[x] += t1;
							rawsatadd[x] += satdiff[2];		

							if (--make[x] == 0) {
								SUMSCORE_X += t1 - falseadd[x];
								falseadd[x] = 0.0;
							}
							else {
								REAL t2;

								t2 = weight[y] * satdiff[1];
								falseadd[x] -= t2;
								SUMSCORE_X += t1 - t2;
							}
						}
					} else {
						for (i = 0; i < nXs[y]; i++) {
							if ((x = Xs[y][i]) == xflip) 
								continue;
							if (Xsgns[y][i] != assign[x]) {
								REAL t1, t2;
								t1 = weight[y] * diffdiff[sat[y]];
								t2 = satadd[x] - t1;
								if (t2 <= 0.0) {
									SUMSCORE_X -= satadd[x];
									satadd[x] = rawsatadd[x] = 0.0;
								}
								else {
									SUMSCORE_X -= t1;
									satadd[x] = t2;
									rawsatadd[x] -= diffdiff[sat[y]];
								}
							} else if (sat[y] == 2) {
								REAL t1;

								t1 = weight[y] * satdiff[2];
								satmsub[x] += t1;
								rawsatmsub[x] += satdiff[2];

								if (--brake[x] == 0) {
									SUMSCORE_X -= (t1 - sat1sub[x]);
									sat1sub[x] = 0.0;
								}
								else {
									REAL t2;
									t2 = weight[y] * satdiff[1];
									sat1sub[x] -= t2;
									SUMSCORE_X -= (t1 - t2);
								}
							} else {
								REAL t1, t2;
								t1 = weight[y] * diffdiff[sat[y]-1];
								t2 = satmsub[x] - t1;
								if (t2 <= 0.0) {
									SUMSCORE_X -= satmsub[x];
									satmsub[x] = rawsatmsub[x] = 0.0;
								}
								else {
									SUMSCORE_X += t1;
									satmsub[x] = t2;
									rawsatmsub[x] -= diffdiff[sat[y]-1];
								}
							}
						}
					}
				}
			}
		}                         /********** t **************/
		if (numsat == m) break;
	}                          /********* r ***************/

	if (numsat != m) 
		r--;

	results->numFlips = flips;
	results->numReweights = floods;
	results->numRestarts = r + 1;
	results->numClauses = m;
	results->numSatisfied = numsat;
	results->isSatisfied = numsat == m;

	#include "mobility3.c"
	}

void postSolve(SATParams *params)
	{
	free(cumcount);
	free(satdiff);
	free(diffdiff);
	free(satscore);
	free(sat);
	free(weight);
	free(yfalse);
	free(wherefalse);
	free(floodtime);
	free(betas);
	free(satmeans);
	free(rhobeta);
	free(cummean);
	free(falseadd);
	free(satadd);
	free(sat1sub);
	free(satmsub);
	free(sumscore);
	free(xinfalse);
	free(countxinfalse);
	free(wherexinfalse);

	#if USE_HEAP
		free(scorevar);
		free(varheapind);
	#endif

	free(make);
	free(brake);
	free(rawsatadd);
	free(rawsatmsub);
	free(maxXs);
	}

SATParams *getSATParams()
	{
	SATParams *params;

	params = malloc(sizeof(SATParams));
	params->numIntParams = 5;
	params->intParams = malloc(params->numIntParams * sizeof(int));
	params->intParamNames = malloc(params->numIntParams * sizeof(char*));
	params->intParamDescs = malloc(params->numIntParams * sizeof(char*));
	params->intParamNames[0] = "-mf";
	params->intParamDescs[0] = "max flips before restarting";
	params->intParams[0] = 0;
	params->intParamNames[1] = "-mr";
	params->intParamDescs[1] = "max restarts before aborting";
	params->intParams[1] = 0;
	params->intParamNames[2] = "-pcp";
	params->intParamDescs[2] = "number of floods between partial corrections";
	params->intParams[2] = 0;
	params->intParamNames[3] = "-fcp";
	params->intParamDescs[3] = "number of partial corrections between full corrections";
	params->intParams[3] = 0;
	params->intParamNames[4] = "-lc";
	params->intParamDescs[4] = "number of literals in smoothing";
	params->intParams[4] = 0;


	params->numDoubleParams = 2;
	params->doubleParams = malloc(params->numDoubleParams * sizeof(double));
	params->doubleParamNames = malloc(params->numDoubleParams * sizeof(char*));
	params->doubleParamDescs = malloc(params->numDoubleParams * sizeof(char*));
	params->doubleParamNames[0] = "-del";
	params->doubleParamDescs[0] = "target post-flood gradient";
	params->doubleParams[0] = 0;
	params->doubleParamNames[1] = "-rho";
	params->doubleParamDescs[1] = "rate of weight shrinkage to mean for sat clauses";
	params->doubleParams[1] = 0;

	params->numStringParams = 0;

	return params;
	}

char *getSolverIdentity()
	{
	#if USE_HEAP
		return "$RCSfile: sdfflt.c,v $ $Revision: 1.45 $ $Date: 2000/12/15 06:39:46 $ (using heap)";
	#else
		#if PRIORITIZE_FALSE_CLAUSE_VARS
			return "$RCSfile: sdfflt.c,v $ $Revision: 1.45 $ $Date: 2000/12/15 06:39:46 $ (prioritized false clauses)";
		#else
			return "$RCSfile: sdfflt.c,v $ $Revision: 1.45 $ $Date: 2000/12/15 06:39:46 $";
		#endif
	#endif
	}

#if USE_HEAP

inline void heapify(REAL *valArr, int *indArr, int *invArr, int heapsize, int node)
	{
	int left, right, largest;

	while (1)
		{
		left = 2 * node;
		right = left + 1;

		if (left <= heapsize && valArr[left] > valArr[node])
			largest = left;
		else
			largest = node;
	
		if (right <= heapsize && valArr[right] > valArr[largest])
			largest = right;

		if (largest != node)
			{
			swapint(invArr[indArr[node]], invArr[indArr[largest]]);
			swap(valArr[node], valArr[largest]);
			swapint(indArr[node], indArr[largest]);
			node = largest;
			}
		else
			break;
		}
	}

inline void rebuildHeap(REAL *valArr, int *indArr, int *invArr, int heapsize, int start)
	{
	int node;

	for (node = start; node >= 1; node--)
		{
		heapify(valArr, indArr, invArr, heapsize, node);
		}

	/* checkHeap(valArr, indArr, invArr, heapsize); */
	}

inline void buildHeap(REAL *valArr, int *indArr, int *invArr, int heapsize)
	{
	rebuildHeap(valArr, indArr, invArr, heapsize, heapsize / 2);
	}

#endif