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

/* Additive updates, using our standard hinge penalty */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include "cop.h"
#include "satfront.h"

#define swapint(x,y) {tempint=x; x=y; y=tempint;}
#define swapreal(x,y) {tempreal=x; x=y; y=tempreal;}

typedef double REAL;

int tempint;
REAL tempreal;
int n, m, N, M;
int *X, *A, *nCs, **Cs, **Ccoefs;
int *B, *nXs, **Xs, **Xcoefs;
int MaxCoef, MinB, MaxB;
int steps, dsteps, considered, r;
int T, R, W;


int shift, S, COR;
REAL alpha, rho, oneminusrho;

int *incrsat, *decrsat, *incrvio, *decrvio;
int *xvreducers, *minxs, *wherexvreduce, *cviol, *wherecviol, *V;
int **numup, **numdn;
REAL *Y, *primpenalty, *dualpenalty;
REAL **Yup, **Ydn, **diff1ppen, **diff2ppen, *diffL;



void preSolve(SATParams *params)
{ 
	int minV, maxV, minval, maxval, c, s;
	double walkprob;

	T = params->intParams[0];
	if (T == 0)
		T = 500000;

	R = params->intParams[1];
	if (R == 0)
		R = 1;

	COR = params->intParams[2];
	if (COR == 0)
		COR = 100;

	alpha = params->doubleParams[0];
	if (alpha == 0)
		alpha = 1;

	/*
		 rho = params->doubleParams[1];
		 if (rho == 0)
		 rho = 0.995;

		 oneminusrho = 1 - rho;
	 */

	walkprob = params->doubleParams[2];
	if (walkprob == 0)
		walkprob = 0.002;

	W = RAND_MAX * walkprob;

	/*alpha = 1.0 + (REAL)par1 * n / (REAL)m;*/

	if (MaxCoef > 1) {
		printf("Error: this implementation assumes {-1,0,+1} valued ");
		printf("coefficients\n");
		exit(1);
	}
	minV = INT_MAX;
	maxV = INT_MIN;
	for (c = 0; c < m; c++) {
		if ( (nXs[c] + B[c]) % 2 ) {
			printf("Error: can only handle even constraint violation values\n");
			exit(1);
		}
		minval = (-nXs[c] - B[c]) / 2;
		maxval = (nXs[c] - B[c]) / 2;
		if (minval < minV) minV = minval;
		if (maxval > maxV) maxV = maxval;
	}
	S = maxV - minV + 1;
	shift = -minV;

	incrsat = (int *) malloc(n * sizeof(int));
	decrsat = (int *) malloc(n * sizeof(int));
	incrvio = (int *) malloc(n * sizeof(int));
	decrvio = (int *) malloc(n * sizeof(int));
	xvreducers = (int *) malloc(n * sizeof(int));
	minxs = (int *) malloc(n * sizeof(int));
	wherexvreduce = (int *) malloc(n * sizeof(int));
	cviol = (int *) malloc(m * sizeof(int));
	wherecviol = (int *) malloc(m * sizeof(int));
	V = (int *) malloc(m * sizeof(int));
	Y = (REAL *) malloc(m * sizeof(REAL));
	diff1ppen = (REAL **) malloc(2 * sizeof(REAL *));
	diff2ppen = (REAL **) malloc(2 * sizeof(REAL *));
	diffL = (REAL *) malloc(n * sizeof(REAL));
	primpenalty = (REAL *) malloc(S * sizeof(REAL));
	diff1ppen[0] = (REAL *) malloc(S * sizeof(REAL));
	diff1ppen[1] = (REAL *) malloc(S * sizeof(REAL));
	diff2ppen[0] = (REAL *) malloc(S * sizeof(REAL));
	diff2ppen[1] = (REAL *) malloc(S * sizeof(REAL));
	dualpenalty = (REAL *) malloc(S * sizeof(REAL));
	numup = (int **) malloc(S * sizeof(int *));
	numdn = (int **) malloc(S * sizeof(int *));
	Yup = (REAL **) malloc(S * sizeof(REAL *));
	Ydn = (REAL **) malloc(S * sizeof(REAL *));

	for (s = 0; s < S; s++) {
		numup[s] = (int *) malloc(n * sizeof(int));
		numdn[s] = (int *) malloc(n * sizeof(int));
		Yup[s] = (REAL *) malloc(n * sizeof(REAL));
		Ydn[s] = (REAL *) malloc(n * sizeof(REAL));

		/*** Dual penalty ***/
		dualpenalty[s] = (s <= shift) ? 0.0 : 1.0;
	}

	/*** Primal penalty ***/
	for (s = 0; s < S; s++) primpenalty[s] = (s <= shift) ? 0.0 : 1.0;

	for (s = 0; s < S; s++) {
		if (s > 0)   diff1ppen[1][s] = primpenalty[s-1] - primpenalty[s];
		if (s < S-1) diff1ppen[0][s] = primpenalty[s+1] - primpenalty[s];
		if (s > 0 && s < S - 1) {
			diff2ppen[0][s] = 2*primpenalty[s]-primpenalty[s+1]-primpenalty[s-1];
			diff2ppen[1][s] = -diff2ppen[0][s];
		}
	}
}


