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

/* Subgradient optimizer using multiplicative weight updates.                */

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

#define SHOW_DETAILS 0

typedef long double REAL;

REAL alpha, oneminusrho;
int R, T, W, COR;
REAL rho, cut;
int shift, S;
int Amax, Asum;
REAL gap;
int startempty;

REAL *Y, *AA, *penalty, *mult, *diffA, *diffL;
REAL *diffPup, *diffPdn, *diffP2up, *diffP2dn; 
int *V, *minxs, *num;


void ALLOCATE(double newalpha, double newrho, int newtries, int newflips, double newnoise, int newcorr, int newempty)
{
   int minV, maxV, minval, maxval, c, s, x;
   
   alpha = (REAL)newalpha;
   rho = (REAL)newrho;
   R = newtries;
   T = newflips;
   W = (int)(newnoise * RAND_MAX);
   COR = newcorr;
	 startempty = newempty;
   oneminusrho = (REAL)(1.0 - rho);
   /* cut = (REAL) ( log( (m + 1)/2.0 ) / log((double) m) ); */
   /*cut = 0.0;*/
   cut = 0.5;
   
   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;

   Y = (REAL *) malloc(m * sizeof(REAL));
   V = (int *) malloc(m * sizeof(int));
   AA = (REAL *) malloc(n * sizeof(REAL));
   penalty = (REAL *) malloc(S * sizeof(REAL));
   diffA = (REAL *) malloc(n * sizeof(REAL));
   diffL = (REAL *) malloc(n * sizeof(REAL));
   mult = (REAL *) malloc(S * sizeof(REAL));
   diffPup = (REAL *) malloc(S * sizeof(REAL));
   diffPdn = (REAL *) malloc(S * sizeof(REAL));
   diffP2up = (REAL *) malloc(S * sizeof(REAL));
   diffP2dn = (REAL *) malloc(S * sizeof(REAL));
   minxs = (int *) malloc(n * sizeof(int));
   num = (int *) malloc(S * sizeof(int));

   /*** Penalty ***/


   /* Step (0,1) (good only for SAT) */
   /*for (s = 0; s < S; s++) penalty[s] = (s <= shift) ? 0.0 : 1.0;*/

   /* Step (-.5,+.5) */
   /*for (s = 0; s < S; s++) penalty[s] = (s <= shift) ? -0.5 : 0.5;*/

   /* Exponential */
   /*for (s = 0; s < S; s++) penalty[s] = pow(m,s-shift-cut) - 1.0;
   norm = pow(m-1, 1-cut);
   for (s = 0; s < S; s++) penalty[s] /= norm;*/

#ifdef USE_LINEAR
   /* Linear (very poor for SAT) */
   for (s = 0; s < S; s++) penalty[s] = (REAL) (s - shift);
#else
   /* Hinge */
/*   for (s = 0; s <= shift; s++) penalty[s] = -0.5;
   for (s = shift+1; s < S; s++) penalty[s] = s-shift-0.5; */
   
	 for (s = 0; s <= shift; s++) penalty[s] = -cut;
/*   for (s = shift+1; s < S; s++) penalty[s] = 1.95*(s-shift) - cut; */
   for (s = shift+1; s < S; s++) penalty[s] = s-shift - cut;
#endif
   

   /*** Scaled objective ***/
   gap = penalty[shift+1] - penalty[shift];
   Amax = 0;
   Asum = 0.0;
   for (x = 0; x < n; x++) {
      AA[x] = (REAL) A[x];
      if (abs(A[x]) > Amax) Amax = abs(A[x]);
      Asum += (REAL) A[x];
   }
   for (x = 0; x < n; x++) AA[x] *= gap / Amax / m / 2.0;

   /*** Pre-calc tables ***/
   for (s = 0; s < S; s++) {
      mult[s] = (REAL)pow((double)alpha, (double)penalty[s]);
      if (s < S-1) diffPup[s] = penalty[s+1] - penalty[s];
      if (s > 0) diffPdn[s] = penalty[s-1] - penalty[s];
      if (s > 0 && s < S-1) {
         diffP2up[s] = 2*penalty[s] - penalty[s+1] - penalty[s-1];
         diffP2dn[s] = -diffP2up[s];
      }
   }
}