void runSolver(SATParams *params, SATResults *results)
{
	REAL Yval = 1.0 / m;

	int x, c, s, i, j, t, d, nmxs, xflip, nxvr, numvio = 0;
	int Scurr, Sflip, agreexflip, agreex;
	REAL mindiff, ytotal;

	steps = considered = dsteps = 0;

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

		for (x = 0; x < n; x++) 
			X[x] = (mrand48() < 0) ? -1 : 1;
		steps++;
		considered++;

		numvio = 0;
		for (c = 0; c < m; c++) {
			V[c] = -B[c];
			for (i = 0; i < nXs[c]; i++) 
				V[c] += X[Xs[c][i]] == Xcoefs[c][i] ? 1 : -1;
			if (V[c] > 0) {
				cviol[numvio] = c;
				wherecviol[c] = numvio;
				numvio++;
			}
			V[c] /= 2;
			Y[c] = Yval;
			Scurr = V[c] + shift;
		}
		d = 1;

		nxvr = 0;
		for (x = 0; x < n; x++) {
			incrsat[x] = decrsat[x] = incrvio[x] = decrvio[x] = 0;
			diffL[x] = 0.0;
			for (s = 0; s < S; s++) {
				numup[s][x] = numdn[s][x] = 0;
				Yup[s][x] = Ydn[s][x] = 0.0;
			}
			for (j = 0; j < nCs[x]; j++) {
				c = Cs[x][j];
				Scurr = V[c] + shift;
				if (X[x] == Ccoefs[x][j]) {
					numdn[Scurr][x]++;
					Ydn[Scurr][x] += Y[c];
					diffL[x] += Y[c] * diff1ppen[1][Scurr];
					if (V[c] > 0) {
						if (!decrvio[x]++) {
							xvreducers[nxvr] = x;
							wherexvreduce[x] = nxvr;
							nxvr++;
						}
					} else
						decrsat[x]++;
				} else {
					numup[Scurr][x]++;
					Yup[Scurr][x] += Y[c];
					diffL[x] += Y[c] * diff1ppen[0][Scurr];
					if (V[c] > 0)
						incrvio[x]++;
					else
						incrsat[x]++;
				}
			}
		}

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

			if (numvio == 0) break;

			if (d % COR == 0) {                        /*** CORRECTION ****/


				ytotal = 0.0;
				for (c = 0; c < m; c++) ytotal += Y[c];
				for (c = 0; c < m; c++) {
					Y[c] /= ytotal;
				}
				d = 1;

				for (x = 0; x < n; x++) {
					diffL[x] = 0.0;
					for (s = 0; s < S; s++) Yup[s][x] = Ydn[s][x] = 0.0;
					for (j = 0; j < nCs[x]; j++) {
						c = Cs[x][j];
						Scurr = V[c] + shift;
						if (X[x]==Ccoefs[x][j]) {
							Ydn[Scurr][x] += Y[c];
							diffL[x] += Y[c] * diff1ppen[1][Scurr];
						} else {
							Yup[Scurr][x] += Y[c];
							diffL[x] += Y[c] * diff1ppen[0][Scurr];
						}
					}
				}
			}

			mindiff = DBL_MAX;                       /*** Candidates ****/
			nmxs = 0;
			for (i = 0; i < nxvr; i++) {
				x = xvreducers[i];
				if (diffL[x] < mindiff) {
					mindiff = diffL[x];
					minxs[nmxs = 0] = x;
					nmxs++;
				} else if (diffL[x] == mindiff)
					minxs[nmxs++] = x;
			}

			if (mindiff >= 0.0 && random() < W) {
				nmxs = 1;
				x = minxs[0] = lrand48()%n;
				diffL[x] = 0.0;
				if (decrvio[x] == 0)
					for (s = 0; s < S; s++) diffL[x] += Yup[s][x] * diff1ppen[0][s]
						+ Ydn[s][x] * diff1ppen[1][s];
				mindiff = -1.0;
			}

			if (mindiff >= 0.0) {                      /*** DUAL update ***/

				dsteps++;

				for (c = 0; c < m; c++) Y[c] += alpha * dualpenalty[V[c]+shift];

				for (i = 0; i < nxvr; i++) diffL[xvreducers[i]] = 0.0;
				for (s = 0; s < S; s++) {
					for (x = 0; x < n; x++) {
						Yup[s][x] += numup[s][x] * alpha * dualpenalty[s];
						Ydn[s][x] += numdn[s][x] * alpha * dualpenalty[s];
						if (decrvio[x] > 0) diffL[x] += Yup[s][x] * diff1ppen[0][s] 
							+ Ydn[s][x] * diff1ppen[1][s];
					}
				}
				d++;

			} else {                                   /*** PRIMAL update */

				if (nmxs == 1)
					xflip = minxs[0];
				else
					xflip = minxs[lrand48()%nmxs];

				X[xflip] = -X[xflip];
				diffL[xflip] = -diffL[xflip];
				for(s = 0; s < S-1; s++) {
					swapint(numup[s][xflip], numdn[s+1][xflip]);
					swapreal(Yup[s][xflip], Ydn[s+1][xflip]);
				}
				swapint(incrsat[xflip],decrsat[xflip]);
				swapint(incrvio[xflip],decrvio[xflip]);
				if (incrvio[xflip] > 0 && decrvio[xflip] == 0) {
					nxvr--;
					xvreducers[wherexvreduce[xflip]] = xvreducers[nxvr];
					wherexvreduce[xvreducers[nxvr]] = wherexvreduce[xflip];
				} else if (incrvio[xflip] == 0 && decrvio[xflip] > 0) {
					xvreducers[nxvr] = xflip;
					wherexvreduce[xflip] = nxvr;
					nxvr++;
				}
				steps++;
				considered += n;

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

					c = Cs[xflip][j];
					agreexflip = X[xflip] == Ccoefs[xflip][j];
					Scurr = V[c] + shift;
					Sflip = Scurr + (agreexflip ? 1 : -1);

					if (Sflip > shift && Scurr <= shift) {
						cviol[numvio] = c;
						wherecviol[c] = numvio;
						numvio++;
						for (i = 0; i < nXs[c]; i++) {
							x = Xs[c][i];
							if (X[x] == Xcoefs[c][i]) {
								if (!decrvio[x]++) {
									xvreducers[nxvr] = x;
									wherexvreduce[x] = nxvr;
									nxvr++;
									diffL[x] = 0.0;
									for (s = 0; s < S; s++)
										diffL[x] += Yup[s][x] * diff1ppen[0][s]
											+ Ydn[s][x] * diff1ppen[1][s];
								}
								decrsat[x]--;
							} else {
								incrvio[x]++;
								incrsat[x]--;
							}
						}

					} else if (Sflip <= shift && Scurr > shift) {
						numvio--;
						cviol[wherecviol[c]] = cviol[numvio];
						wherecviol[cviol[numvio]] = wherecviol[c];
						for (i = 0; i < nXs[c]; i++) {
							x = Xs[c][i];
							if (X[x] == Xcoefs[c][i]) {
								if (!--decrvio[x]) {
									nxvr--;
									xvreducers[wherexvreduce[x]] = xvreducers[nxvr];
									wherexvreduce[xvreducers[nxvr]] = wherexvreduce[x];
								}
								decrsat[x]++;
							} else {
								incrvio[x]--;
								incrsat[x]++;
							}
						}
					}

					for (i = 0; i < nXs[c]; i++)
						if ((x = Xs[c][i]) != xflip) {
							if ((agreex = X[x]) == Xcoefs[c][i]) {
								numdn[Scurr][x]--;
								numdn[Sflip][x]++;
								Ydn[Scurr][x] -= Y[c];
								Ydn[Sflip][x] += Y[c];
							} else {
								numup[Scurr][x]--;
								numup[Sflip][x]++;
								Yup[Scurr][x] -= Y[c];
								Yup[Sflip][x] += Y[c];
							}
							if (decrvio[x] > 0)
								diffL[x] += Y[c] * ( agreex == agreexflip ?
										diff2ppen[0][Scurr] :
										diff2ppen[1][Sflip]);
						};

					V[c] = Sflip - shift;
				}
			}                                          /*******************/
		}                                      /***************************/
		if (numvio == 0) break;
	}                          /*******************************************/

	if (numvio > 0) 
		r--;

	results->numFlips = steps;
	results->numReweights = dsteps;
	results->numRestarts = r + 1;
	results->numSatisfied = m - numvio;
	results->isSatisfied = numvio == 0;
	#include "copysoln.c"
}


void postSolve(SATParams *params)
{
	int s;

	free(Y);
	free(V);
	free(primpenalty);
	free(dualpenalty);
	free(diffL);
	free(minxs);
	free(incrsat);
	free(decrsat);
	free(incrvio);
	free(decrvio);
	free(xvreducers);
	free(wherexvreduce);
	free(cviol);
	free(wherecviol);
	free(diff1ppen[0]);
	free(diff1ppen[1]);
	free(diff1ppen);
	free(diff2ppen[0]);
	free(diff2ppen[1]);
	free(diff2ppen);
	for (s = 0; s < S; s++) {
		free(numup[s]);
		free(numdn[s]);
		free(Yup[s]);
		free(Ydn[s]);
	}
	free(numup);
	free(numdn);
	free(Yup);
	free(Ydn);

}
SATParams *getSATParams()
{
	SATParams *params;

	params = malloc(sizeof(SATParams));
	params->numIntParams = 3;
	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] = "-cp";
	params->intParamDescs[2] = "number of reweights between corrections";
	params->intParams[2] = 0;


	params->numDoubleParams = 3;
	params->doubleParams = malloc(params->numDoubleParams * sizeof(double));
	params->doubleParamNames = malloc(params->numDoubleParams * sizeof(char*));
	params->doubleParamDescs = malloc(params->numDoubleParams * sizeof(char*));
	params->doubleParamNames[0] = "-alpha";
	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->doubleParamNames[2] = "-noise";
	params->doubleParamDescs[2] = "probability of random walk when stuck";
	params->doubleParams[2] = 0;

	params->numStringParams = 0;

	return params;
}


char *getSolverIdentity()
{
	return "$RCSfile: asg.c,v $ $Revision: 1.8 $ $Date: 2001/11/23 06:35:04 $";
}