void SEARCH(double newalpha, double newrho, int newtries, int newflips, double newnoise, int newcorr, int newempty)
{ 
   int x, c, s, i, j, d, nmxs, xflip, numvio;
   int curobj, Scurr, Sflip, agreexflip;
   REAL curval, minfeasval, mindiff, ytotal, oneminusrhomeanY;
   int *Xsc, *Xcoefsc;
   
   minfeasval = DBL_MAX;
   minfeasobj = INT_MAX;
   steps = considered = dsteps = 0;
   
   
   for (r = 0; r < R; r++) {  /*******************************************/
      if (startempty)
	       for (x = 0; x < n; x++) X[x] = -1;
			else
         for (x = 0; x < n; x++) X[x] = mrand48() < 0 ? -1 : 1;

			considered++;

      numvio = 0;
      for (s = 0; s < S; s++) num[s] = 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) numvio++;
         V[c] /= 2;
         Scurr = V[c] + shift;
         num[Scurr]++;
         Y[c] = 1.0 / m;
      }
      d = 1;

      curval = 0.0;
      curobj = 0;
      for (x = 0; x < n; x++) {
         if (X[x] > 0) {
            diffA[x] = -2*AA[x];
            curval += AA[x];
            curobj += A[x];
         } else {
            diffA[x] = 2*AA[x];
            curval -= AA[x];
         }
         diffL[x] = diffA[x];
         for (j = 0; j < nCs[x]; j++) {
            c = Cs[x][j];
            if (X[x] == Ccoefs[x][j])
               diffL[x] += Y[c] * diffPdn[ V[c]+shift ];
            else
               diffL[x] += Y[c] * diffPup[ V[c]+shift ];
         }
      }

      if (!numvio) {
         minfeasval = curval;
         minfeasobj = curobj;
      }

      while (steps < T) {              /***************************/
   
         if (minfeasobj <= opt) break;

         if (d % COR == 0) {                        /*** Correction ****/
            for (x = 0; x < n; x++) {
               diffL[x] = diffA[x];
               for (j = 0; j < nCs[x]; j++) {
                  c = Cs[x][j];
                  if (X[x] == Ccoefs[x][j])
                     diffL[x] += Y[c] * diffPdn[ V[c]+shift ];
                  else
                     diffL[x] += Y[c] * diffPup[ V[c]+shift ];
               }
            }
         }
         d++;

         mindiff = (REAL)DBL_MAX;                 /*** Candidates ****/
         nmxs = 0;
         for (x = 0; x < n; x++) {
            if (diffL[x] <= mindiff) {
               if (diffL[x] < mindiff) {
                  mindiff = diffL[x];
                  minxs[0] = x;
                  nmxs = 1;
               } else
                  minxs[nmxs++] = x;
            }
         }
         if (mindiff >= 0.0 && random() < W) {
            minxs[0] = random()%n;
            nmxs = 1;
            mindiff = -1.0;
         }

         if (mindiff >= 0.0) {                      /*** Dual step *****/

            dsteps++;
            d = 0;

            ytotal = 0.0;
            for (c = 0; c < m; c++) {
               Y[c] *= mult[ V[c] + shift ];
               ytotal += Y[c];
            }
            oneminusrhomeanY = oneminusrho * ytotal / m;
            for (c = 0; c < m; c++) 
               Y[c] = rho * Y[c] + oneminusrhomeanY;

         } else {                                   /*** Primal step ***/

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

            X[xflip] = -X[xflip];
            diffL[xflip] = -diffL[xflip];
            diffA[xflip] = -diffA[xflip];

            if (X[xflip] > 0) {
               curval += 2*AA[xflip];
               curobj += A[xflip];
            } else {
               curval -= 2*AA[xflip];
               curobj -= A[xflip];
            }
            
            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);
               Xsc = Xs[c];
               Xcoefsc = Xcoefs[c];

               num[Scurr]--;
               num[Sflip]++;

               if (Sflip > shift && Scurr <= shift) 
                  numvio++;
               else if (Sflip <= shift && Scurr > shift)
                  numvio--;

               for (i = 0; i < nXs[c]; i++)
                  if ((x = Xsc[i]) != xflip) {
                     if (X[x] == Xcoefsc[i])
                        diffL[x] += Y[c] * (agreexflip ? diffP2up[Scurr] 
                                                       : diffP2dn[Sflip]);
                     else
                        diffL[x] += Y[c] * (agreexflip ? diffP2dn[Sflip] 
                                                       : diffP2up[Scurr]);
                  };
               V[c] = Sflip - shift;
            }

            if (!numvio && curval < minfeasval) {
               minfeasval = curval;
               minfeasobj = curobj;

/*
printf("r%d t%d mv%d: ", r, t, minfeasobj);
for (x = 0; x < n; x++) if (X[x] > 0) printf("x%d ", x);
printf("\n");
*/

            }
         }                                          /*******************/
      }                                      /***************************/
      if (minfeasobj <= opt) break;
   }                          /*******************************************/
	 
	 if (minfeasobj == INT_MAX)
			minfeasobj = 0;
}



void CLEANUP()
{
   free(Y);
   free(AA);
   free(diffA);
   free(diffL);
   free(penalty);
   free(mult);
   free(diffPup);
   free(diffPdn);
   free(diffP2up);
   free(diffP2dn);
   free(V);
   free(minxs);
   free(num);
}